Elastic neutron scattering models for NCrystal

The NCrystal library provides a range of models for simulation of both elastic and inelastic scattering of thermal neutrons in a range of material structures. This article presents the available models for elastic scattering, and includes detailed discussion of their theoretical background, their implementation, and in particular their validation. The lineup includes a model for Bragg diffraction in crystal powders as well as one for incoherent elastic scattering, but the main focus is given to models of Bragg diffraction in ideally imperfect single crystals: both for the most widely applicable model of isotropic Gaussian mosaicity, and for a more specific model of layered single crystals which is relevant for materials such as pyrolytic graphite. Although these single crystal models are utilising computationally efficient approximations where appropriate, attention is given to the provision of precise and trustworthy results also for the extreme cases of back-scattering, forward-scattering, and crystals with very large mosaic spreads. Together with NCrystal's other features for crystal structure initialisation and inelastic physics, the presented models enable realistic modelling of components at neutron scattering instruments in frameworks like Geant4 and McStas, including monochromators, analysers, filters, support materials, shielding, and many kinds of samples. As a byproduct of the work, an improved formula for approximating cross sections in isotropic single crystals with Gaussian mosaicity is provided.

in crystal powders as well as one for incoherent elastic scattering, but the main focus is given to models of Bragg diffraction in ideally imperfect single crystals: both for the most widely applicable model of isotropic Gaussian mosaicity, and for a more specific model of layered single crystals which is relevant for materials such as pyrolytic graphite. Although these single crystal models are utilising computationally efficient approximations where appropriate, attention is given to the provision of precise and trustworthy results also for the extreme cases of back-scattering, forward-scattering, and crystals with very large mosaic spreads.
Together with NCrystal's other features for crystal structure initialisation and inelastic physics, the presented models enable realistic modelling of components at neutron scattering instruments in frameworks like Geant4 and McStas, including monochromators, analysers, filters, support materials, shielding, and many kinds of samples. As a byproduct of the work, an improved formula for approximating cross sections in isotropic single crystals with Gaussian mosaicity is provided.

Introduction
The software package NCrystal [1] is an open source software package, with capabilities for modelling of thermal neutron transport in a variety of materials.
The code is freely available [2] and released under a liberal open source license [3].
The initial release, v1.0.0 [1], focused on Bragg diffraction in crystals, while release v2.0.0 [4] introduced realistic modelling of both inelastic and incoherent elastic scattering processes, as well as the support for certain non-crystalline materials like liquids. The following releases up until v2.2.1 [5], were mostly focused on technical changes not affecting the physics modelling (the notable exception being that atomic data definitions were made more flexible, allowing for atoms not simply being natural elements).
The publication [1] provided an in-depth overview of the NCrystal framework, including issues of configuration and material definition, and how users might interact with the library: through standalone tools, interfaces for C++, C and Python programming languages, or indirectly by using NCrystal to enhance the physics capabilities of simulation packages such as Geant4 [6,7,8] and McStas [9,10]. It also provided a review of the relevant neutron scattering theory, and focused in particular on the initialisation of crystal structures and form factors, and included benchmarks versus existing crystallographic software and powder diffraction data. It did not, however, go into detailed discussions of the actual implementation of particular neutron scattering models. The present article improves upon this situation by providing a thorough review of the models for elastic scattering presently available in NCrystal, including details of implementation, performance, and validation.
Although not discussed further in the present article, inelastic scattering does contribute to a few of the validation plots presented, and it is therefore worthwhile to briefly note that the improved support for inelastic scattering introduced in NCrystal v2.0.0 is based on sampling of scattering kernels [11], which can be either provided by the users or calculated on-the-fly by NCrystal from phonon density of state (DOS) curves. This latter calculation relies on an independent implementation of a method due to A. Sjölander [12], which is also used for a similar purpose in well established applications like LEAPR [13]. The combination of these new inelastic models and the elastic models described in the present article, enables highly realistic simulations of typical components at neutron scattering instruments, including beam filters, monochromators, analysers, detectors, support materials and many kinds of samples. The realism is increased to unprecedented levels when NCrystal is used as a backend in an application Geant4, where the injection of realistic thermal neutron scattering models nicely complements the pre-existing features for support of complicated geometries and both nuclear-and electromagnetic physics.
It should be noted that while inelastic models improved greatly between After establishing some common formalism and terminology, and reviewing common features of Bragg diffraction in Section 2, the next three sections are devoted to Bragg diffraction under particular models of crystallite distributions: powders (Section 3), single crystals with isotropic Gaussian mosaicity (Section 4), and single crystals with a layered structure similar to pyrolytic graphite (Section 5). Along the way, issues such as theory, details of implementation, and validation are considered. A noteworthy byproduct of the work on single crystals is the derivation of a formula in Section 4.4 (Eq. 35) which can be used to approximate scattering cross sections. This formula is an improvement over the widely used Gaussian approximation, in that it more accurately incorporates effects related to large mosaic spreads or configurations approaching back-scattering. Next, incoherent elastic scattering is discussed in Section 6.
This section also considers the validity of the commonly assumed isotropic angular distribution of this process. Finally, the computational efficiency of all the presented models is discussed in Section 7, and possible future developments are considered in Section 8.

Theoretical background and features of Bragg diffraction
The theory of scattering by thermal neutrons in crystals is reviewed in [1,Sec. 2]. There, it is described how neutron scattering in crystals can be treated under the Born approximation, under certain general conditions. Notably this includes the assumption that considered crystal systems are "ideally imperfect" in the sense that the crystal structure can essentially be assumed to be perfect inside tiny crystal grains, or "crystallites", but that interference between interactions in separate grains can be neglected due to their imperfect mutual alignments. Formulas are also developed under the assumption that thermal movements of individual atoms around their nominal positions in the lattice are essentially harmonic and isotropic. Based on this, thermal fluctuation effects are incorporated via Debye-Waller factors, that are themselves estimated based on mean-squared-displacements estimated via the Debye model and a phenomenological parameter, the Debye temperature. Finally, it is discussed in [1,Sec. 4] how NCrystal is able to derive lists of all relevant reflection planes in a given crystal, along with their associated parameters like form factors, d-spacings, and plane normals. For more details, readers are referred to [1] and references therein, in particular [14,15,16]. The contribution of a given lattice plane hkl to the differential coherent elastic cross section for scattering the incident neutron with wavevector k i and energy E i into the final state with wavevector k f and energy E f , is given by [1,Sec. 2.3]: Which is normalised to the number of atoms in the target. Here, V uc is the unit cell volume, n a the number of atoms per unit cell, and Q ≡ k f − k i the momentum transfer. The factor of δ(∆E) ≡ δ(E f − E i ) imposes energy con- Where λ = 2π/k i is the neutron wavelength. When the condition is fulfilled, the scattering angle θ will satisfy the Bragg equation: λ = 2d hkl sin(θ/2) = 2d hkl sin θ B Where also the Bragg angle, θ B ≡ θ/2, has been introduced. For the purposes of the present article, the complementary angle α is additionally defined as: With this definition, α is given as the angle between − k i andn when the alignment satisfies Eq. 3. Thus, Bragg diffraction can take place only when the plane normaln is located on a cone around − k i of opening angle α. The intersection of this cone with the unit sphere defines a circle of radius sin α, which in the present article will be referred to as the Bragg circle.
As the crystalline systems considered here are at the macroscopic scale composed of a number of independently oriented microscopic crystallites, one must average Eq. 1 over the distribution of microscopic crystallite orientations in order to derive macroscopic cross sections. Denoting such a mosaicity distribution of crystallite orientations with W , the total macroscopic cross sections due to a particular hkl plane is thus given by an integral over neutron final states and crystallite orientations: Where Ωn represents the solid angles covered by the direction of the normal n ≡ τ hkl /τ hkl in specific crystallites, and W (n) is the density of crystallites expressed as a function ofn. Denoting the cosine of the angle betweenn and − k i with µ and a corresponding azimuthal angle with t, W (n) can be written as W (µ, t) and Eq. 5 becomes: If Eq. 2 is not satisfied, this trivially evaluates to zero, as Q can then never equal τ hkl . Otherwise it can be evaluated using the identity dΩδ( r − a) = 2a −1 δ(r 2 − a 2 ): Where it was used that λ = 2π/k i , cos α = λ/2d hkl = τ hkl /2k i (since τ hkl = 2π/d hkl ), and sin 2α = 2 sin α cos α. Thus, the resulting cross section can be factorised into two parts: Here, q hkl represents an intrinsic strength of the interaction: And g hkl is a geometrical factor depending on the Bragg angle, mosaicity distribution, and direction of the incoming neutron. It is given by an integration of the crystallite densities over the directions where Bragg diffraction is possible (i.e. along the Bragg circle which has a radius of sin α): In addition to providing cross section values, implementations of Bragg diffraction models in NCrystal must also provide Monte Carlo-based sampling of the direction of k f in case of a scattering event. When multiple reflection planes contribute, one is first trivially selected at random, with a probability given by its relative contribution to the total cross section. Next, scattering on the chosen hkl plane can proceed using whatever method is best suited to the specific mosaicity distribution. One such method proceeds as follows: first an actual (as opposed to nominal) direction of the plane normal,n, must be sampled according to the contribution to the cross section. This corresponds to selecting a specific orientation of the crystallite in which the neutron scatters.
Subsequently, the neutron must undergo a specular reflection on the plane defined by that normal, resulting in k f = k i + τ hkln . The plane normal sampling can be carried out by sampling a value of t according to the contribution to the integral in Eq. 10, or in other words P (t) ∝ W (µ = cos α, t).

Crystal powders
The simplest model for Bragg diffraction in NCrystal implements the socalled powder approximation, in which individual crystallite orientations are assumed to be completely independent and uniformly distributed over all solid angles. This approximation is not only suitable for modelling actual crystal powders, but can additionally 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. Due to its simplicity and usefulness, implementations of Bragg diffraction under the powder assumption are widely available (for instance [17,18,19,20,21]). The implementation in NCrystal has a particular focus on efficiency, which is necessary since the aim is to ensure realistic cross sections also at shorter wavelengths, potentially requiring very long lists of reflection planes (cf. [1,Sec. 4]).
The mosaicity distribution in the powder approximation is a constant, W = 1/4π, which when inserted in Eq. 10 yields: and thus: Unless Eq. 2 is not satisfied, in which case σ hkl (λ) vanishes. The complete cross section for Bragg diffraction in a crystal powder can thus be written as the following sum over all lattice plane families satisfying Eq. 2: Although it is straightforward to evaluate Eq. 13, a naive implementation using it with a list of N hkl lattice planes to provide the cross section for a particular neutron, would in the worst scenarios result in evaluation times scaling as O(N hkl ). As shown in [1,Fig. 4], N hkl might be higher than 10 6 when includ-   used in the second array to immediately evaluate the sum in Eq. 13. Additionally, the arrays are kept as short as possible by merging entries with identical d-spacing values. Figure 1 shows a typical example of Bragg diffraction cross section in a crystal powder, as provided by NCrystal.
In case of a scattering event, the relative probability for it to involve a plane with a particular d-spacing value will depend on the total contribution at that d-spacing to the sum in Eq. 13. The Monte Carlo-based selection of such a d-spacing value is carried out by another binary search, this time in the array with cumulative sums of d|F | 2 . The search target is the cumulative value found at the index i but multiplied with a pseudo random number from the unit interval, thus ensuring the desired weight when selecting different d-spacings.
Specifically, the index j resulting from this second binary search can be used to access a value in the array of h 2 /8md 2 values, which can be used to extract the cosine of the scattering angle: It trivially follows from the fact that W = 1/4π is a constant, that the azimuthal scattering angle must be uniformly distributed between −π and π. Now, if the calling code is querying the non-oriented NCrystal interface method (cf. [1, Table 1]), the value of θ will be provided to the user after an arccos evaluation. If instead (as will be most typical), the calling code uses the oriented vector interface, providingk i and expecting a value ofk f in return, an alternative procedure is used. This procedure does not require any expensive trigonometric function calls, and is as fast as the non-oriented one, despite the fact that it must sample the azimuthal scattering angle and deal with full directional vectors. First, a vector is sampled isotropically on the unit sphere using an efficient algorithm [23], and the sampling is redone in the rare case of providing a vector almost co-linear withk i . Next, a cross product between the sampled vector and k i yields a vector v orthogonal to k i but with uniformly distributed azimuthal scattering angle. Finally,k f is provided as: Figure 2 shows an example of the sharply defined scattering angles in a crystal powder, as provided by NCrystal.  ability of NCrystal to initialise crystal structures in [1,Sec. 4.4], to which interested readers are referred.

Single crystals with isotropic Gaussian mosaicity
As is often the case in physics, a Gaussian distribution lies at the heart of the most widely used model for the distributions of crystallites in single crystal materials utilised for monochromators, analysers, or filters at neutron scattering facilities. In this model, introduced by C. G. Darwin [24], the deviation of crystallite orientations from some reference is described by a Gaussian distribution of the angular deviation. The term mosaic spread, or simply mosaicity, is commonly used to indicate either the standard deviation or FWHM value of this angular distribution, values of which in typical materials range from tens of arc seconds to several degrees. Such a Gaussian mosaicity distribution is considered isotropic in the sense that the resulting distribution of crystallite orientations is independent of the choice of reference orientation -a feature which is for instance not the case for the distribution discussed in Section 5. Obviously it is just the mosaicity distribution which is considered isotropic in this sense, as derived properties like cross sections will show a very strong dependency on the direction of the neutron.
The required integrations of angular Gaussian mosaicity distributions in Eq. 10 are non-trivial and cumbersome at best, and as a consequence software [19,20,25]  2 The authors of [25] even goes so far as mistakenly concluding a breakdown of the theory in the presence of back-scattering due to the factor of 1/ sin 2θ B in Eq. 9. However, the divergence from this factor for θ B → π/2 is actually cancelled by the factor of sin α = sin(π/2 − θ B ) in Eq. 10 (the divergence for θ B → 0 is cancelled by the factor of λ 3 in Eq. 9).

Definition of the Gaussian mosaicity distribution
It is straightforward to simply define a Gaussian mosaicity distribution as: Where the angle δ signifies the deviation from the nominal reference orientation, and σ the mosaicity as a standard deviation value (the FWHM mosaicity value is then approximately 2.3548σ). However, the distribution given in Eq. 16 yields non-zero densities for all δ ∈ [0, π], and therefore implies a finite scattering cross section for any direction of the incoming neutron. While faithful to the concept of a Gaussian having infinite tails, such an idealised definition implies that the calculation of scattering cross sections for a particular incident neutron always has to involve the estimation of partial cross sections for scattering on all lattice planes satisfying Eq. 2. In particular at shorter neutron wavelengths where many reflection planes must be considered, this can imply very long evaluation times -with most effort being spent integrating through truly negligible (and thus completely uninteresting) parts of the Gaussian tails.
Instead, it is desirable to cleanly truncate the tails at some truncation angle, τ , which when defined at some reasonable value such as τ ≡ 5σ, allows for very significant algorithmic speedups with minimal degradation of realism or precision. Of course, the actual value of τ /σ could be a configurable parameter, depending on the precision-to-performance trade-off needed for a particular use case. However, for consistency and in order to keep configuration as simple as possible, only a single overall precision parameter, , is exposed to users, via the configuration parameter named mosprec (cf.  instead of numerical integration, as will be discussed further in Section 4.4. For the case of τ /σ, it is adjusted so that a Gaussian in a two-dimensional plane will have 1 − of its volume inside a radius of τ . Adding also a safety factor of 10% and requiring τ /σ to be at least 3, the relationship becomes: Which is shown in Figure 3. For the default value of = 10 −3 , this implies a value of τ ≈ 4.0886σ. Conversely, a value of τ ≈ 5σ is achieved by setting = 3.262 × 10 −5 . In order to simplify later calculations, it is explicitly required that τ < π/2, which guarantees that cos τ is positive. At = 10 −7 this implies that σ can at most be 14.4 • , corresponding to a FWHM mosaicity value of 33.9 • .
Fortunately, this value is far above even the largest mosaic spreads encountered in practice.
Having thus determined both the mosaicity σ and the truncation angle τ , the Gaussian mosaicity distribution will in NCrystal be defined as: Where the angle δ as before signifies the deviation from the nominal reference orientation, Θ is the Heaviside step function enforcing the truncation, and N is a normalisation factor which can be determined by the condition: During initialisation of a particular single crystal model, NCrystal determines N by numerically evaluating this integral to a relative error of less than O(10 −12 ).

Formulating the geometrical integral
The value of the integral in Eq. 10 with the Gaussian mosaicity distribution in Eq. 18 only depends on α, and the relative angle between the incoming neutron and the nominal normal of the lattice plane in question,n ≡ τ hkl /τ hkl .
For convenience this latter angle, γ, will be defined as the angle betweenn and − k i and a coordinate system specific to each neutron direction and hkl normal will be adopted, as illustrated in Figure 4: − k i lies along the positive z-axis, and the nominal normal is located in the xz-plane at: Bragg diffraction is only possible when plane normals are located at an angle of α from − k i , defining a circle as indicated in red in Figure 4, which can be parameterised as: As mentioned previously, this circle is in the present article referred to as the Bragg circle. With δ(t) indicating the angular separation between the nominal normal position (n) andp(t), Eq. 10 becomes: Taking the dot product ofn andp(t) yields the relation: cos δ(t) = sin α sin γ cos t + cos α cos γ (23) is indicated with blue shaded areas. The angle between −k i and plane normals satisfying the condition for Bragg diffraction is given as α = π/2 − θ B , defining the Bragg circle (red). The remaining variables,p(t), t, t , and δ(t), concern a parameterisation of the Bragg circle, which is discussed in the main text.
As sin α and sin γ are both non-negative, it is straightforward to confirm (with trigonometric addition formulas) what one might also intuitively infer from Figure 4: that the minimal value of δ(t) is always attained at t = 0 and is δ 0 ≡ δ(0) = |α − γ|. Likewise, a maximal value of |α + γ| is attained at t = ±π. Now, with the chosen definitions, δ(−t) = δ(t), so: where t is the point in [0, π] satisfying δ(t ) = τ , unless the Bragg circle is fully contained within the region of non-zero mosaic density (i.e. |α + γ| < τ ) in which case t = π. Of course, if the Bragg circle does not intersect the region of non-zero mosaic density at all (i.e. |α − γ| > τ ), g hkl is trivially vanishing.

Pre-selecting contributing lattice planes
Before proceeding to evaluate Eq. 24, it is important to note that in many typical scenarios, the vast majority of lattice planes fulfilling the Bragg condition 2d hkl < λ, will fail the geometrical requirement |α − γ| < τ , and therefore not contribute to the scattering at all. To avoid wasting resources on evaluating Eq. 24 for such planes, it is therefore crucial to be able to discard such noncontributing lattice planes with minimal effort. In practice this pre-filtering is often where most CPU resources are spent during a simulation of thermal neutrons in single crystals, with even small optimisations here benefiting overall simulation time.
In the implementation described here, lattice plane data is first of all grouped into "families" of planes having identical values of d-spacing and form factors, differing only in directions of their plane normals. 3 For a given incident neutron, the value of cos α only needs to be calculated once per family, and this can be done inexpensively using cos α = λ × (1/2d), which costs a mere multiplication when evaluating for a given neutron state. Additionally, lattice planes (h, k, l) and (−h, −k, −l) form natural pairs with identical form factors and d-spacings, but opposite normals (loosely speaking they can be thought of as representing the front and back side of the same plane). Thus, both planes in such a pair will reside in the same family, and for the sake of evaluating Eq. 24 they will only differ in having complementary values of γ: if one plane has a certain value of γ, the other will have the value π − γ. It is thus advantageous to only cache one of two such paired normals in memory and deal with both paired planes simultaneously. Now, the requirement for a plane to have a non-zero contribution is: As sin α and sin γ are both non-negative, Eq. 25 is always satisfied when the right hand side is negative. Thus, the condition for a neutron to have non-zero scattering cross section with a given normal can be rewritten as: Which is efficient to evaluate, since it involves only cheap non-branching operations at each side of the inequality. The only trigonometric factor particular to each pair of normals, cos γ, is readily available via a dot product between the cached normal andk i , but both + cos γ and − cos γ must be considered in order to deal with both of the paired planes being investigated jointly. However, the plane whose normal points into the same hemisphere ask i (i.e. cos γ < 0) can only ever contribute if the plane whose normal points into the hemisphere opposite tok i (i.e. cos γ > 0) contributes. Formally this follows from the fact that cos α is non-negative, so cos α cos γ ≤ cos α |cos γ|, and therefore neither plane can contribute unless: Only when Eq. 27 is valid is it possible for either of the two paired normals to contribute, in which case Eq. 26 can be used to test them individually. In this respect, note that in case of extreme forward scattering, it is entirely possible for both of the two paired planes to contribute to the scattering cross section for a particular neutron. 4

Evaluating the geometrical integral
The evaluation of Eq. 24 is either performed through direct numerical integration or, when possible, through the application of an appropriate approximation. In both cases, values of sin α and sin γ are computed relatively inexpensively using sin = √ 1 − cos 2 . In case direct numerical integration is required, the next step is the determination of the upper limit of integration, t . When sin γ sin α = 0, one can insert t = t and δ(t ) = τ in Eq. 23 and get: cos t = cos(τ ) − cos α cos γ sin γ sin α The degenerate cases where sin γ or sin α are vanishing must be handled separately, but fortunately Eq. 24 is particularly simple to evaluate in those scenarios: the right hand side vanishes when α = 0, and the integrand is independent of t when γ = 0. In the non-degenerate cases, Eq. 28 can be solved for t through a call to arccos when the right hand side leads to a value in [−1, 1]. Given the pre-selection described in Section 4.3, this only fails when the Bragg circle is fully contained within the region of non-zero mosaic density, i.e.: In either case, the correct value of t is given by: Precise numerical integration of Eq. 24 is potentially very expensive, as it requires the integrand to be evaluated repeatedly, needing calls to mathematical functions exp, arccos, and cos at each sampled value of t. The exp and arccos calls are avoided by generating, at initialisation time, a cubic spline represen- , which must then be evaluated with x(t) = sin α sin γ cos t + cos α cos γ.
The need to evaluate cos t for multiple values of t is eliminated by choosing a numerical integration algorithm which always evaluates the integrand at a large number of equidistant points along t. As the points are equidistant, cosine and sine need only be evaluated at the leftmost point and the interpoint distance, while the cosine values for all other points can then be generated from those using inexpensive angular addition formulas. Additionally, it is obviously desirable to choose an integration algorithm with a fast convergence, in order to keep the number of evaluations reasonable. The chosen algorithm is that of Romberg integration [26], but customised so as to always start by evaluating the integrated function at 17 equidistant points in [0, t ] in one go, and immediately constructing an accurate estimate from those. Thus, the first 4 layers of the traditional Romberg algorithm are combined in order to avoid evaluating the integrand at too few points at a time. In the rare cases where convergence is slow, additional function evaluations will be requested as the Romberg algorithm dives in to its 5th layer and beyond: first 16 additional equidistant points, then 32, etc. It is important to note that the very fast convergence rate of Romberg integration is realised only if the integrand is sufficiently smooth, which is why the lookup table for f (x) has to be implemented with a cubic spline, rather than a simpler scheme involving linear interpolation.
While the numerical integration of Eq. 24 is thus performed with utmost attention to efficiency, it is still preferable to carry out the evaluation with a cheaper closed-form expression. Although such an expression is not generally available, it is possible to find one which is a good approximation under certain very common conditions. Consider again Eq. 23: when δ(t) and t are both reasonably small, their cosines can be approximated with second order Taylor expansions: The magnitude of δ 0 ≡ δ(0) = |α − γ| is strictly smaller than or equal to δ(t), so we can expand the remaining cosine as well, and get: Using this, Eq. 24 can be approximated by: Performing a variable change u(t) ≡ √ t sin α sin γ/σ gives: With the error function, erf(x) ≡ 2π −1/2 x 0 exp(−s 2 )ds. Using δ(t ) = τ and Eq. 32, it follows that u(t ) = (τ 2 − δ 2 0 )/σ 2 , and thus: This result can be readily interpreted if one first considers the non-spherical case of integrating a two-dimensional non-truncated Gaussian field in a planar geometry along an infinite straight line. The result will depend only on the distance of closest approach, δ 0 , between the line and the centre of the Gaussian, and is trivially found to be exactly ( This factor, which we shall refer to as the simple Gaussian approximation, is present in Eq. 35, along with three correction factors for the spherical geometry and tail truncation: sin α/ sin γ is the lowest order correction for the curvature of the path of integration, the error function corrects for the truncation, and the factor of 2πσ 2 N represents a correction to the normalisation influenced by both of these effects. Of the three correction factors, the most interesting is arguably sin α/ sin γ, as it introduces asymmetries in rocking curves and shifts their peaks slightly away from the expected positions at γ = α. It is interesting to note that, using a different approach, J. Wuttke [27] also derived an approximation formula analogous to Eq. 35. This formula [27, Eq. 37] has a different functional form and does not properly account for normalisation (which is evident from [27, Fig. 5]), but it predicts an asymmetry factor which to lowest order in η = γ − α expands to the same result as the sin α/ sin γ factor derived here, namely 1 − 1 2 η cot α. Although somewhat complicated, Eq. 35 can be evaluated rather efficiently since everything except the factor of sin α/ sin γ depends only on δ 2 0 , which can be calculated from cos δ 0 . Thus, with S representing a suitable pre-calculated 1-dimensional cubic spline-based lookup table: As previously mentioned, cos α and cos γ are cheaply available for a particular incoming neutron via a multiplication and a dot product respectively, and their corresponding sinus values are calculated at the expense of a square-root evaluation. Likewise, cos δ 0 = cos(α − γ) = cos α cos γ + sin α sin γ, so the entire cross section evaluation through Eq. 36 consists of a few basic arithmetic manipulations, a cubic spline evaluation, and three square-root evaluations.
As Eq. 36 represents an approximation, the modelling code must choose when it should be employed, and when it is necessary to fall back to the full numerical integration of Eq. 24. The two approximations invoked above were both related to the expansion of the cosine function: Demanding that these are accurate to a certain precision, amounts to putting a limit, denoted ψ, on the magnitudes of the cosine arguments. If for instance a precision of = 10 −3 is required in all Taylor expansions of the cosines, the arguments must all have magnitudes less than ψ ≈ 22.6°. 5 The first approximation in Eq. 37 is valid at the desired precision for all 0 ≤ t < t only when τ < ψ, the validity of which is a global and fixed property of a given crystal.
Fortunately, τ < ψ is satisfied for most typical NCrystal configurations. On the other hand, it is always possible to find back-scattering scenarios for which the second approximation in Eq. 37 breaks down: when values of α and γ are low enough that the Bragg circle will be entirely located within an angle of τ fromn, t becomes equal to π. To be more precise, it follows from Eq. 28 that the condition t < ψ is equivalent to: cos ψ sin γ sin α + cos α cos γ < cos τ Which can indeed never be satisfied for the extreme back-scattering case of 5 Typical errors on final cross sections will actually in general be smaller than those indicated by . This is because the shape of the Gaussian mosaic density ensures that the largest contributions to the cross section will come from areas where the magnitudes of t and δ(t) are relatively smaller, and the Taylor expansions more accurate.
vanishing α and γ. So in conclusion, when τ < ψ and Eq. 38 is satisfied, the code will evaluate the cross section using Eq. 36, and otherwise it will fall back to the Romberg-based numerical integration outlined previously. In all cases providing results that are accurate at the desired level, defaulting to = 10 −3 .
As mentioned in Section 4.1, users of NCrystal can modify this default value through the configuration parameter named mosprec.

Sampling scattering events
In order to sample an outgoing neutron direction (k f ) in case of a scattering event, a particular hkl plane is first sampled randomly among those contributing, with the sampling based upon their contributions to the total cross section.
These contributions are usually already available, as they were determined and cached during a previous cross section calculation.
Having thus selected the particular hkl plane involved in the scattering, a crystallite orientation must then be sampled in order to provide an actual (as opposed to nominal) normal of the plane on which to scatter. This is first done in the coordinate system of Figure 4 by sampling a value of t in [−t , t ] according to the value of the integrand in Eq. 24 and using the resulting normal given by Eq. 21 to perform a specular reflection ofk i intok f . Finally, a suitable rotation matrix is applied to the resultingk f vector in order to rotate it into the coordinate system used by the calling code. The actual sampling of a t value proceeds via acceptance-rejection sampling [28], relying on the fact that the density function, the integrand of Eq. 24, attains its largest value at t = 0. This maximal value can thus be used as a constant overlay function for the rejection sampling. As splined lookup tables of the integrand were already prepared for the purpose of numerical integration of Eq. 24, these are reused during the sampling, which in practice is found to have acceptance ratios from 15%-35%, depending on the scenario. The resulting sampling speeds are acceptable, and usually insignificant, as will be discussed further in Section 7.

Effective model at very short d-spacing
As described in [1], great care is taken to construct lists of reflection planes is typical for most crystal structures, is that while the isotropic approximation ends up being applied to the vast majority of planes, the affected planes are all so weakly scattering that even their combined contribution only amounts to a tiny fraction of the cross section. As will be shown in Section 7, this in practice means that single crystal modelling is made 1-2 order of magnitudes faster for short wavelength neutrons, at very little cost of realism. Even in the most pessimistic cases, like the corresponding plot for a Corundum powder shown in Figure 5b, the trade-off between speed and accuracy provided by the approximation is likely to be appropriate for all practical use cases.

Validation
The work done to validate the single crystal model, as implemented in Secondly, it is verified in Section 4.7.3 that it is able to reproduce results from accepted existing models, in the domains of their validity. In addition to this, benchmark numbers for computational efficiency are presented and discussed in  attention was given to using not only "easy" configurations, but also those with back-scattering (α comparable to σ) and forward-scattering (π − α comparable to σ).
To verify first the evaluation of Eq. 24 a series of rocking curves were constructed for various scenarios, showing the effect of varying the neutron incidence angle γ, on the estimated cross sections through g hkl . Note that in order to remove trivial differences by design due to the -dependency of the relative mosaic truncation angle (cf. Figure 3), a fixed truncation angle of τ = 5σ was used for all curves. Three configurations will be discussed in the following, but figures with the results of other configurations are available in the appendix in Reference Gaussian approximation Improved approximation NCrystal default NCrystal = 10 −7     as ( √ 2πσ) −1 exp − 1 2 δ 2 0 /σ 2 is not able to provide very accurate results, while the improved approximation in Eq. 35 is able to provide results accurate to 3-4 digits. It is clear from the overlapping curves in Figure 6 that NCrystal with default settings is utilising this approximation. There is, however, a minor breakdown in accuracy around the edges at γ ≈ 20 • and γ ≈ 30 • , which is due to artefacts caused by the limitations of the cubic spline used in Eq. 36 to implement the evaluation of Eq. 35. As the decrease in accuracy happens only at the edges, where g hkl ≈ 0, this is an acceptable price to pay for the increased computational efficiency. When NCrystal is configured for increased accuracy, = 10 −7 , the code automatically falls back to numerical integration, with 9 significant digits correctness at all angles.
Next, Figure 7 shows the same scenario but with a much smaller mosaicity of σ = 1 . Here, both approximations have more ideal conditions, and the improved approximation of Eq. 35 is able to provide 7 significant digits over the entire range of incidence angles. Consequently, even when NCrystal is configured for increased accuracy, = 10 −7 , the approximation is utilised, completely foregoing any fall back to numerical integration. Again some artefacts appear Reference Gaussian approximation Improved approximation NCrystal default NCrystal = 10 −7 Finally, Figure 8 shows a back-scattering scenario. Although the improved approximation still outperforms the simple Gaussian approximation in this case, neither approximation is able to adequately reproduce the reference results.
Consequently, NCrystal falls back to numerical integration with highly accurate results. Again the usage of cubic splines, this time for the integrand of Eq. 24, leads to acceptable inaccuracies at the edge where g hkl ≈ 0.
In order to evaluate the scattering event sampling code described in Sec-   Finally, Figure 11 investigates a scenario with forward-scattering, using an even higher mosaicity of σ = 3 • for visualisation purposes. Scattering onn again produces a bell-curve around 0, but this time there is a non-zero contribution for scattering on −n as well, giving rise to events scattered to ϕ ≈ ±π. Also in this case, the NCrystal simulations perfectly reproduce the calculated reference.

General consistency checks
Irrespective of mosaic distribution, any model of Bragg diffraction which is implemented in a self-consistent manner and supports neutrons at any incidence angle, must be able to pass a few generic consistency checks.
First of all, the fundamental symmetry between (h, k, l) and (−h, −k, −l) planes means that a neutron which just scattered in a Bragg diffraction process on the (h, k, l) plane will have non-zero cross section for a subsequent scattering on the (−h, −k, −l) plane, bringing it back to its original direction. Thus, in the absence of other reflection planes or physics processes, a neutron in an infinite material will keep reflecting back and forth between the two planes in a "zig-zag" or "ping-pong" walk between the two sides of the planes. 7 It was verified for a 7 As was noted in Section 4.3, it is in scenarios involving extreme forward scattering possible for both (h, k, l) and (−h, −k, −l) to simultaneously contribute to the scattering cross section.
In such cases, the "zig-zag" walk will strictly speaking also include the occasional "zig-zig" or range of scenarios that this "zig-zag" walk is reproduced in NCrystal simulations by letting neutrons scatter 10 10 times in a wide range of scenarios, covering both large and small mosaicities as well as both forward-and back-scattering.  This is acceptable as it is still an order of magnitude more precise than the target O(10 −3 ) precision with that setting. Increasing the requested level of precision to = 10 −5 makes this discrepancy disappear, thus verifying that the source "zag-zag". Based on the Darwin-Hamilton equations [30], Sears [31] derived analytical closed-form solutions to the reflectivity under the restriction that wavevectors never scatter out of the plane spanned by the plane normal and initial neutron direction -essentially amounting to a requirement that the azimuthal scatter angle always be strictly 0 in the interactions. This is certainly an approximation at best, and one which completely breaks down for back-scattering and scattering at grazing incidence (cf. Figures 10 and 11). J. Wuttke [27] generalised the Darwin-Hamilton equations in order to lift this azimuthal restriction of the scattering direction, and under certain approximations developed analytical results for reflections on a slab. He also provided a Monte Carlo code able to perform numerical integrations needed to support also large mosaicities and provide a reference for his analytical expressions. That code does, however, also explicitly exclude back-scattering and scattering at grazing incidence.
In order to compare the outcome of NCrystal-based simulations with these existing solutions, a Monte Carlo stepping simulation was carefully set up, emulating the restrictions of the references, in order to predict idealised reflection probabilities: inelastic and incoherent elastic scattering events were treated in the same manner as absorption events, all but a single reflection plane was removed from the simulation, the normal of which was made to coincide with the surface normal of the slab. Finally, scenarios with back-scattering or grazingincidence scattering were avoided, and in particular values of neutron incidence and Bragg-angles which were also used in [27] were favoured -on the assumption that they might have been particularly well validated. In order to extract reference results, Sears' formulas were directly evaluated, while Wuttke's Monte Carlo programme was downloaded and executed. Where relevant, parameters were set to emulate scattering on a Germanium-511 plane.
First, the rocking curves in Figure 14 show the resulting reflection probabilities as a function of neutron incidence angle for a FWHM mosaicity of 0.   As a complement to the rocking curves, Figure 17 shows the reflection probabilities as a function of slab thickness, for neutron incidence and θ B both 80 • and for 3 different mosaicities. Although the discrepancies between the different models are less pronounced than in the rocking curves, the qualitative conclu-

Layered single crystals for pyrolytic graphite
The isotropic Gaussian mosaicity distribution discussed in Section 4 provides an appropriate description for modelling of a wide range of single crystal materials encountered at neutron scattering facilities and elsewhere. However, exceptions do exist, with some materials exhibiting radically different crystallite distributions. One such material which is of particular importance for instrumentation at neutron scattering facilities is pyrolytic graphite, with important and frequent deployment as both monochromators and analysers targeted at selecting longer wavelengths [32,33], as well as tunable beam filters intended to remove short wavelength contamination in monochromated spectra [34,35].
The structure of graphite can be described as stacks of two-dimensional graphene sheets, in which carbon atoms are bound strongly into hexagonal grids, with a weaker binding between the sheets. This layout at the microscopic level ultimately gives rise to anisotropic features in the macroscopic mosaicity distributions. With suitable manufacturing methods [36,37], the axis orthogonal to the graphene sheets, the lattice c-axis, will be distributed along a preferred direction, suitable for description with a Gaussian mosaicity distribution like the one discussed in Section 4.1. The rotation around this axis, defined for instance by the direction of the lattice a-axis, will however be uncorrelated between different crystallites. Thus, the associated rotation angle will be uniformly distributed in [−π, π], giving rise to features in scattering interactions akin to those observed with crystal powders.
In order for NCrystal to model pyrolytic graphite, a specialised model of layered single crystals is implemented as described in Sections 5.1-5.3. An alternative reference model for validation work is described in Section 5.4, and actual validation work is then carried out in Section 5.5.
It is worth noting that while the presented models were developed in order to support studies involving pyrolytic graphite [38], it could in principle also be used to model other materials like hexagonal boron nitride [39] which exhibits similar layered structure. Under certain conditions relating to the neutron velocity and rotational speed, it could even be used to describe an isotropic Gaussian single crystal which is placed on a spinning sample holder.

Geometrical integral and definitions
NCrystal effectively models a layered single crystal as consisting of a large number of small crystals with isotropic Gaussian mosaicity distributions, but Specifically, if g SC hkl (α, γ) represents the geometrical weight for scattering on a given crystal plane in a crystal with an isotropic Gaussian mosaicity distribution given by Eq. 24, then the equivalent factor in a layered single crystal is defined to be: Where ϕ represents the rotation aroundL, the offset of which is in principle arbitrary and can be defined in whatever manner is most convenient. For instance, one could define it so that the lattice a-axis in pyrolytic graphite would coincide with a certain direction in the laboratory frame when ϕ = 0. However, the present discussion will instead adopt a convention relative to the direction of the neutron, k i , and specific to each hkl-plane, which minimises the distance between the plane normal and − k i at ϕ = 0. More specifically, the particular coordinate system shown in Figure 18 is adopted:L is placed on the positive z-axis, andk i is placed in the xz-plane with coordinates: Where θ k is the angle between −k i andL. If θ n is the angle betweenn andL, the nominal coordinates of the normals are then given by: n(ϕ) = (sin θ n cos ϕ, sin θ n sin ϕ, cos θ n ) With these definitions, γ is straightforward to determine from the dot product of −k i andn(ϕ): cos γ = sin θ k sin θ n cos ϕ + cos θ n cos θ k Note that, due to the symmetries between the (h, k, l) and (−h, −k, −l) reflection planes, the physics should be unchanged under the transformation L → −L, and as long as the code is written so as to deal with the two sets of equivalent planes at θ n and π − θ n jointly, it can therefore always be assumed for simplicity that both θ k and θ n lie in [0, π/2]. Additionally, all hkl-planes with similar θ n and d-spacing can be treated concurrently, as a single group of planes whose combined scattering strength is simply the sum of the contributing planes q hkl values, thus effectively reducing both memory usage and time spent processing planes during simulations. Additionally, the joint processing of the groups of planes at θ n and π − θ n benefits from the fact that the group of planes at π − θ n can only ever contribute if the group θ n contributes. For simplicity, the discussions will from this point onward mostly ignore direct mention of the planes at π − θ n , as the required modifications to the equations are mostly trivial. Similarly, the discussions will ignore the degenerate cases θ n ≈ 0 and θ k ≈ 0, as they are trivially solvable: the θ n ≈ 0 case (relevant for 00l planes in pyrolytic graphite) can be treated with the model developed in Section 4, and when θ k ≈ 0 the integrand in Eq. 39 becomes independent of ϕ.
The particular coordinate system defined in Figure 18  angle between −k i and plane normals satisfying the condition for Bragg diffraction is given as α = π/2 − θ B , defining the Bragg circle (red). Finally, θ k is the angle between −k i andL, and γ is the angle betweenn(ϕ) and − k i . restrict ϕ to the interval [0, π] and use instead: Additionally, either the integrand in Eq. 43 will be vanishing everywhere, or it will be non-vanishing exactly on a single sub-interval, [ϕ min , ϕ max ], which will be identified in Section 5.2. A final advantage of the chosen coordinate system is that it allows an efficient pre-check of whether or not a given set of planes will contribute to the scattering of a particular neutron, by determining whether or not the z-range of the Bragg circle overlaps with the z-range of the mosaicity bands. The latter can be trivially pre-calculated at initialisation time and the former is given by: Which can be calculated cheaply at simulation time by usage of cosine addition formulas, using cos α = λ × (1/2d) and sin α = √ 1 − cos 2 α. In addition to a few basic arithmetic manipulations, the cost of the pre-check is therefore just a simple square-root evaluation. To further speed up this critical section of the code, the check is carried out first using a Taylor expansion approximation of the square-root evaluation.

Evaluating the geometrical integral
In principle Eq. 43 can be evaluated directly via numerical integration, using the algorithms described in Section 4 to evaluate g SC hkl at each ϕ point. However, for non-vanishing α and realistic mosaic spreads, this integrand will be vanishing at most ϕ values, with non-zero contributions arising at most in a single narrow region. For numerical integration to proceed reliably and efficiently, it is therefore necessary to predetermine the ϕ-regions in which the integrand is non-zero analytically, and apply the numerical integration algorithm only to these regions. In this context it is also worth noting the importance of being able to evaluate the integrand of Eq. 43 with the efficient method of Eq. 36, thus most of the time avoiding the computationally demanding scenario of having a numerical integration of an integrand which is itself evaluated with numerical integration. Now, a neutron has a non-zero cross section to scatter on a small crystal with rotation ϕ if and only if: Since γ ≥ 0 by definition, this is equivalent to: Now, all three parts of this double inequality resides in [0, π], since by definition γ ≤ π, α ≤ π/2, and τ ≤ π/2 (for the case of τ this is a deliberate restriction, as discussed in Section 4.1). On this domain, the cosine function is one-to-one and with non-positive derivative, so Eq. 46 is equivalent to: Using Eq. 42 this becomes: cos(max(0, α − τ )) ≥ sin θ k sin θ n cos ϕ + cos θ k cos θ n ≥ cos(α + τ ) (48) As the degenerate cases θ n ≈ 0 and θ k ≈ 0 are dealt with separately, it is safe to divide with sin θ k sin θ n : cos(α + τ ) − cos θ k cos θ n sin θ k sin θ n ≤ cos ϕ ≤ cos(max(0, α − τ )) − cos θ k cos θ n sin θ k sin θ n Since ϕ ∈ [0, π], the range of phi values representing crystallite orientations with non-zero contributions to the scattering cross section is thus ϕ min < ϕ < ϕ max , with: ϕ min ≡ arccos min 1, cos(max(0, α − τ )) − cos θ k cos θ n sin θ k sin θ n ϕ max ≡ arccos max −1, cos(α + τ ) − cos θ k cos θ n sin θ k sin θ n (50) Obviously, g LC hkl (α,k i ) = 0 unless ϕ max > ϕ min , in which case Eq. 43 can be replaced by: Now, the code developed in Section 4 for evaluation of g SC hkl (α, γ(θ k , θ n , ϕ)) works directly from cos γ rather than γ for reasons of efficiency, which with Eq. 42 translates into a need to provide a value of cos ϕ at each evaluated ϕ point.
As in Section 4.4, such values are cheaply generated via trigonometric addition formulas, and the customised implementation of Romberg integration discussed in Section 4.4 is reused for efficient and accurate integration of Eq. 51 as well. Figure 19 shows the resulting Bragg diffraction cross sections in pyrolytic graphite as a function of neutron wavelength and θ k , with distinctive structures corresponding to particular groups of scattering planes. In order to understand how these structures relate to actual scattering planes, it is instructive to consider a hypothetical layered single crystal with just a single reflection plane having θ n = 60 • . If both "sides" of this reflection plane are included, i.e. both (h, k, l) and (−h, −k, −l) are considered (as is always the case in NCrystal), the resulting cross section structure is shown in Figure 20. The decomposition into contributions from scattering on (h, k, l) (θ n ≤ π/2) and (−h, −k, −l) (θ n ≥ π/2) is shown in Figure 21. Considering first Figure 21a, its structure can be understood (at least in the limit of small mosaicities) by carefully considering Figure 18 while recalling that cos α = λ/2d hkl . The present discussion will limit itself to consideration of a few important features. Firstly, when θ k ≈ 0, the neutron is aligned alongL and the only contribution comes when α ≈ θ n , which means λ/2d hkl ≈ cos θ n , translating to λ ≈ 0.5 Å. Next, in case of back-scattering, λ → 2d hkl and the Bragg circle contracts to a point, the only contribution comes when θ k ≈ θ n . Keeping θ k ≈ θ n but lowering λ, the Bragg circle will grow and the integral along it through the mosaic density, g hkl , will quickly tend to a constant, which is reminiscent of -but not identical tothe formula for a crystal powder in Eq. 11. A constant g hkl inserted into Eq. 8 implies σ hkl ∝ λ 2 / 1 − (λ/2d hkl ) 2 , which can be compared with the equiva-  lent prediction in a powder, σ hkl ∝ λ 2 . When θ k < π/2 − θ n , the situation is a bit different: here, −k i lies so much "to the north" of the mosaicity band in Figure 18, that when the Bragg circle is at its largest, it will lie "to the south" of the mosaicity band everywhere, which explains the vanishing cross sections in the bottom left of Figure 21a. The ridge along this triangle shows enhanced cross sections, which can be understood from an intermittent enhanced overlap between the Bragg circle and the mosaicity band for rotations around ϕ = π, just before the Bragg circle extends too far south. Finally, it might be worth noticing that at θ k = 90 • , cross sections will vanish unless the Bragg circle is large enough to reach the mosaicity band, which happens when α = π/2 − θ n .
I.e, the right-most ridge in Figure 21a intersects the axis at θ k = 90 • at the point where λ = 2d hkl sin θ n , which in this case corresponds to λ ≈ 0.866 Å.
Moving on to Figure 21b and the contribution from scattering on the mosaicity band at π − θ n , such scattering is impossible unless θ k is large enough that the Bragg circle can reach the mosaicity band on the "southern hemisphere" (note that this band is not explicitly shown in Figure 18). In other words, contributions only exist in the region defined by θ k > θ n − α. The edges of this region intersect the axes at λ = 0 and θ k = 90 • in the same locations as the ridges in Figure 21b, thus completing the distorted triangular structure visible in Figure 20. The ridge-like structure along the edge in Figure 21b can again be understood as coming from an enhanced overlap of the Bragg circle and the mosaicity band near the edge. Away from the edge in Figure 21b and inside the region in the upper left corner, the cross section will again fall off as Qualitatively, the distorted triangular structure seen in Figure 21   The structures seen in Figure 19 can now be understood as arising from the various values of (θ n , d hkl ) encountered in pyrolytic graphite. In Section 5.5.3, the connection to specific hkl values will be explored further.

Sampling scattering events
In case of a scattering event in a layered single crystal, the sampling of an outgoing neutron direction (k f ) starts of in the same vein as in the case of nonlayered single crystals discussed in Section 4.5: a set of hkl planes with common scenarios, and additionally the code is written so as to throw an exception in case an invalid overlay height is ever encountered (obviously no reports of this has surfaced so far). The acceptance rate seems to be around 30-40%, which implies that the sampling requires roughly 11-12 evaluations of the integrand of Eq. 39: 9 to construct the overlay curve, and around 2-3 to carry out the sampling.
Having sampled a particular value of ϕ and thereby determined the exact orientation of the hypothetical small Gaussian mosaic crystal on which to scatter, the code discussed in Section 4.5 is used to generate the actual scatteringstill in the coordinate system of Figure 18. Finally, the resulting neutron direc- (a)  (b) Figure 24: Examples of overlay and acceptance functions as in Figure 23, but for examples involving situations with both small (Figure (a)) and large (Figure (b)) mosaic spreads.
tion,k f must be rotated into the laboratory frame, which involves the careful determination of an appropriate rotation matrix based on the values ofL and k i in the laboratory frame.

Alternative reference implementation
For reasons of efficiency, the implementation described in Sections 5.1-5.3 groups hkl planes according to d-spacing and θ n , and then treats each such group separately and in its own coordinate system. This approach takes advantage of rotational symmetries in these coordinate systems, exploits symmetries between the (h, k, l) and (−h, −k, −l) reflection planes, and ensures that numerical integration algorithms will only be deployed on smoothly varying and non-vanishing functions.
As an alternative to this efficient but admittedly complicated implementation, a simpler reference model is useful for validation and comparison purposes -even if it is prohibitively inefficient for practical usage. Fortunately, such a model is readily implemented, based again on the notion of modelling a layered single crystal as consisting of a large number of small non-layered single crystals, rotated uniformly aroundL. In the reference implementation presented here, the small crystals are instead treated directly with the code for non-layered single crystals described in Section 4. Thus, with {ϕ j } representing a set of N j rotations, ideally very large in number and distributed uniformly over [−π, π], the cross section for coherent elastic scattering of a neutron with wavevector k i on a layered single crystal is modelled as: Where R ϕj is a rotation operator implementing a rotation of ϕ j aroundL, with ϕ = 0 meaning no rotation, and σ SC ( k i ) is the single crystal cross section. In order to sample the outcome of a scattering event, a particular ϕ j is first sampled according to the contribution of the corresponding term in Eq. 52, and the nonlayered single crystal model is then used to scatter R ϕj k i into a final state k f , which is subsequently transformed back to the original frame to get the actual outcome of the interaction: k f = R −1 ϕj k f . For simplicity, and to avoid issues of numerical stability, the rotation operator and its inverse are implemented using Rodriguez' rotation formula [40].
The only remaining issue is how to construct the set of rotational values, {ϕ j }. One approach is the random sampling of n uniformly distributed values in [−π, π] -with a different set generated for each incoming neutron to reduce statistical bias. Alternatively, the n values can simply be distributed equidistantly over the interval [−π, π]. In the latter case, it is clear that for the model to yield reliable results, n must be large enough that π/n will be much smaller than the mosaicity of the crystal. In the former case of randomised ϕ values, n will likely need to be much higher since the sampled ϕ j points tend to cluster in some regions, leaving other regions with reduced density. On the other hand, one might worry that the structure inherent to a non-randomised set of ϕ j values might be reflected as an unwanted artefact in the final distributions of k f values, but in principle it should not be a significant problem as long as n is high enough that angular effects at the order of π/n are washed out by the non-vanishing mosaicity effects in the non-layered single crystal models.
Nevertheless, NCrystal supports both models through the lcmode configuration parameter. If left at the default setting, lcmode=0, the efficient model described in Sections 5.1-5.3 is used. If set to a positive value, lcmode=n, the reference model described in the present section will be used with n values of ϕ, evenly spaced in [−π, π]. If set to a negative value, lcmode=-n, the reference model will again be used, but now with the n values of ϕ sampled randomly and independently for each neutron, as described above. It is interesting to note that the Single_crystal component [19] of McStas since 2015 supports an experimental mode for pyrolytic graphite modelling, which essentially corresponds to NCrystal's lcmode=-1. As should be clear from the present discussions, and which will be verified in section Section 5.5.3, |n| = 1 is much too low to provide reliable results.
In general, the validation work in Section 5.5.1 will employ an evenly spaced distribution of ϕ j values, with a high value of n. Requiring 2π/n to be equal to 5% of the FWHM mosaicity of a crystal, one finds n = 7.2 × 10 3 for a mosaicity of 1 • , n = 4.32 × 10 5 for a mosaicity of 1 , and n = 2.592 × 10 6 for a mosaicity of 10 . In practice this limits the usage of the reference models, even for validation work, to crystals with large mosaicities.

Validation
As was the case for the non-layered single crystal model in Section 4.7, the work done to validate the layered single crystal model first consists of a verification in Sections 5.5.1 and 5.5.2 that the implementation described in

Comparison with reference implementation
Due to the double-integration involved, it is unfortunately impractical, in terms of computational time requirements, to proceed in a similar vein as in Section 4.7.1 and provide a high-precision reference implementation for layered single crystals with mpmath. Instead, the alternative model described in Section 5.4, available in NCrystal via the lcmode configuration parameter, will be used to provide reference results. These are then used to verify the implementation of the primary model described in Sections 5.1-5.3.
First, Figure 25 shows the scattering cross section as a function of wavelength for a large FWHM mosaicity of 3 • , and a neutron incidence angle of θ k = 40 • .
As is clear from the relative difference curves, the results are in good agreement with the reference, at a level which is everywhere equal to or exceeding the requested precision. The only exception is a degradation for the = 10 −3 curve near the edges of the truncated mosaicity -in particular around 3.7Å. As was the case in Section 4.7.1, this degradation is due to artefacts related to the usage of cubic splines. Also in this case it is deemed to be an acceptable price to pay  for the increased computational efficiency, given that it happens only in regions where the cross section is almost negligible anyway.
Next, Figure 26 shows the same curves, but this time for a much smaller FWHM mosaic spread of 0.01 • . As expected, the structure in the curves is more sharply defined at the reduced mosaicity, and the level of precision in the curves are still at an acceptable level. The degradation in precision near the edges seems less significant than in Figure 25  In order to evaluate the scattering event sampling code described in Section 5.3, both the default NCrystal model and the reference model (with 10 4 evenly distributed angular rotations) are used to sample scattering angles in pyrolytic graphite with a FWHM mosaic spread of 1 • , for neutron wavelengths uniformly distributed between 2.5Å and 3.5Å, and a fixed neutron incidence of θ k = 65 • . As can be inferred from Figure 19, this setup exercises all parts of the sampling procedure, as it involves scattering on multiple reflection planes, with both θ n = 0, 0 < θ n < π/2, and π/2 < θ n < π (corresponding qualitatively to Figures 22a, 21a, and 21b respectively), and involves reflections both with and without back-scattering. The resulting distribution of scattered directions by the default NCrystal model is shown in Figure 28, which was confirmed visually to be indistinguishable from the same distribution created by the reference model. For a more quantitative comparison, Figures 29 and 30

Comparison with existing results
Unlike the situation concerning non-layered single crystals (cf. Section 4.7.3), no existing code or analytical models available to the authors provides reliable and precise results for neutron scattering on layered crystals like pyrolytic graphite. However, E. Frikkee [35] provides analytical predictions of the cross section structure, parameterising the peak positions visible in Figure 19. to scattering on specific hkl reflection planes, E. Frikkee was not able to perform a more quantitative analysis of the data. However, with NCrystal it is now straightforward to reproduce the setup in a dedicated simulation. The exercise benefits from the fact that E. Frikkee corrected the observed transmission spectrum for contaminations due to higher-order λ/n reflections in the upstream monochromator, and had collimators placed before the slab and between the slab and detectors. The result is shown in Figure 33, with a remarkable agreement between the data points and the predictions of NCrystal. It is worth noting that, as can be inferred from Figure Where the summation index j runs over all atomic positions in the unit cell, the constant σ inc j is the incoherent scattering cross section of the atom in question, Q ≡ k f − k i is the momentum transfer, and W j ( Q) the Debye-Waller function.
True to the nature of incoherent scattering, Eq. 53 shows no interference terms between atoms at different positions. For isotropic and harmonic displacements, the Debye-Waller function is given by 1 2 Q 2 δ 2 j , where δ 2 j is the mean-squared displacement of the jth atom from its nominal position in the crystal, as seen over a large ensemble of unit cells. All in all, the contribution of atoms occupying the jth position in the unit cell, is simply given as: To get the cross section at a given neutron wavelength, λ = 2π/k, one must integrate Eq. 54 over all outgoing directions of k f . Designating the scattering angle as θ, and defining µ = cos θ, one finds: Now, dΩ = sin θdθdϕ = dµdϕ and integration over dϕ merely yields a factor of 2π, so: Where the parameter t, was introduced: For reasons of numerical efficiency and stability, Eq. 56 is in NCrystal evaluated with a third-order Taylor expansion when t < 0.01, and with the limiting expression σ inc /t when t > 24. In order to sample a value of µ for a particular scattering, a probability density proportional to the integrand in Eq. 56 must be used. This implies that a µ value must be sampled in [−1, 1] according to the distribution: With the normalisation factor N t = t/(4 sinh(t/2)). For reasons of numerical stability and computational efficiency, when t < 0.02 the implementation in NCrystal samples µ via straightforward acceptance-rejection sampling using a constant overlay function and a Taylor expansion for evaluating Eq. 58.
For larger values of t the transformation method is used instead, yielding µ = −1 + 2t −1 log(1 + (e t − 1)R) where R is a pseudo-random number distributed uniformly over the unit interval. As is the case for the implementation of Bragg diffraction in a crystal powder discussed in Section 3, the incoherent elastic model is also fast enough that there is potentially a non-neglible overhead from the construction of complete directional vectors in NCrystal's vector interface for scattering event sampling. Consequently, this interface is implemented using the same efficient methods for vector construction as those discussed in Section 3.
Coming back to the significance of the parameter t: when λ 4πδ, t approaches 0 and incoherent elastic scattering becomes isotropic with a wavelengthindependent cross section -which is indeed a widely used approximation used to model incoherent elastic scattering. The approximation is, however, worse for materials with large atomic displacements which could for instance be a result of high material temperature, and it will always eventually break down when the wavelength of the incoming neutron is small enough. Figure     shows the importance of incoherent elastic scattering for Vanadium and Nickel, and illustrates how the model as implemented in NCrystal is essential for reproducing the measured total cross sections. Due to the general low interest in incoherent elastic scattering in itself, no suitable measurements were found with which a corresponding comparison of scattering angles could be performed.
In particular, the search for scattering angle data sensitive to deviations from isotropicity with vanadium samples is complicated by the fact that many neutron scattering instruments are calibrated with vanadium samples.

Computational efficiency
The computational cost of using NCrystal to provide information about neutron scatterings is of course highly dependent on the use case: specific materials, distribution of neutrons, and configuration options can all influence this greatly. In order to nonetheless gauge this, the present discussion will use a benchmark in which an isotropic source of neutrons at a given wavelength interacts in a material which is either a powder, a single crystal, or a layered single crystal. The direction-dependent interfaces are used in all cases, as this is what is typically used when NCrystal is embedded in a generic simulation framework like Geant4 or McStas. To enable more meaningful comparisons, the crystal structure of the material will in all cases be that of pyrolytic graphite. This is important, since the comparisons would otherwise be influenced by materialspecific details such as the number of reflection planes and their d-spacings (see for instance [1, Fig. 4]). The timings were in all cases carried out using an otherwise unoccupied node at the ESS-DMSC computing cluster in Copenhagen on a 2.40 GHz Intel ® Xeon ® Processor (E5-2680 v4). This of course represents a worst case estimation of the impact of scattering Bragg diffraction in a powder. This is not surprising since, as mentioned in Section 6, they employ the same algorithm for finding the directional vector of the outgoing neutron on the basis of the sampled scattering angle. It is, however, not inconceivable that a faster algorithm could be adopted for longer wavelengths where incoherent elastic scattering is almost isotropic. Given that the incoherent elastic sampling rate is already high, this was not deemed a high priority so far -but such an improvement might be included in a future NCrystal release.
The curve for Bragg diffraction in a powder in Figure 37 shows what is essentially constant performance below the Bragg cutoff of 6.71 Å, corresponding to a rate of O(50 MHz). This is not surprising, given the implementation discussed in Section 3, whose main overhead is a binary search in a single contiguous array. As already mentioned, the sampling speed of O(10 MHz) is dominated by the construction of the directional vector of the outgoing neutron.
Turning next to the single crystal models, layered or not, the situation is more complicated. Starting at the Bragg cutoff at 6.71 Å, more and more reflection planes are able to satisfy the Bragg condition as the neutron wavelength decreases. Unlike the case of a powder, each (group of) reflection planes must be checked separately, leading to a growth all the way down to 0.8 Å, where the isotropic approximation described in Section 4.6 prevents additional planes from being considered. The dashed lines in Figure 37 show the effect of foregoing this approximation, leading to a continued growth of evaluation times down to 0.2 Å. The final plateau below 0.2 Å is due to the other threshold in NCrystal, which completely leaves out reflection planes with d-spacing below 0.1 Å. It is worth noting that this second threshold does not in general represent a trade-off between realism and computational efficiency, as it was shown in [1, Fig. 6] that the reflection planes removed by it do not actually contribute significantly to the cross sections.
Next, the differences between the various single crystal curves in Figure 37 can mainly be attributed to two effects. The first of these is that layered single crystals are about an order of magnitude slower to process than non-layered single crystals, which is hardly surprising given that the former must essentially If deemed a reasonable trade-off for a specific use case, it is of course possible to increase the threshold at which the isotropic approximation is applied, through the sccutoff configuration parameter, and thus improve the computational speed at lower wavelengths. Concerning the speed of scattering event sampling, it is in general seen to be almost negligible compared to the cross section evaluation time. This is especially true at shorter wavelengths where more planes must be considered, and is hardly surprising since the code sampling scattering events can benefit from the work already done when calculating the cross sections, in order to quickly select just a single plane with which to proceed. As already mentioned, the treatment of layered single crystals is roughly an order of magnitude slower than that for single crystals. The exception in Figure 37 is for neutron wavelengths above 4.27 Å, where only the 002 plane contributes. As this plane is aligned with the lattice c-axis, it has θ n = 0, and is in practice treated with code similar to that used in non-layered single crystals.
Finally, it should be noted that the alternative implementations of layered single crystal code discussed in Section 5.4 would be at the very least 3-4 orders of magnitude slower than the model used in Figure 37, i.e. the one described in Sections 5.1-5.3, thus rendering it unusable for practical simulation work.

Summary and outlook
The models presented in Sections 3-6 cover the capabilities for elastic thermal neutron scattering in NCrystal releases v1.0.0-v2.2.1. Along with the scattering kernel-based inelastic scattering models introduced in NCrystal release v2.0.0, the elastic models enable realistic representation of the most important crystalline materials found at neutron scattering facilities. The implementations in particular focus on the primary intended usage of NCrystal: to serve as a physics backend in generic simulation frameworks like Geant4. Consequently, they focus on ease of configuration, computational efficiency, and in particular robustness -which means that the code provides accurate results for all supported configurations and neutron state parameters. Thus, the code supports running with a very large number of reflection planes at shorter wavelengths, both small and large mosaic spreads in single crystals, and both extreme backwards or forward scattering. Failure to support some of these cases correctly could lead to surprises and misleading results for casual users, whose simulations might just happen to involve one or more of them, which would be unacceptable.
Obviously, the presented models do not account for all effects which one might find in real crystalline materials, and additional features might in the future be incorporated in order to reproduce specific characteristics of some materials -in particular concerning detailed simulations of objects placed in the sample position at a neutron scattering instrument, which by design can be more sensitive to small effects. Examples of such effects might include the ability to account for effects of the finite crystallite grain sizes, fluctuations in lattice parameters, site occupancies, chemical disorder (e.g. in doped crystals), material strain, longer-scale fluctuations in scattering length densities (SANS models), and texture effects in polycrystals. Additionally, while support for multi-phase materials is in principle already included, it currently relies on user-code for handling the multiple phases. A future release of NCrystal should make it possible to define multiple phases directly in the NCrystal configuration strings -thus ensuring easy and consistent setup for all users. Finally, it is of course also desirable to continue to expand NCrystal's library of validated material files, and to improve the performance of existing models.
In addition to developments in capabilities for physics modelling, several potential improvements to NCrystal of a more technical nature are also desirable.
Firstly, it is planned in the near future to make the framework multi-thread safe and provide optional interfaces suited for parallel computing. Secondly, it is obviously very desirable for NCrystal to be integrated for use as a backend in more simulation packages than the existing, Geant4 and McStas. Obvious candidates for which some interest has already been expressed in the community are MCNP [42,43,44,45,46], OpenMC [47], PHITS [48], and VITESS [20,49].
Which of these many additional features will be available in the future, will depend on the needs of the community as well as the availability of manpower.
In this respect it might be interesting to note that NCrystal release v2.2.0 introduced support for user-contributed plugins, which will hopefully facilitate contributions to the development of new physics models from a wider community.