Monte Carlo simulation of neutron scattering by a textured polycrystal

Laliena, V.; Vicente-Álvarez, M.Á.; Campo, J.

doi:10.1107/S1600576720002290

research papers

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 53| Part 2| April 2020| Pages 512-529

https://doi.org/10.1107/S1600576720002290

Open

access

Monte Carlo simulation of neutron scattering by a textured polycrystal

Victor Laliena,^a ^* Miguel Ángel Vicente-Álvarez ^b and Javier Campo ^a

^aInstituto de Ciencia de Materiales de Aragón (CSIC – Universidad de Zaragoza) and Departamento de Física de Materia Condensada, Universidad de Zaragoza, C/Pedro Cerbuna 12, E-50009 Zaragoza, Spain, and ^bDepartamento Física de Neutrones, LAHN project, Centro Atómico Bariloche, CNEA/CONICET, SC de Bariloche, Argentina
^*Correspondence e-mail: laliena@unizar.es

Edited by G. J. McIntyre, Australian Nuclear Science and Technology Organisation, Lucas Heights, Australia (Received 9 September 2019; accepted 18 February 2020; online 30 March 2020)

A method of simulating the neutron scattering by a textured polycrystal is presented. It is based on an expansion of the scattering cross sections in terms of the spherical harmonics of the incident and scattering directions, which is derived from the generalized Fourier expansion of the polycrystal orientation distribution function. The method has been implemented in a Monte Carlo code as a component of the McStas software package, and it has been validated by computing some pole figures of a Zircaloy-4 plate and a Zr–2.5Nb pressure tube, and by simulating an ideal transmission experiment. The code can be used to estimate the background generated by components of neutron instruments such as pressure cells, whose walls are made of alloys with significant crystallographic texture. As a first application, the effect of texture on the signal-to-noise ratio was studied in a simple model of a diffraction experiment, in which a sample is placed inside a pressure cell made of a zirconium alloy. With this setting, the results of two simulations were compared: one in which the pressure-cell wall has a uniform distribution of grain orientations, and another in which the pressure cell has the texture of a Zr–2.5Nb pressure tube. The results showed that the effect of the texture of the pressure cell on the noise of a diffractogram is very important. Thus, the signal-to-noise ratio can be controlled by appropriate choice of the texture of the pressure-cell walls.

Keywords: neutron diffraction; crystallographic texture; Monte Carlo simulations.

1. Introduction

Neutron and X-ray scattering are extensively used in materials science for many purposes, in particular to analyse the structure of phases, quantifying their volume fractions and determining the state of stress and the crystallographic texture. The continuous demand for these techniques by the technological and scientific community gave rise to the construction of dedicated instruments at neutron and synchrotron facilities. Because of the low flux of neutrons compared with X-rays, in neutron laboratories the instruments are optimized for a particular set of scientific applications, which implies looking for the highest flux on the sample while keeping the resolution required by the technique to ensure a reasonable signal-to-noise ratio (SNR). The experimental setup determined by the optimization defines the characteristics of the beam impinging on the sample, which in turn influences the measurements, for instance the shape and position of the diffraction peaks (Mikula et al., 1997 ; Stoica et al., 2001 ).

Knowing how the instrument configuration affects the measurements is important not only during the design process of the instrument but also during operation, to interpret the bias of the experimental observations. The large number of variables that define the instrument configuration gives rise to an increasing use of Monte Carlo simulations. In these models, the neutron travels from the source to the detector and in its path interacts with the different components of the instrument and, eventually, with the sample. Examples of such software are the McStas package (Lefmann & Nielsen, 1999 ), the VITESS project (Zsigmond et al., 2002 ), IB (Zhao, 2011 ) and IDEAS (Lee & Wang, 2002 ), among others. Monte Carlo engines have also been added to some analysis programs like RESTRAX (Šaroun & Kulda, 1997 ) and have been used to estimate the corrections needed to extract physical quantities from experimental measurements (Vickery et al., 2013 ). A good example is the estimation of pseudo-stresses in neutron diffraction experiments (Šaroun et al., 2013 ). Attempts at realistic simulations by combining detailed instrument and sample modelling were presented by Farhi et al. (2009 ) and by Lin et al. (2016 ).

Monte Carlo simulations are particularly important for neutron instruments due to the large gauge volume necessary to have a significant signal, the reason being the low brightness of neutron sources compared with synchrotron or even laboratory X-ray instruments. This large volume brings about unwanted spatial resolution effects called pseudo-strains, which are caused by perturbation of the instrumental gauge volume due to the heterogeneous distribution of the scattering probability in the sample. The surface effect when the gauge volume is only partially immersed in the material is a well known special case. In general, any heterogeneity or beam extinction mechanism which causes significant variation of the scattering probability on a distance comparable to the gauge size may give rise to pseudo-strains, for example gradients in phase composition and texture, or a strong variation in beam attenuation with wavelength near a Bragg edge. The pseudo-strains are often of the same magnitude as the measured lattice strain and need to be properly treated. Monte Carlo models proved to be useful for this objective since they can account for beam attenuation, multiple scattering, divergence effects etc.

Another important application of Monte Carlo modelling is related to estimation of the SNR. In some cases, for example when using sample environment devices like furnaces or pressure cells, the neutron travels through the device wall before reaching the sample and/or the detector. In its path, it may suffer multiple scattering, either elastic or inelastic, increasing the instrument background. This is particularly important in high-pressure neutron instruments, where the pressure cell has a thick wall (Rodríguez-Velamazán et al., 2011 ; Rodríguez-Velamazán & Noguera, 2011 ). To reduce the background as much as possible, the selection of materials and their fabrication processes are critical. Alloys that minimize the background, such as TiAlV or CuBe, have been proposed (Kibble et al., 2019 ). However, to lower the background further, the crystallographic texture can be considered a design variable.

The scattering of neutrons by textured polycrystals, including a detailed description of texture, has not yet been fully incorporated into the available Monte Carlo programs. The nxs library to compute the neutron total scattering cross sections (Boin, 2012 ), which uses the March–Dollase model (Dollase, 1986 ) to include the effect of preferred grain orientations in the amplitude of Bragg edges, was implemented in McStas Release 2.5. Concerning the analysis of transmission (Bragg edge) spectra in textured materials, the total coherent cross sections have been implemented in terms of integration of the pole figures (Santisteban et al., 2012 ; Malamud et al., 2014 ), but these are not suitable for implementation in an efficient Monte Carlo code due to the demanding computational cost. Other tools for analysing transmission data through polycrystalline samples which implement approximations to the total cross section are Sinpol (Dessieux et al., 2018 , 2019 ) and RITS (Sato et al., 2011 ), and earlier work includes that of Vogel (1999 ).

In this work, we present expressions for the differential and total elastic coherent cross sections in terms of the generalized Fourier coefficients of the orientation distribution function (ODF), which are suitable for implementation in Monte Carlo programs. A closed expression for the total cross section derived here allows a time-efficient evaluation of this quantity, a necessary condition for its use in Monte Carlo simulations. As mentioned above, other expressions for this quantity were obtained earlier (Santisteban et al., 2006 , 2012). In our case, the truncation of the generalized Fourier series of the ODF renders the Monte Carlo simulations feasible, although they are computationally much more expensive than the standard simulations with single crystals or powder materials. The efficiency can be greatly improved by using variance reduction techniques. These developments have been implemented in a Monte Carlo code as a new component of the McStas package.

The paper is organized as follows: in Section 2 we give a brief description of the ODF, to state clearly the conventions used in this work; in Sections 3 and 4 we present, respectively, the expressions used to compute the differential and total neutron cross sections for coherent elastic scattering by a polycrystalline material; in Section 5 we describe in detail how the method is implemented in the McStas Monte Carlo code; in Section 6 we analyse the effects of the truncation of the Fourier series; in Section 7 we present the results of simulations performed to validate the code; and in Section 8 we discuss, as a first application, an estimation of the SNR of an experiment involving a pressure cell with a sharp texture, comparing it with the SNR associated with a pressure cell of the same characteristics and a uniform texture. Finally, the methods and results are summarized in Section 9. Some details of the computations and other useful information are provided in the appendices.

2. Orientation distribution function

The crystallographic texture of a polycrystalline sample is characterized by its ODF, which gives the relative number of crystal grains that have a particular orientation. The neutron scattering cross section can be computed from the ODF, under some approximations to be discussed in the next section. Let us recall here the basic properties of the ODF, which serves also to fix the notation.

Let $[\{ {\hat x}, {\hat y}, {\hat z}\}]$ be a right-handed orthonormal triad defining a reference frame attached to the sample, and {a, b, c} a system of three independent crystal lattice vectors that generate the whole lattice, oriented in some fixed specified way with respect to the sample frame. The vectors, G, of the reciprocal lattice are determined by the Miller indices hkl through

$[{\bf G} = h{{2\pi} \over {v_0}}{\bf b}\times{\bf c} + k{{2\pi} \over {v_0}}{\bf c}\times{\bf a} + l{{2\pi} \over {v_0}}{\bf a}\times{\bf b}, \eqno(1)]$

where v₀ = a · (b × c) is the volume of the crystal unit cell.

A polycrystal is a material composed of crystalline grains with different orientations at different sample points. The orientation of a grain at a point r in the sample is described by a rotation g(r) ∈ SO(3) (the three-dimensional rotation group), so that the crystal orientation at a point r is given by the triad

$[\{g({\bf r})\,{\bf a}, g({\bf r})\,{\bf b}, g({\bf r})\,{\bf c}\}, \eqno(2)]$

and the vectors of the corresponding reciprocal lattice are given by g(r)G.

The ODF of the polycrystal is a real function $[f: {\rm SO}(3)\rightarrow {\bb R}]$ that gives the volume fraction of grains having an orientation with respect to the sample determined by the rotation g (Bunge, 1993 ). The ODF satisfies f(g) ≥ 0 and is normalized so that

$[\textstyle\int\limits_{{\rm SO}(3)} \, {\rm d}g \ f(g) = 1, \eqno(3)]$

where dg is the Haar (invariant) measure on SO(3), normalized so that

$[\textstyle\int\limits_{{\rm SO}(3)} \, {\rm d}g = 1. \eqno(4)]$

A rotation g can be expressed in terms of the three Euler angles (α, β, γ) as

$[g = g_{\hat z} (\alpha) \, g_{\hat y} (\beta) \, g_{\hat z} (\gamma), \eqno(5)]$

where $[g_{\hat n} (\varphi)]$ denotes the rotation by an angle φ about the $[{\hat n}]$ axis. Note that the Euler angles are defined here in terms of rotations about the fixed sample axes $[\{ {\hat x}, {\hat y}, {\hat z}\}]$ , and α and γ take values in [0, 2π] and β in [0, π]. In terms of the Euler angles, the invariant measure has the form

$[{\rm d}g = {{1} \over {8\pi^2}} \sin\beta \, {\rm d}\beta \, {\rm d}\alpha \, {\rm d}\gamma. \eqno(6)]$

The ODF is the key point of the present work, as it uniquely determines the neutron scattering cross sections in a polycrystalline material, within reasonable assumptions (see next section). In neutron and X-ray diffraction experiments, the ODF is not directly measurable and has to be computed from measurements of related quantities like pole figures. The mathematical problem of extracting the ODF from pole figure measurements is called the pole figure inversion problem and was first addressed in the pioneering work of Bunge (1965 ) and Roe (1965 ). Since then, several methods have been proposed and perfected by several authors (Pospiech & Jura, 1974 ; Jura et al., 1974 , 1976 ; Matthies & Pospiech, 1980 ; Pospiech et al., 1981 ; Houtte, 1983 ; Imhof, 1983 ; Pawlik, 1986 ; Schaeben, 1988 ; Matthies, 1988 ; Helming & Eschner, 1990 ; Houtte, 1991 ; Vadon & Heizmann, 1991 ; van den Boogaart et al., 2007 ; Bernier et al., 2006 ; Hielscher & Schaeben, 2008 ).

The ODF can be expanded in a generalized Fourier series as (Bunge, 1993)

$[f(g) = \textstyle\sum\limits_{l = 0}^\infty \sum\limits_{m = -l}^l \sum\limits_{n = -l}^l C_l^{mn} D_{mn}^l (g), \eqno(7)]$

where D_l^mn(g) are the Wigner D matrices and

$[C_l^{mn} = (2l+1) \textstyle\int\limits_{{\rm SO}(3)} \, {\rm d}g \, D_{mn}^{l^*}(g) \, f(g). \eqno(8)]$

The star superscript stands for complex conjugation. For conciseness, here we call C_l^mn the Fourier coefficients and equation (7) the Fourier series of the ODF, although it is an abuse of language. The relation

$[C_l^{mn^*} = (-1)^{m-n}C_l^{-m-n} \eqno(9)]$

holds by virtue of the reality of the ODF. In terms of the Euler angles, the Wigner matrices are given by

$[D^l_{mn}(\alpha,\beta,\gamma) = \exp{(-i m \alpha)} \, d^{\,l}_{mn} (\cos\beta) \, \exp{(-i n \gamma)}, \eqno(10)]$

where d^l_mn(x) are the Wigner d functions, an explicit expression of which is given in Appendix A. Given an ODF measured on a discrete mesh of SO(3), its Fourier coefficients can be computed with texture analysis software, such as MTEX (Hielscher & Schaeben, 2008).

The Fourier expansion of the ODF is currently used in some Rietveld refinement programs that deal with crystallographic texture, for instance MAUD (Lutterotti et al., 1997 , 1999 ; Wenk et al., 2010 ).

3. Neutron scattering differential cross section

Let us obtain the coherent elastic scattering differential cross section of a neutron propagating through a polycrystalline material. We use the following notation: N_c is the number of unit cells in a crystal, v₀ the volume of the unit cell, G a reciprocal-lattice vector attached to the fixed crystal frame {a, b, c} and F_G the corresponding structure factor. The wavevectors of the incident and scattered neutrons are k and k′, respectively, and the scattering vector is q = k − k′. We deal only with elastic scattering, so that k′ = k.

The coherent elastic differential cross section for the scattering by a perfect single crystal, small enough that the kinematical approximation (disregarding primary extinction) holds, is given by (Squires, 1996 )

$[\left ( {{{\rm d}\sigma} \over {{\rm d}\Omega_{{\bf k}^\prime}}} \right )_{{\bf k} \rightarrow {\bf k}^\prime}^{\rm el,coh} = N_{\rm c} {{(2\pi)^3} \over {v_0}} \, \sum_{\bf G} \left | F_{\bf G} \right |^2 \delta ({\bf k} - {\bf k}^\prime - {\bf G}) . \eqno(11)]$

In a polycrystal there is no interference between the scattering produced by different grains, since they are very large in comparison with the neutron wavelength and highly disoriented, and thus the cross section is merely the sum of the cross sections due to the individual grains (Sears, 1989 ). Furthermore, the grains can be considered as perfect single crystals, since the effect of mosaicity is completely masked by the effect of orientation disorder and can in principle be neglected. Taking all this into account, and given that the number of unit cells with orientation determined by g is N_cf(g)dg, the cross section can be written as

$[\left ( {{{\rm d}\sigma} \over {{\rm d}\Omega_{{\bf k}^\prime}}} \right )_{{\bf k} \rightarrow {\bf k}^\prime}^{\rm el,coh} = \int\limits_{{\rm SO}(3)} {\rm d}g \, f(g) \, N _{\rm c} {{(2\pi)^3} \over {v_0}} \, \sum_{\bf G} \left | F_{\bf G} \right |^2 \delta ({\bf k} - {\bf k}^\prime - g{\bf G}) . \eqno(12)]$

The integral over g can be performed before the sum over G, and thus we have to compute

$[I_{\bf G} ({\bf q}) = \textstyle\int\limits_{{\rm SO}(3)} {\rm d}g \, f(g) \, \delta ({\bf q} - g{\bf G}) . \eqno(13)]$

An explicit expression for I_G(q) can be obtained by using the Fourier expansion of the ODF given by equation (7). Details of the computations are given in Appendix B. The result is the following expression for the differential cross section:

$[\left ( {{{\rm d}\sigma} \over {{\rm d}\Omega_{{\bf k}^\prime}}} \right )_{{\bf k} \rightarrow {\bf k}^\prime}^{\rm el,coh} = N_{\rm c} {{(2\pi)^3} \over {v_0}} \, \sum_{\bf G} {{\left | F_{\bf G} \right |^2} \over {G^2}} \, {{1} \over {4\pi}} \, \delta (q - G) \, \Theta ({\hat G}, {\hat q}) . \eqno(14)]$

Here, $[{\hat q}]$ = q/q is the unit vector along the scattering vector direction, $[{\hat G}]$ = G/G and

$[\Theta( {\hat G}, {\hat q}) = \sum_{l = 0}^\infty {{4\pi} \over {2l + 1}} \sum_{m = -l}^l \sum_{n = -l}^l C_l^{nm} {\rm Y}_l^n ({\hat G}) \, {\rm Y}_l^{m^*} ({\hat q}) , \eqno(15)]$

with $[{\rm Y}_l^m ({\hat n})]$ being the spherical harmonics evaluated at the point $[{\hat n}]$ on the unit sphere. Some properties of the spherical harmonics and the conventions used in this work are summarized in Appendix A.

For a uniform ODF (a `powder') $[\Theta ({\hat G}, {\hat q})]$ = 1, since the only non-vanishing Fourier coefficient is C₀⁰⁰ = 1, and the well known expression for the scattering by a powder is recovered (Sears, 1989). Note that $[\Theta ({\hat G}, {\hat q})]$ is proportional to the corresponding pole function (Bunge, 1993): they differ only by the 4π factor entering equation (15). We prefer this normalization because in this way the $[\Theta ({\hat G}, {\hat q})]$ factor is the modulation in $[{\hat q}]$ of the powder scattering cross section originating from the texture.

Formulas (14) and (15) for the differential scattering cross section were used previously in models for Rietveld refinement programs for textured polycrystals (Popa, 1992 ).

4. Total cross section

The total elastic coherent cross section is obtained by integrating the differential cross section over the direction of the scattered neutron, $[{\hat k^\prime}]$ :

$[\sigma_{\rm el}^{\rm coh} ({\bf k}) = \int {\rm d}\Omega_{\hat k^\prime} \left ( {{{\rm d}\sigma} \over {{\rm d}\Omega_{{\bf k}^\prime}}} \right )_{{\bf k} \rightarrow {\bf k}^\prime}^{\rm el,coh} . \eqno(16)]$

For a textured polycrystal the differential cross section is given by (12). Hence, equation (16) actually involves a double integration over dg and $[{\rm d}\Omega_{\hat k^\prime}]$ . To get an explicit expression in terms of the Fourier coefficients of the ODF we find it convenient to perform the integral over $[{\rm d}\Omega_{\hat k^\prime}]$ first, and then the integral over g, instead of using equation (14). Details of the computations are given in Appendix C. The resulting expression is

$[\sigma_{\rm el}^{\rm coh} ({\bf k}) = N_{\rm c} {{(2\pi)^3} \over {v_0k^3}} \sum_{\bf G} \left | F_{\bf G} \right |^2 {{k} \over {2G}} \, {\rm H} \left ( 1 -G/2k \right ) \, \Upsilon ({\bf G},{\bf k}) , \eqno(17)]$

where H(x) is the Heaviside step function, which is 0 for x < 0 and 1 for x > 0, and

$[\Upsilon ({\bf k},{\bf G}) = \sum_{l = 0}^\infty {{4\pi} \over {2l+1}} \, {\rm P}_l \left ( G/2k \right ) \sum_{m = -l}^l \sum_{n = -l}^l C_l^{nm} \, {\rm Y}_l^n ({\hat G}) \, {\rm Y}_l^{m^*} ({\hat k}) , \eqno(18)]$

where P_l(x) is the Legendre polynomial of order l. Again, for a uniform ODF ϒ(G, k) = 1, and the well known total elastic coherent cross section for the scattering by a powder is recovered (Sears, 1989). To our knowledge, the above expressions for the total cross section in terms of the Fourier coefficients of the ODF [equations (17) and (18)] have not been derived before, although a similar expression for the angular distribution function is given in the book by Bunge (1993) when deriving an expression for the inverse pole figure in terms of the series expansion.

Expressions (17) and (18) can be very useful for analysing transmission experiments involving polycrystalline materials. In particular, in neutron imaging experiments with energy resolution, the elastic coherent term contributes to the appearance of the Bragg edges. The position and shape of these edges will depend on the spacing between the diffraction planes of the grains and on the crystallographic texture, respectively. From equation (17), it is clear that the edge for a particular plane G starts to contribute to the total cross section when the Heaviside function becomes nonzero, i.e. when G = 2k. Thus, in principle, provided the instrument has sufficient energy resolution, the position of the edge will serve to determine the state of strain of those grains whose reciprocal vector G is parallel to the direction of incidence. The shape of the Bragg edge as a function of the incident energy is controlled by the product of three factors: the square of the structure factor, |F_G|², the term k/2G and ϒ given by equation (18). This last term depends on the direction of the incident beam and the scattering plane G, and carries all the information regarding the crystallographic texture through the Fourier coefficients, C_l^mn. In principle, from a mathematical point of view, the summation over l, m and n in equation (18) prevents the possibility of reconstructing the full ODF of a material from a single transmission experiment, even if it is done with energy resolution. However, it is also clear that, from the combined analysis of a set of transmission experiments with different k, some of the C_l^mn can be approximated by inverting equation (18). This can be useful to obtain from imaging experiments integrated quantities that depend on texture, such as Kearns factors for hexagonal crystals (Kearns, 2001 ) or average elastic constants, which only depend on C_l^mn with low l.

The total cross section computed from pole figures has been used to analyse transmission experiment data (Santisteban et al., 2006; Malamud et al., 2014), but the Fourier expansion presented here may have some advantages. First, for neutron wavelengths of the order of 2 Å or less, several planes contribute to the total cross section, as indicated by the Heaviside function of equation (17). In these cases, integration of the pole figure demands the evaluation of several pole figures, while in equation (17) the contributions of all planes are obtained from the same Fourier coefficients. Second, a good description of the total cross section can be obtained using only a few terms in the case of materials with a soft texture, improving computing times. As a by-product, equation (18) can be used as the basis for ODF inversion problems from transmission experiments. This will be studied in more detail in future work.

Finally, it is worth stressing that both expressions for the differential and total cross section [equations (14) and (17)] are valid for any crystal and sample symmetry. All this information is conveyed by the Fourier coefficients.

5. Code implementation

We have developed code to simulate the scattering of thermal neutrons by a polycrystal using expressions (14) and (17). In order to make it available to the community, we have implemented it in the widely used McStas software package (Lefmann & Nielsen, 1999; Willendrup et al., 2004 , 2019 ). The McStas component, called Texture.comp, uses the Union development of McStas (Bertelsen, 2017 ), which is very convenient as it allows the separation of physical processes and geometry.

The strategy for simulating the scattering of a neutron of wavevector k in a material comprises three steps: (i) sampling the neutron free path to get an interaction point; (ii) sampling the interaction process according to a probability proportional to the scattering cross section of the process; and (iii) sampling the wavevector k′ of the outgoing neutron as determined by the differential cross section corresponding to the selected interaction process.

In a homogeneous material, the free path, ξ, is distributed according to an exponential function, $[P_{\rm nfp} (\xi)]$ = $[\mu\exp(-\mu\xi)]$ , where μ is the linear attenuation coefficient (or macroscopic cross section), given by

$[\mu({\bf k}) = {{\rho} \over {A}} \, {{1} \over {N_{\rm c}}} \, \sigma_{\rm tot} ({\bf k}) . \eqno(19)]$

Here, ρ is the density of the material and A the mass contained in the crystal unit cell, so that ρ/A is the number of unit cells per unit volume, and σ_tot/N_c is the total cross section per unit cell. The density can be written as ρ = pA/v₀, where p ∈ [0, 1] is the packing factor, which can be used instead of ρ/A. Thus, for step (i) only μ is required, which, in general, depends on k.

The Monte Carlo simulation requires that the polycrystal be statistically homogeneous (i.e. homogeneous after averaging over grain disorder). If it is not, it has to be divided into statistically homogeneous pieces. Moreover, only the value of μ averaged over the grain disorder enters the free path distribution, P_nfp(ξ). This is one further approximation that amounts to neglect of the spatial correlations of the grain orientations. Note, however, that this approximation is not specific to the polycrystals with non-uniform texture considered here: it is also used for the simulation of powder samples in the current Monte Carlo codes, although spatial correlations in powders are not expected to be very important.

The linear attenuation coefficient receives additive contributions from the different interaction processes available to the neutron (incoherent elastic, coherent elastic, inelastic etc.). Step (ii) selects the interaction process according to the relative probabilities given by the fractional contribution of each process to μ. Thus, for step (ii) only the relative contributions to μ of the available processes are necessary. At step (iii), $[{\bf k}^\prime]$ is sampled according to the differential cross section of the interaction process selected at step (ii). The strategy for this sampling depends strongly on the form of the corresponding differential cross section.

In the Union development of McStas, the geometry and the interaction processes are separated into different components, and multiple scattering is taken into account automatically by the union master component, which calls the functions of the components that deal with geometry to perform the ray tracing, and the functions of the components that deal with the interaction processes to sample the free path, the interaction process and $[{\bf k}^\prime]$ . Two functions provide the interface of an interaction process component, such as Texture.comp, with the McStas union master. One receives k as input and returns the contribution of the interaction process to μ. The other one, which is called if and only if the interaction process described by the component is selected by the master at step (ii), again receives k as input and returns $[{\bf k}^\prime]$ .

Let us describe Texture.comp in some detail. It has to compute the functions ϒ and Θ, for which a cut-off, l_max, on l has to be used, so that the sum over l runs from 0 to l_max. Although the number of terms in the sum is (l_max + 1)², the computation speed does not depend crucially on l_max, since the main ingredients necessary to obtain Θ and ϒ are precomputed in dense two-dimensional grids of $[{\hat k}]$ and $[{\hat q}]$ and interpolated as needed in the course of the simulation. This is one of several optimization strategies implemented in the code.

The evaluation of ϒ and Θ is computationally expensive and some strategies to improve the efficiency have been developed. All the terms that do not depend on k or q can be precomputed and stored in data structures for use in the simulations. The sums entering equations (15) and (18) can be reordered as follows. First we define

$[V_l^n ({\hat G}) = \textstyle\sum\limits_{m = -l}^l C_l^{mn} \, {\tilde {\rm P}}_l^m (\cos\theta_{\hat G}) \exp{(i m \varphi_{\hat G})} . \eqno(20)]$

Here, in general $[\theta_{\hat n}]$ and $[\varphi_{\hat n}]$ denote the polar coordinates of the unit vector $[{\hat n}]$ in the sample reference system and

$[{\tilde {\rm P}}_l^n (x) = \left [ {{(l - n)!} \over {(l + n)!}} \right ]^{1/2} {\rm P}_l^n(x) , \eqno(21)]$

where P_lⁿ(x) are the associated Legendre functions defined in Appendix A. It is convenient to work with the polar components of V_lⁿ, defined by

$[V_l^n ({\hat G}) = R_l^n ({\hat G}) \, \exp [ i \phi_l^n ({\hat G})] , \eqno(22)]$

where $[R_l^n ({\hat G})]$ and $[\phi_l^n ({\hat G})]$ are the modulus and the argument of the complex number $[V_l^n ({\hat G})]$ , respectively. Then, equations (15) and (18) can be expressed as

$[\Theta ({\hat G}, {\hat q}) = \textstyle\sum\limits_{l = 0}^{l_{\rm max}} U_l ({\hat G}, {\hat q}) \eqno(23)]$

and

$[\Upsilon ({\bf G},{\bf k}) =\textstyle \sum\limits_{l = 0}^{l_{\rm max}} {\rm P}_l \left( G/2k \right ) \, U_l ({\hat G}, {\hat k}) , \eqno(24)]$

where

$[U_l ({\hat G}, {\hat n}) =\textstyle \sum\limits_{m = -l}^l R_l^m ({\hat G}) \, {\tilde {\rm P}}_l^m (\cos\theta_{\hat n}) \, \cos \left [ \phi_l^m ({\hat G}) - m \varphi_{\hat n} \right ] . \eqno(25)]$

To derive the above equations we have used the fact that the cross sections are real numbers. The vanishing of the imaginary parts of Θ and ϒ can be readily proved using relation (9) and was used as a test for the code, since an expression similar to (25), with $[\cos[\ldots]]$ replaced by $[\sin[\ldots]]$ , has to vanish.

The computation of Θ and ϒ is too expensive if l_max is large, even using expressions (23) and (24), with precomputed $[R_l^m({\bf G})]$ and $[\phi_l^m({\bf G})]$ . As mentioned above, to speed up the computations $[U_l({\bf G}, {\hat n})]$ is precomputed, for each G and l, on a dense two-dimensional grid of $[\cos\theta_{\hat n} \in [-1,1]]$ and $[\varphi_{\hat n} \in [0,2\pi]]$ . The Legendre polynomials entering (24) are also precomputed on a dense grid of [−1, 1]. Thus, the computation time for ϒ scales as l_max, since it amounts to performing the sum over l with the values of the terms obtained by interpolation. In the case of Θ, the whole sum over l can be precomputed, and thus its computation time is independent of l_max. That is, $[\Theta ({\hat G}, {\hat q})]$ is precomputed, for each G, on a dense two-dimensional grid of $[\cos \theta_{\hat q}]$ and $[\varphi_{\hat q}]$ .

The contribution to μ of the coherent elastic scattering by the textured polycrystal, denoted here by μ_coh, is given by μ_coh = $[( p/v_0 N_{\rm c} ) \sigma_{\rm el}^{\rm coh}]$ , and thus it is numerically computed from equations (17) and (24). The union master component uses it to obtain the interaction point and to sample the interaction process. If the coherent scattering by the polycrystal is selected, Texture.comp has to sample the value of the scattered wavevector, $[{\bf k}^\prime]$ . This sampling is explained in what follows.

First, note that, according to equation (17), the probability that the scattering is due to the set of lattice planes perpendicular to G is

$[{\cal P}_{\rm S} ({\bf G}) = {{1} \over {N_{\rm S}}} \left | F_{\bf G} \right |^2 {{k} \over {2G}} \, {\rm H} \left ( 1 - G/2k \right ) \, \Upsilon ({\bf G}, {\bf k}) \eqno(26)]$

where

$[N_{\rm S} = \sum_{\bf G} \left | F_{\bf G} \right |^2 \, {{k} \over {2G}} \, {\rm H} \left ( 1 - G/2k \right ) \, \Upsilon ({\bf G}, {\bf k}) \eqno(27)]$

is the normalization factor. The values of $[{\cal P}_{\rm S} ({\bf G})]$ are computed and stored when calculating the contribution of the coherent elastic scattering to μ and need not be computed again. A vector G is selected according to the probability $[{\cal P}_{\rm S} ({\bf G})]$ . This sampling is standard, since the set of G that satisfy the condition G < 2k imposed by the Heaviside function (the Bragg cut-off) is finite. For the selected G, the delta function in equation (14) determines the scattering angle, θ_s, which is given by

$[\cos \theta_{\rm s} = {\hat k} \cdot {\hat k^\prime} = 1 - G^2/2k^2 . \eqno(28)]$

Thus, as is well known, $[{\hat k^\prime}]$ lies on the surface of a cone whose axis is given by $[ {\hat k}]$ and whose angle is θ_s (the Debye–Scherrer cone). It only remains to sample the azimuthal angle, $[\varphi^\prime \in [0, 2\pi]]$ , around the cone axis. Introducing two unit vectors $[{\hat t}_1]$ and $[{\hat t}_2]$ so that $[\{ {\hat t}_1, {\hat t}_2, {\hat k} \}]$ forms a right-handed orthonormal triad, we have

$[{\hat k^\prime} = \sin \theta_{\rm s} \cos \varphi^\prime\,{\hat t}_1 + \sin \theta_{\rm s} \sin \varphi^\prime \,{\hat t}_2 + \cos \theta_{\rm s}\, {\hat k} \eqno(29)]$

and

$[{\hat q} = {{k} \over {G}} \left [ -\sin \theta_{\rm s} \cos \varphi^\prime \,{\hat t}_1 - \sin \theta_{\rm s} \sin \varphi^\prime\, {\hat t}_2 + (1 - \cos \theta_{\rm s}) \, {\hat k} \right] . \eqno(30)]$

Let us denote by $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ the probability density of $[\varphi^\prime]$ . Note that this probability density gives the modulation of intensity of the diffracted beam along the Debye–Scherrer cones. According to equation (14), it is obtained, except for a normalization factor, by substituting the above expression for $[{\hat q}]$ into $[\Theta ({\hat G}, {\hat q})]$ , which is computed numerically from equation (23). Hence, given G and k, the first step is to obtain $[\cos \theta_{\rm s}]$ , $[\sin \theta_{\rm s}]$ , and the two vectors $[{\hat t}_1]$ and $[{\hat t}_2]$ . Then, $[\varphi^\prime]$ is sampled according to its probability distribution, for which a simple rejection method is convenient. However, rejection methods need an upper bound for the probability density maximum, and their efficiency is worse the higher the upper bound. In our case, global upper bounds can be obtained from equation (23) but, although they work reasonably well in most instances, they are so bad in some cases that the rejection method becomes useless. The solution, although not very efficient, is to compute $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ in a sufficiently dense grid in [0, 2π] using equation (23). Its maximum is obtained from the discrete values. This computation has to be performed each time $[\varphi^\prime]$ is sampled. The simple rejection method is thus straightforward and works as follows. An angle $[\varphi^\prime]$ is uniformly selected in [0, 2π], and $[{\hat q}]$ is computed according to equation (30). Then $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ is obtained by linear interpolation on the grid. The ratio of this probability density to the maximum probability density is compared with a random number selected uniformly on [0, 1]. If the ratio is smaller than the random number, the value of $[\varphi^\prime]$ is accepted, $[{\hat k^\prime}]$ is computed from equation (29) and $[{\bf k}^\prime]$ = $[k {\hat k^\prime}]$ . Otherwise, another value of $[\varphi^\prime]$ is selected uniformly in [0, 2π] and the process is repeated.

Clearly, the sampling of $[\varphi^\prime]$ is the bottleneck in the simulation. An improvement in the sampling procedure might dramatically increase the code efficiency. The main problem is that, although we can compute $[\Theta ({\hat G}, {\hat q})]$ with reasonably efficiency (by interpolating precomputed values), this is not enough. For an importance sampling method, like rejection, we need the maximum in $[\varphi^\prime]$ of $[\Theta ({\hat G}, {\hat q})]$ , with $[{\hat q}]$ given by equation (30), which depends on the wavevector of the incident neutron, k. To estimate the maximum, $[\Theta ({\hat G}, {\hat q})]$ has to be evaluated many times. Alternatively, we might consider using the weight factor transformation,¹ frequently used in McStas (Willendrup et al., 2019, 2018 ). A similar problem arises in this case, however, since the weight is given by the normalized probability density, $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ , which is proportional to $[\Theta ({\hat G}, {\hat q})]$ , the proportionality factor being the inverse of the integral of $[\Theta ({\hat G}, {\hat q})]$ over $[\varphi^\prime]$ . To compute the normalization factor, which again depends on k, one has to evaluate $[\Theta ({\hat G}, {\hat q})]$ many times. A simple brute-force possibility would be to precompute either the maximum in $[\varphi^\prime]$ of $[\Theta ({\hat G}, {\hat q})]$ or the normalization factor of $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ on a three-dimensional grid in k space. However, the precomputation would be very time consuming and the three-dimensional table very large, and we prefer to avoid interpolation on three-dimensional grids. Hence, we discarded this possibility. In any case, it is clear that there is ample room for improvement at this point.

The need to compute $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ on the $[\varphi^\prime]$ grid each time that $[\varphi^\prime]$ is sampled would make the simulations unfeasible if Θ and ϒ_l had to be computed from expressions (23) and (25). The precomputation of Θ is crucial. In comparison, the use of the precomputed U_l in the evaluation of ϒ is less important: it greatly improves the efficiency of the simulation, but simulations would still be feasible without it.

The efficiency of the simulation can be greatly improved by using variance reduction techniques (stratified sampling) provided by the McStas kernel, especially the use of the SPLIT keyword (Willendrup et al., 2019, 2018). Using the SPLIT n keyword, each incoming neutron is reused n times. The values of μ_coh and $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ computed on the grid are saved and reused when another identical neutron enters the component.

Summarizing, the user has to provide the McStas Texture.comp component with the necessary information through eight input parameters:

(a) The paths to three files: one which contains the coordinates of the crystal reference frame, $[\{ {\hat a}, {\hat b}, {\hat c} \}]$ , in the sample frame; another one which contains the crystallographic information in Lazy/ICSD format (Yvon et al., 1977 ); and a third which contains the Fourier coefficients of the ODF.

(b) Four integer numbers: the cut-off l_max; the sizes in each dimension of the 2D grid in $[(\cos \theta_{\hat n}, \varphi_{\hat n})]$ space where $[\Theta ({\hat G}, {\hat q})]$ and $[U_l ({\hat G}, {\hat k})]$ are precomputed, n_ct × n_φ; and the size of the grid used to sample $[\varphi^\prime]$ , $[n_{\varphi^\prime}]$ .

The program obtains the Miller indices and the corresponding structure factors from the crystallographic file. The user has to guarantee consistency between the sample frame, the crystal reference frame, the crystallographic information and the Fourier coefficients of the ODF.

6. Cut-off effects

The feasibility of Monte Carlo simulations using the method proposed in this work relies on the truncation of the expansion, restricting the sum in l to l ≤ l_max, with l_max sufficiently small. To investigate the cut-off effects, we considered the textures of a Zircaloy-4 plate and a Zr–2.5Nb pressure tube that were obtained experimentally by Malamud et al. (2018 ). The coefficients C^mn_l are computed from the ODFs reported in this reference, using the MTEX software. Fig. 1 displays C_l = $[\sum_{mn} |C_l^{mn}|/(2l + 1)]$ as a function of l for the two materials. The Zr–2.5Nb pressure tube has a sharper texture than the Zircaloy-4 plate, which is reflected in the slower vanishing of the Fourier coefficients as l → ∞.

Figure 1
Fourier coefficients of (a) the Zircaloy-4 plate ODF and (b) the Zr–2.5Nb pressure tube ODF. What is actually plotted is C_l = $[\sum_{m,n}|C_{mn}^l|/(2l+1)]$ .

The dependence on l_max of the pole figures $[P_{\bf G} ({\hat r})]$ associated with the lattice planes perpendicular to the reciprocal-lattice vector G is obtained by truncating at l_max the sum in l in the well known relation (Bunge, 1993)

$[P_{\bf G} ({\hat r}) = \sum_{l = 0}^{\infty} \sum_{m = -l}^l \sum_{n = -l}^l {{C_l^{mn}} \over {2l+1}} \, {\rm Y}_l^n ({\hat G}) \, {\rm Y}_l^{m^*} ({\hat r}) , \eqno(31)]$

where $[{\hat r}]$ represents a direction relative to the sample reference frame.

The cut-off effects on pole figures in the Zircaloy-4 plate and Zr–2.5Nb pressure tube are displayed in Figs. 2 and 3, respectively. The columns, from left to right, correspond to the $[(10{\overline 1}0)]$ , (0002) and $[(11{\overline 2}0)]$ crystal planes, respectively. The top panels display the pole figures computed directly from the experimental ODF [cf. Figs. 4 and 6 of Malamud et al. (2018)]. The lower panels display the absolute difference between the pole figures shown in the top panels and the pole figures computed using equation (31), with cut-offs l_max of 30, 35 and 20 in the Zircaloy-4 case, and of 40, 35 and 30 in the Zr–2.5Nb pressure tube case. The pointwise convergence of the Fourier expansion can be appreciated by looking at the scale set by `Max' in the figures.

Figure 2
Pole figures of the Zircaloy-4 plate corresponding to the crystal planes (left) $[(10{\overline 1}0)]$ , (middle) (0002) and (right) $[(11{\overline 2}0)]$ . In the first row they are computed from the experimental ODF of Malamud et al. (2018

). The second, third and fourth rows display the differences between the pole figures shown in the first row and those computed using equation (31)

with cut-off l_max = 30, 25 and 20, respectively, and Fourier coefficients obtained from the ODF of Malamud et al. (2018

Figure 3
The same as Fig. 2

but for the Zr–2.5Nb pressure tube. The second, third and fourth rows correspond to l_max = 40, 35 and 30, respectively.

The dependence of the cross sections on l_max is also interesting. Fig. 4 displays, as a function of l_max, the contribution of different hkl planes to the total scattering cross section for a neutron of λ = 3.1 Å propagating along a direction given by polar and azimuthal angles of 80 and 60°, respectively, with respect to the sample reference frame. These angles were chosen arbitrarily and correspond to an impinging direction 80° to the normal direction and 60° to the rolling direction in the case of the Zircaloy-4 plate, and 80° to the radial direction and 60° to the axial direction in the case of the Zr–2.5Nb tube. The left- and right-hand panels correspond to the Zircaloy-4 plate and the Zr–2.5Nb pressure tube, respectively. The total cross section is also displayed (black line). To appreciate the convergence towards the l_max → ∞ limit, the values are normalized by those with the largest cut-off (30 and 40 for the Zircaloy-4 plate and the Zr–2.5Nb pressure tube, respectively). The insets magnify the region of larger l_max. We see that the uncertainties introduced by the finite value of l_max are very small if l_max is large enough.

Figure 4
The contributions of different crystallographic planes to the total cross section of (a) the Zircaloy-4 plate and (b) the Zr–2.5Nb pressure tube, as a function of the cut-off l_max. To appreciate the convergence to the l_max → ∞ limit, they are normalized to the value at the largest l_max.

Fig. 5 displays the probability density, $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ , of $[\varphi^\prime]$ for the scattering of a neutron with wavevector k described in the preceding paragraph by various crystal planes of the Zircaloy-4 plate. The probability density is normalized by its maximum. Remember that this function corresponds to the modulations of intensity around the corresponding Debye–Scherrer cone diffracted by a small sample. Each panel corresponds to a different set of crystal planes, whose Miller indices are shown. The different curves correspond to different values of l_max, displayed in the legend. For l = 0, the curves are constant, as they have to be for a uniform texture, with no preferred orientation. Note that the differences decrease considerably as l_max is increased and are large only for l_max < 10.

Figure 5
The probability density, $[{\cal P}_{\bf G} ({\bf k}, \varphi^\prime)]$ , of $[\varphi^\prime]$ for the Zircaloy-4 plate for different crystallographic planes and a fixed k (see main text), for the values of l_max displayed in the legends.

The analogous data for the Zr–2.5Nb pressure tube are displayed in Fig. 6. While for the Zircaloy-4 plate the curves converge for l > 30, for the Zr–2.5Nb texture convergence occurs for l > 40. This is consistent with the higher texture sharpness of the Zr–2.5Nb tube, as discussed by Malamud et al. (2018) and observed in Figs. 2 and 3. In the Zr–2.5Nb case there are appreciable differences for l_max ≤ 15, but for l_max ≥ 20 the differences are small. Note that for l_max ≤ 15 the probability density even becomes negative in some regions. In these regions, however, the true probability is small. The code deals with this problem just by replacing the negative probabilities by zero. This is reasonable since, for instance, some useful although not accurate estimation of the background generated by a Zr–2.5Nb pressure tube might be obtained by Monte Carlo simulations using l_max = 15. Nevertheless, it is advisable to increase the value of l_max since, due to the optimization implemented in the code, this will not significantly affect the efficiency of the simulations.

Figure 6
The same as Fig. 5

but for the Zr–2.5Nb pressure tube.

7. Code validation

To validate the code we performed Monte Carlo simulations to obtain the pole figures of the Zircaloy-4 plate and the Zr–2.5Nb pressure tube, and compared the results with the exact pole figures displayed in the upper panels of Figs. 2 and 3. By exact we mean that these are the pole figures that correspond to the set of Fourier coefficients used in this work, which obviously suffer from the uncertainties associated with their experimental and theoretical determination. The Monte Carlo simulations were performed with the same Fourier coefficients and cut-off, and therefore have to reproduce them with high accuracy. All simulations described in this paper (in this and the next section) were performed by precomputing $[U_l ({\hat G}, {\hat n})]$ and $[\Theta ({\hat G}, {\hat q})]$ on a 201 × 601 uniform grid in $[(\cos \theta, \varphi)]$ space. The grid could probably be made coarser with no great loss of accuracy, but we have not systematically studied the trade-off between simulation accuracy and grid density.

To avoid systematic uncertainties, the simulations are performed in an almost ideal case: the beam, highly collimated (with negligible divergence) and perfectly monochromatic, with λ = 3.1 Å, is scattered by a small spherical sample of 1 mm radius, and multiple scattering is forbidden. The area of the detector is chosen to be small enough that its influence on the results is negligible. The statistical uncertainties are kept low by simulating a high number of neutron histories. Another source of uncertainty is introduced by the approximations made in the code to optimize the computations (for instance, interpolations of precomputed quantities). Although not convenient in real neutron experiments, in this virtual experiment the Schulz setup (Schulz, 1949 ) is used, as shown in Fig. 7(a): the direction of the incident beam and the position of the detector are chosen so that the scattering vector is always directed along the $[{\hat y}_{\rm L}]$ direction in the laboratory reference frame, given by the orthonormal triad $[\{ {\hat x}_{\rm L}, {\hat y}_{\rm L}, {\hat z}_{\rm L}\}]$ . The sample is initially positioned so that $[{\hat x}]$ = $[{\hat x}_{\rm L}]$ , $[{\hat y}]$ = $[- {\hat z}_{\rm L}]$ and $[{\hat z}]$ = $[{\hat y}_{\rm L}]$ . The vectors of the sample reference frame, $[\{ {\hat x}, {\hat y}, {\hat z}\}]$ , are identified with, respectively, the rolling direction (RD), the transverse direction (TD) and the normal direction (ND) in the Zircaloy-4 plate case, and with the axial direction (AD), the hoop direction (HD) and the radial direction (RD) in the Zr–2.5Nb pressure tube case. Then, the sample is rotated by an angle θ about $[{\hat x}]$ and subsequently by an angle φ about $[{\hat z}]$ , and a Monte Carlo simulation is performed. This process is repeated in steps of 5° in both θ and φ, starting from θ = 0 and φ = 0. Fig. 7(b) displays the projection of the scattering vector onto the pole figure as the two rotations on θ and φ are performed. It is clear that, with the set of rotations proposed, a full coverage of the pole figure is achieved.

Figure 7
(a) A scheme of the experimental setup used for the simulation of the pole figure measurement. The sample coordinates correspond to (x, y, z) = (RD, TD, ND) for the Zircaloy-4 plate and (x, y, z) = (AD, HD, RD) for the Zr–2.5Nb pressure tube. (b) A pole figure scan with the experimental setup shown in panel (a).

The pole figures obtained from the Monte Carlo simulation are displayed in Figs. 8 and 9 for the Zircaloy-4 with l_max = 20 and for the Zr–2.5Nb pressure tube with l_max = 30, respectively. The lower panels show the absolute differences relative to the exact result; the scale of the figures indicates that they are small. To appreciate better the quality of the simulations, Fig. 10 displays several cuts at constant θ of the $[(11{\overline 2}0)]$ pole figure. These curves represent the variation in intensity in the pole figure along the circle centred at z with radius θ, as shown in Fig. 7(b). The symbols are the results of a high-statistics Monte Carlo simulation and the continuous red line is the exact result, obtained with the appropriate truncation of equation (31). The error bars signalling the statistical uncertainties of the simulations, smaller than the symbols, are barely visible. The left- and right-hand panels correspond to, respectively, the Zircaloy-4 plate, with l_max = 20, and the Zr–2.5Nb pressure tube, with l_max = 30.

Figure 8
Pole figures of the Zircaloy-4 plate for the crystal planes (left) $[(10{\overline 1}0)]$ , (middle) (0002) and (right) $[(11{\overline 2}0)]$ from a Monte Carlo simulation with l_max = 20. The bottom panels display the differences relative to the exact result given by equation (31)

Figure 9
The same as Fig. 8

but for the Zr–2.5Nb pressure tube with l_max = 30.

Figure 10
Cuts at constant θ of the $[(11{\overline 2}0)]$ pole figure of (a) a Zircaloy-4 plate computed with l_max = 20 and (b) a Zr–2.5Nb pressure tube with l_max = 30. The values of θ are displayed in the legends. The points are the results of a Monte Carlo simulation and the continuous red line the exact result obtained from equation (31)

with the corresponding cut-off.

We also performed simulations to validate the implementation of the linear attenuation coefficient. This implementation, however, is much easier than the implementation of the scattering process, which, as described in Section 5, has to sample $[{\bf k}^\prime]$ according to the proper probability distribution. For the linear attenuation coefficient we only had to implement in the McStas code the computation of μ_coh using equations (17) and (18). The McStas union master, which has been validated elsewhere (Bertelsen, 2017), takes care of sampling the interaction point and the interaction process.

The simulations setup is as follows. A very small rectangular sample, with dimensions 0.1 × 0.1 × 1 mm, is irradiated with an almost perfectly collimated beam, with divergence smaller than 1.0 × 10⁻⁴ °, and with a uniformly distributed wavelength, λ, between 2 and 6 Å. Two detectors with wavelength resolution are placed in front of and behind the sample. To avoid uncertainties we force McStas to absorb neutrons that suffer interaction. In this way, the detector behind the sample collects the neutrons that traverse the sample without interaction, while the detector in front of the sample merely counts the number of incident neutrons. If I₀(λ) and I₁(λ) are the intensities recorded by the detectors in front of and behind the sample, respectively, the simulated linear attenuation coefficient is given by -ln(I₁/I₀)/L, where L = 1 mm is the transmitted neutron path length through the sample. The exact value is computed independently from equations (17) and (18) without relying on McStas. Fig. 11 displays the results. The red and green lines are, respectively, the results of the simulation and the exact values. The panels, from left to right and from top to bottom, correspond to, respectively, beams propagating along the hoop, the axial and the radial directions through a sample with the texture of the Zr–2.5Nb pressure tube, and along the Z direction through a Zr sample with uniform ODF. Note the perfect agreement (within the simulation noise) between the simulation and the exact results.

Figure 11
Simulation of neutron transmission through a small sample under ideal conditions. The upper panels and the bottom left-hand panel correspond to a Zr sample with the texture of the Zr–2.5Nb pressure tube reported by Malamud et al. (2018

). In each case the beam propagates along the sample direction displayed in the figure. The bottom right-hand panel corresponds to a Zr sample with uniform ODF. The red lines are the results of the simulation and the green lines the exact result. In the first three panels the tiny differences are due to noise in the simulation results. In the last panel (bottom right) the differences between the simulation and the exact result cannot be appreciated on the scale of the figure.

Note also the large difference between the attenuation coefficient of a textured material and another with a uniform ODF. The results displayed in Fig. 11 are similar to those presented in Fig. 8 of Santisteban et al. (2012). In that work, experimental values for a Zr–2.5Nb pressure tube with similar texture to the one considered in this work were compared with theoretical evaluations obtained using the technique of integrating pole figures to compute the total coherent elastic cross section.

The perfect agreement between the Monte Carlo results and the exact pole figures and attenuation coefficients is a strong indication that the code works properly. In the case of the Zr–2.5Nb pole figures, very small but sizable (larger than three standard deviations of the statistical uncertainty) discrepancies between the Monte Carlo results and the exact values are evident for θ = 0° and θ = 30°. They are caused by the unavoidable systematic effects of the simulation such as, for instance, the size of the detector, which is small but finite, and the systematic approximations made in the algorithm (e.g. the interpolation of precomputed values and the finite size of the grids).

8. Example: signal-to-noise ratio in a simplified model of a pressure cell

As an example, we simulated the SNR in a simple experiment in which a powder sample of Na₂Ca₃Al₂F₁₄ (Courbion & Ferey, 1988 ) is located inside a cylindrical container that simulates a pressure cell, with a wall made of Zr with the Zr–2.5Nb texture.

Note that, generically, the separation of the detector readings into signal and noise components is not universally defined, but depends on the experimental goals: what in one experiment is part of the signal might be considered noise in a different experiment. Since we are not concerned here with a particular experiment, but with the background generated by the instrument, we consider as signal all neutrons scattered only by the sample, and as noise neutrons scattered at least once by the container. Hence, the SNR is defined here as the ratio between the number of neutrons that reach the detectors after having been scattered only by the sample (one or more times) and the number of neutrons that reach the detectors after having been scattered at least once by the container (and perhaps also by the sample). In some experiments, however, neutrons scattered by the sample incoherently or more than once would be considered noise generated by the sample.

Fig. 12 displays the setup. The container is a hollow cylinder with a diameter of 18 mm, height of 150 mm and thickness of 3 mm, so that inside there is an empty cylindrical space of 12 mm in diameter. The sample has cylindrical geometry, 6 mm in diameter and 10 mm in height. The system is irradiated with a neutron beam that has a Gaussian-distributed wavelength of 3 ± 0.015 Å and a Gaussian distribution of divergence with a standard deviation of 0.4°. The beam is limited by a 6.6 × 110 mm slit, a bit larger than the sample, located 30 cm before it. The scattered neutrons are collected by two area detectors with the geometry of cylindrical sectors of 1 m radius, centred at the sample position, which are 4 m high (vertical direction) and cover angular intervals with respect to the beam direction from 5 to 170° and from 190 to 355°, respectively. By symmetry, the intensity collected by both detectors is essentially the same, and we only show the results for the first detector.

Figure 12
A scheme of the simulated experiment to estimate the SNR. A cylindrical sample is located inside a hollow cylinder modelling a pressure cell made of a Zr alloy. Two detectors of cylindrical shape are placed to have almost 2π coverage (only parts of them are shown). Neutrons (black lines) may suffer from multiple scattering by the different components before reaching the detector.

The texture of the Zr–2.5Nb pressure tube is actually the texture of a small cylindrical sector cut from the tube. That means a whole cylinder is composed of many small cylindrical sectors, with the texture of each sector oriented according to the corresponding AD, HD and RD directions. Hence, we cannot simply simulate a whole cylinder with the same texture, relative to the laboratory frame, at any point. Rather, we have to divide the cylinder into small sectors and assign to each sector the texture oriented according to the local AD, HD and RD directions. In practice we divide the cylinder into 24 sectors of 15°, which causes a big increase in the simulation complexity.

Fig. 13 displays the results. Panel (a) shows the total neutron intensity, in arbitrary units, collected by the detectors; panel (b) displays the intensity of the background, i.e. of neutrons that have been scattered at least once by the container (and some of them also by the sample); panel (c) displays the signal, i.e. the intensity of neutrons that have been scattered only by the sample, showing the intersection with the detectors of the corresponding Debye–Scherrer cones; and panel (d) displays the total intensity and its components, signal and background, along the equatorial plane of the detector system, as a function of the angular position. This provides a typical diffractogram. Note that at some points the background is much higher than the signal. The lines seen in panel (b) correspond to the intersection of the detector surface with Debye–Scherrer cones originated by the wall material, which is not at the centre of the detector system. Therefore, several factors contribute to the modulations along the rings: the crystallographic texture, the differences in neutron path length through the various materials, which causes differences in attenuation, and the differences between the solid angles subtended by the detectors and interaction points.

Figure 13
The intensity collected by the detectors in the simulation of an experiment with the simplified model of a pressure cell described in Section 8

, with the cell walls made of Zr with the texture of the Zr–2.5Nb pressure tube reported by Malamud et al. (2018

). (a) Total intensity, (b) intensity of neutrons scattered at least once by the cell walls (background), (c) intensity of neutrons scattered only by the sample (signal) and (d) intensity of neutrons collected by the detectors in the equatorial plane, providing a typical diffractogram: total (black), signal (red) and background (blue).

To analyse the effect of the cell wall texture, the simulation has been repeated considering a Zr wall with uniform ODF. The results are displayed in Fig. 14. Panels (a) and (b) show the intensity of neutrons scattered only by the cell wall, with the Zr–2.5Nb texture and with the uniform ODF, respectively. Fig. 14(a) is essentially indistinguishable from Fig. 13(b), which means that the intensity of neutrons scattered both by the cell walls and by the sample is small. The different modulations of the intensity along the Debye–Scherrer cones in panels (a) and (b) of Fig. 14 are due to the texture, which has a big influence on the background, as seen in panel (c), where the difference between the intensities of neutrons scattered by both types of cell wall is displayed. The results in the equatorial plane, as a function of the angular position, are shown in panel (d).

Figure 14
The intensity collected at the detectors of neutrons scattered only by the cell walls, (a) with the texture of the Zr–2.5Nb pressure tubes and (b) with a uniform texture. The difference is displayed in (c) and the results in the equatorial plane in (d).

The texture has an important influence on the SNR, as expected. The left-hand panel of Fig. 15 displays the ratio SNR₁/SNR₂, where SNR₁ and SNR₂ stand for the SNR of a Zr pressure cell with the Zr–2.5Nb texture and with uniform ODF, respectively. The right-hand panel shows the data along the equatorial plane of the detector system. There are points at which the SNRs differ by a factor higher than five. Note that, although the SNR depends crucially on the sample, which provides the signal, the ratio of SNRs in the case of two different containers is nearly independent of the sample, since the influence of the container on the signal is rather small. Thus, the ratio of SNRs is, to high accuracy, the ratio of the background produced by the containers, which, in turn, depends very little on the sample. This means that Fig. 15(a) can be visualized, to a very good approximation, as the point-to-point ratio of Figs. 14(b) and 14(a).

Figure 15
(a) The ratio of SNR₁/SNR₂ of the SNR of a pressure cell with Zr walls with the texture of the Zr–2.5Nb pressure tubes (SNR₁) and with a uniform texture (SNR₂). (b) A cut of panel (a) along the equatorial plane of the detector system, height = 0.

9. Summary and conclusions

We have developed a method of simulating the transport of thermal neutrons through polycrystalline materials. It is based on the generalized Fourier expansion, in terms of Wigner D matrices, of the orientation distribution function, which leads to an expansion of the differential and total cross sections in terms of the Fourier coefficients. These expansions are suitable for Monte Carlo codes. As expected, the expression for the differential cross section associated with a crystal plane is proportional to the well known analogous expansion of the corresponding pole figure (Bunge, 1993). Although alternative expressions are currently used, to our knowledge the expression for the total cross section given here has not been derived before. In some cases, for instance in Monte Carlo codes, it has advantages over other expressions.

The method has been implemented in a McStas component code called Texture.comp. It has been validated by computing the pole figures of a Zircaloy-4 plate and a Zr–2.5Nb pressure tube through Monte Carlo simulations of an ideal neutron diffraction experiment, where the sample is rotated about two axes, and by simulating a transmission experiment under ideal conditions. As a first application, we estimated the signal-to-noise ratio of a diffraction experiment in which a small sample is placed inside a cylindrical pressure cell made of a Zr alloy with the texture of the Zr–2.5Nb pressure tube obtained by Malamud et al. (2018). To see the effect of texture, the simulation was repeated considering a Zr alloy with uniform texture. We found that texture has a deep impact on the SNR: at some points the two SNRs differ by a factor greater than five.

The computational cost of simulating thermal neutron transport through textured polycrystals is obviously much higher than that through a polycrystal with a uniform ODF. The cost depends strongly on the complexity of the problem. The higher the complexity, the higher the relative cost of the problem with non-trivial texture. For the simplest problem, in which neutrons are scattered only by a small sample, so that multiple scattering is very unlikely, the computing time in the non-trivial texture case is only three and a half times longer than that in the uniform ODF case, if the McStas SPLIT keyword is used heavily. Without using SPLIT, it is 12 times longer. For complex problems the SPLIT keyword is not as effective. For instance, in the problem described in Section 8, which is rather complex since the cylindrical container with non-trivial texture was divided into 24 sectors, the simulations (using SPLIT) with the Zr–2.5Nb texture were 32 times longer than with the uniform ODF. This is, however, not a big problem, given (i) the power of current computation resources, (ii) that the simulations are trivially parallelizable and (iii) that McStas is extremely fast at simulating powder materials.

The generalized Fourier expansion of the ODF is useful only if the texture is sufficiently mild. For very sharp textures, the ODF can be split into a smooth component and some sharp peaks. The smooth part can be simulated with the software presented here and the sharp peaks with the methods proposed by Laliena et al. (2019 ).

The software will be used to compute the background generated by components like pressure cells in neutron scattering instruments, which depends strongly on the texture of the device materials. The final goal is to assist in the design of neutron instruments for extreme conditions by estimating, through Monte Carlo simulations, the SNR in different configurations. The software, however, has a much broader scope, and may also be used for the analysis of experiments involving samples of polycrystalline materials, like pole figures and residual stresses in alloys. Interestingly, the expression for the total cross section can be used to analyse data from transmission experiments with polycrystalline materials (Vicente-Álvarez et al., work in progress).

The component Texture.comp will appear in the next McStas release, so that it will be available to the community. It can also be obtained in advance from the authors upon request.

APPENDIX A

Some properties of SO(3)

The Hermitian infinitesimal generators of SO(3) satisfy the algebra

$[[L_x ,L_y] = i L_z, \quad [L_x, L_z] = -i L_y, \quad [L_y, L_z] = i L_x. \eqno(32)]$

The irreducible representations of SO(3) are characterized by the non-negative integers l associated with the eigenvalues, l(l + 1), of L² = L_x² + L_y² + L_z², and have dimension 2l + 1. An orthonormal basis of the representation space is given by the eigenvectors $[| lm \rangle]$ of L² and L_z, with the integer m bounded by −l ≤ m ≤ l. This basis is uniquely determined by choosing a phase convention. Here we adhere to the Condon & Shortley (1957 ) choice of phases, which is fixed by the relation

$[L_\pm | lm \rangle = [(l \pm m) \, (l \pm m + 1)]^{1/2} | lm \pm 1 \rangle , \eqno(33)]$

where $[L_\pm = L_x \pm {\rm i} L_y]$ .

For each g ∈ SO(3), the Wigner D matrix is defined by

$[D^l_{mn} (g) = \langle lm | R_l(g) | ln \rangle , \eqno(34)]$

where R_l(g) is the operator that implements the action of g in the l representation, which implies the relation

$[D^l_{mn} (g_1 g_2) = \textstyle\sum\limits_{m^\prime = -l}^l D^l_{mm^\prime} (g_1) \, D^l_{m^\prime n} (g_2) . \eqno(35)]$

In addition, the unitarity of the representation implies that the Wigner matrices are unitary and

$[D^l_{mn} (g^{-1}) = D^{l^*}_{nm} (g) . \eqno(36)]$

In terms of Euler angles the Wigner matrices read

$[D^l_{mn} (\alpha, \beta, \gamma) = \langle lm| \!\exp{(-i \alpha L_z)} \exp{(-i \beta L_y)} \exp{(-i \gamma L_z)} | ln \rangle , \eqno(37)]$

so that

$[D^l_{mn} (\alpha, \beta, \gamma) = \exp{(-i m \alpha)} \, d^{\,l}_{mn} (\cos \beta) \exp{(-i n \gamma)} , \eqno(38)]$

where

$[ d^{\,l}_{mn} (\cos \beta) = \langle lm | \!\exp{(-i \beta L_y)} | ln \rangle \eqno(39)]$

is called the Wigner d function. We have the following explicit formula (Galindo & Pascual, 1990 ):

$[\eqalignno{d^{\,l}_{mn} (x) = & \, {{(-1)^{l-n}} \over {2^l}} \left [ {{(l+n)!} \over {(l-n)!(l+m)!(l-m)!}} \right ]^{1/2} \cr & \, \times \left [ {{(1-x)^{m-n}} \over {(1+x)^{m+n}}} \right ]^{1/2} {{{\rm d}^{\,l-n}} \over {{\rm d}x^{l-n}}} \left [ (1-x)^{l-m} (1+x)^{l+m} \right ] . &(40)}]$

Below we will use the following relation (Galindo & Pascual, 1990):

$[D_{0m}^l (0, \beta, \gamma) = \left ( {{4\pi} \over {2l+1}} \right )^{1/2} {\rm Y}_l^{-m} (\beta, \gamma), \eqno(41)]$

where $[{\rm Y}_l^m (\theta, \varphi)]$ are the spherical harmonics, defined by

$[{\rm Y}_l^m (\theta, \varphi) = \left [ {{2l+1} \over {4\pi}} {{(l-m)!} \over {(l+m)!}} \right ]^{1/2} {\rm P}_l^m (\cos \theta) \, \exp{(i m \varphi)} , \eqno(42)]$

with 0 ≤ θ ≤ π and 0 ≤ φ < 2π, and P_l^m (x) are the Legendre associated functions, given by

$[{\rm P}_l^m (x) = (1-x^2)^{m/2} {{(-1)^m} \over {l!2^l}} {{d^{\,l+m}} \over {dx^{l+m}}} \, (x^2 - 1)^l . \eqno(43)]$

Note that the Legendre polynomial P_l(x) of order l is the associated Legendre function of order l and m = 0. The arguments of the spherical harmonics define a unit vector, $[{\hat n}]$ , given by the polar angle θ and the azimuthal angle φ, so that we can use the notation $[{\rm Y}_l^m ({\hat n})]$ . Other useful properties of the Wigner matrices and the representations of SO(3) can be found in many books, for instance in Appendix B of Galindo & Pascual (1990).

APPENDIX B

Computation of the differential cross section

To compute the integral I_G(q) of equation (13) we use the representation of the three-dimensional Dirac delta function in polar coordinates:

$[\delta ({\bf x} - {\bf y}) = {{1} \over {x^2}} \, \delta (x - y) \, \delta_{\rm S} ({\hat x}, {\hat y}) , \eqno(44)]$

where $[\delta_{\rm S} ({\hat x}, {\hat y})]$ is the Dirac delta function on the sphere, which in polar coordinates θ and φ reads

$[\delta_{\rm S} ({\hat x}, {\hat y}) = \delta (\cos \theta_{\hat y} - \cos \theta_{\hat x}) \, \delta_{\rm P} (\varphi_{\hat y} - \varphi_{\hat x}) , \eqno(45)]$

and δ_P(φ) is the periodic Dirac delta function. Then, introducing the Fourier representation of the ODF, f(g), we have

$[I_{\bf G} ({\bf q}) = {{1} \over {G^2}} \delta ( | {\bf k} - {\bf k}^\prime | - G) \sum_{l = 0}^\infty \sum_{m = -l}^l \sum_{n = -l}^l C_l^{mn} I^l_{mn} , \eqno(46)]$

where

$[I^l_{mn} = \textstyle\int\limits_{{\rm SO}(3)} {\rm d}g \, D^l_{mn} (g) \, \delta_{\rm S} ({\hat q}, g {\hat G}) , \eqno(47)]$

with $[{\hat q} = ({\bf k} - {\bf k}^\prime)/ | {\bf k} - {\bf k}^\prime |]$ . Let $[g_{\hat n}]$ denote a rotation that brings a unit vector $[{\hat n}]$ to $[{\hat z}]$ , so that $[g_{\hat n} {\hat n} = {\hat z}]$ . Since the Dirac delta function on the sphere is rotationally invariant, we have

$[\delta_{\rm S} ({\hat q}, g {\hat G}) = \delta_{\rm S} ({\hat z}, g_{\hat q} g g_{\hat G}^{-1} {\hat z}) . \eqno(48)]$

Using the invariance of the measure we have

$[I^l_{mn} = \textstyle\int\limits_{{\rm SO}(3)} {\rm d} g \, D^l_{mn} (g_{\hat q}^{-1} g g_{\hat G}) \, \delta_{\rm S} ({\hat z},g {\hat z}) . \eqno(49)]$

Since the Wigner D functions support a unitary representation of SO(3) [equation (35)], we have

$[I^l_{mn} = \textstyle\sum\limits_{m^\prime} \sum\limits_{n^\prime} D^{l^*}_{m^\prime m} (g_{\hat q}) \, D^l_{n^\prime n} (g_{\hat G}) \int\limits_{{\rm SO}(3)} {\rm d}g \, D^l_{m^\prime n^\prime} (g) \, \delta_{\rm S} ({\hat z}, g {\hat z}) . \eqno(50)]$

The integral can be readily evaluated in terms of the Euler angles, (α, β, γ), that characterize g, since

$[g {\hat z} = \sin \beta \cos \alpha \,{\hat x} + \sin \beta \sin \alpha \,{\hat y} + \cos \beta \,{\hat z} , \eqno(51)]$

so that

$[\delta_{\rm S} ({\hat z}, g {\hat z}) = \delta (\cos \beta - 1) \, \delta_{\rm P} (\alpha) . \eqno(52)]$

Taking into account that $[d^{\,l}_{mn}(1) = \delta_{mn}]$ , we obtain

$[\int\limits_{{\rm SO}(3)} {\rm d}g \, D^l_{m^\prime n^\prime} (g) \, \delta_{\rm S} ({\hat z}, g {\hat z}) = {{1} \over {4\pi}} \delta_{m^\prime 0} \delta_{n^\prime 0}, \eqno(53)]$

and therefore

$[I^l_{mn} = {{1} \over {4\pi}} D^{l^*}_{0m} (g_{\hat q}) \, D^l_{0n} (g_{\hat G}) . \eqno(54)]$

The Euler angles corresponding to a rotation $[g_{\hat n}]$ that brings $[{\hat n}]$ to $[{\hat z}]$ are α = 0, $[\beta = \theta_{\hat n}]$ and $[\gamma = \pi - \varphi_{\hat n}]$ , where $[\theta_{\hat n}]$ and $[\varphi_{\hat n}]$ are the polar coordinates of $[{\hat n}]$ , and therefore, using relation (41), we have

$[D^l_{0m} (g_{\hat n}) = \left ( {{4\pi} \over {2l+1}} \right )^{1/2} {\rm Y}_l^m ({\hat n}) . \eqno(55)]$

Using the above expression, equation (54) becomes

$[I^l_{mn} = {{1} \over {2l+1}} {\rm Y}_l^n ({\hat G}) \, {\rm Y}_l^{m ^*} ({\hat q}) , \eqno(56)]$

and inserting this expression into equation (46) we obtain the desired result.

APPENDIX C

Computation of the total cross section

The total cross section can be obtained by integrating over the solid angle $[{\rm d}\Omega_{{\bf k}^\prime}]$ the differential cross section given by equation (14). However, it is easier to integrate equation (12) over $[{\rm d}\Omega_{{\bf k}^\prime}]$ before performing the integral over dg. Since the scattering is elastic, we have $[{\bf k}^\prime = k{\hat k^\prime}]$ , so that

$[\eqalignno{\delta ({\bf k} - {\bf k}^\prime - g{\bf G}) = & \, \delta ({\bf k} - g{\bf G} - k{\hat k^\prime}) \cr = & \, {{1} \over {k^2}} \delta (|{\bf k} - g{\bf G}| - k) \, \delta_{\rm S} \left ( {\hat k} - {{1} \over {k}} g {\bf G}, {\hat k^\prime} \right ) . &(57)}]$

The integral over $[{\rm d}\Omega_{\hat k^\prime}]$ of the above expression is 1, since the first Dirac delta function forces $[{\hat k} - g {\bf G}/k]$ to lie on the unit sphere. Therefore we have

$[\sigma_{\rm el}^{\rm coh} = N {{(2\pi)^3} \over {v_0k^3}} \sum_{\bf G} |F_{\bf G}|^2 \int {\rm d}g \, f(g) \, k \delta (|{\bf k} - g{\bf G}| - k) . \eqno(58)]$

Let us remember that we defined $[g_{\hat n}]$ as the rotation that brings the unit vector $[{\hat n}]$ to $[{\hat z}]$ , so that $[g_{\hat n} {\hat n} = {\hat z}]$ . Then we have

$[|{\bf k} - g{\bf G}| = |kg_{\hat k}^{-1} {\hat z} - Ggg_{\hat G}^{-1} {\hat z}| = |k {\hat z} - Gg_{\hat k} gg_{\hat G}^{-1} {\hat z}| . \eqno(59)]$

Using the invariance of the measure we can write

$[\textstyle\int \!{\rm d}g \, f(g) \, k \delta (| {\bf k} - g{\bf G}| \!- k) = \!\int \!{\rm d}g \, f(g_{\bf k}^{-1} gg_{\bf G}) \, \delta (|k {\hat z} - Gg {\hat z}| \!- k) , \eqno(60)]$

and, using the Fourier expansion of the ODF and the group properties of the Wigner matrices, we have

$[\eqalignno{& \textstyle\int\limits {\rm d}g \, f(g) \, k\delta (|{\bf k} - g{\bf G}| - k) = \cr & \quad \quad \textstyle\sum\limits_{l = 0}^\infty \sum\limits_{m} \sum\limits_{n} \sum\limits_{m^\prime} \sum\limits_{n^\prime} C_l^{mn} D_{m^\prime m}^{l^*} (g_{\hat k}) \, D_{n^\prime n}^l (g_{\hat G}) \, S_{m^\prime n^\prime}^l (G, k) , &(61)}]$

where

$[S_{m^\prime n^\prime}^l (G, k) = \textstyle\int\limits {\rm d}g \, D_{m^\prime n^\prime}^l (g) \, k\delta (|k {\hat z} - Gg {\hat z}| - k) . \eqno(62)]$

Using equation (51) we have

$[|k {\hat z} - Gg {\hat z}| = \left ( k^2 + G^2 - 2kG \cos \beta \right )^{1/2} , \eqno(63)]$

which is independent of α and γ. Therefore, the integrals over these angles give $[4 \pi^2 \delta_{m^\prime 0} \delta_{n^\prime 0}]$ . Making the change of variable z = $[\cos\beta]$ and taking into account that d₀₀^l(z) = P_l(z), we have

$[\eqalignno{& S_{m^\prime n^\prime}^l (G ,k) = \cr & \delta_{m^\prime 0} \delta_{n^\prime 0} {{1} \over {2}} \int\limits_{-1}^{1} \, {\rm d}z \, {\rm P}_l (z) \, k \delta \left [ \left ( k^2 + G^2 - 2kG \cos \beta \right )^{1/2} - k \right ] . &(64)}]$

The argument of the delta function vanishes for z = G/2k, and its derivative at this point is −G. Then we have

$[S_{m^\prime n^\prime}^l (G,k) = \delta_{m^\prime 0} \delta_{n^\prime 0} {{k} \over {2G}} {\rm P}_l \left ( {{G} \over {2k}} \right ) {\rm H} \left ( 1 - {{G} \over {2k}} \right ) . \eqno(65)]$

Inserting the above result into (61) and using equation (55) we get

$[\int {\rm d}g \, f(g) \, k \,\delta \,(|{\bf k} - g{\bf G}| - k) = {{k} \over {2G}} \Upsilon ({\bf G}, {\bf k}) . \eqno(66)]$

Footnotes

¹We are grateful to one of the referees for bringing the weight factor transformation to our attention. It might have some advantages over importance sampling methods, opening a way for improving the code efficiency, although the problems discussed in the text have to be overcome first.

Acknowledgements

VL thanks Mads Bertelsen for his help with McStas and its Union component set, and E. Lelièvre-Berna for discussions and encouragement.

Funding information

Funding for this research was provided by the European Union's Horizon 2020 (grant No. 654000); Ministerio de Ciencia, Innovación y Universidades (Spain) (grant No. PGC-2018-099024-B- I00); Consejo Superior de Investigaciones Científicas (Spain) (grant No. I-COOP B20319).

References

Bernier, J. V., Miller, M. P. & Boyce, D. E. (2006). J. Appl. Cryst. 39, 697–713. Web of Science CrossRef CAS IUCr Journals Google Scholar
Bertelsen, M. (2017). PhD thesis, The Niels Bohr Institute, Faculty of Science, University of Copenhagen, Denmark. Google Scholar
Boin, M. (2012). J. Appl. Cryst. 45, 603–607. Web of Science CrossRef CAS IUCr Journals Google Scholar
Boogaart, K. G. van den, Hielscher, R., Prestin, J. & Schaeben, H. (2007). J. Comput. Appl. Math. 199, 122–140. Google Scholar
Bunge, H. J. (1965). Z. Metallkd. 56, 872–874. CAS Google Scholar
Bunge, H.-J. (1993). Texture Analysis in Materials Science. Göttingen: Cuvillier Verlag. Google Scholar
Condon, E. U. & Shortley, G. H. (1957). The Theory of Atomic Spectra. Cambridge University Press. Google Scholar
Courbion, G. & Ferey, G. (1988). J. Solid State Chem. 76, 426–431. CrossRef ICSD CAS Web of Science Google Scholar
Dessieux, L. L., Stoica, A. D. & Bingham, P. R. (2018). Rev. Sci. Instrum. 89, 025–103. Google Scholar
Dessieux, L. L., Stoica, A. D., Bingham, P. R., An, K., Frost, M. J. & Bilheux, H. Z. (2019). Nucl. Instrum. Methods Phys. Res. B, 459, 166–178. Google Scholar
Dollase, W. A. (1986). J. Appl. Cryst. 19, 267–272. CrossRef CAS Web of Science IUCr Journals Google Scholar
Farhi, E., Hugouvieux, V., Johnson, M. & Kob, W. (2009). J. Comput. Phys. 228, 5251–5261. Web of Science CrossRef CAS Google Scholar
Galindo, A. & Pascual, P. (1990). Quantum Mechanics I. Berlin, Heidelberg, New York: Springer-Verlag. Google Scholar
Helming, K. & Eschner, T. (1990). Cryst. Res. Technol. 25, K203–K208. CrossRef CAS Web of Science Google Scholar
Hielscher, R. & Schaeben, H. (2008). J. Appl. Cryst. 41, 1024–1037. Web of Science CrossRef CAS IUCr Journals Google Scholar
Houtte, P. V. (1983). Textures Microstruct. 6, 1–19. Google Scholar
Houtte, P. V. (1991). Textures Microstruct. 13, 199–212. Google Scholar
Imhof, J. (1983). Phys. Status Solidi A, 75, K187–K189. CrossRef CAS Google Scholar
Jura, J., Pospiech, J. & Bunge, H. J. (1974). Texture, 1, 201–203. Google Scholar
Jura, J., Pospiech, J. & Bunge, H. J. (1976). Metallurgia, 24, 111–176. Google Scholar
Kearns, J. (2001). J. Nucl. Mater. 299, 171–174. Google Scholar
Kibble, M. G., Laliena, V., Goodway, C. M., Lelièvre-Berna, E., Kamenev, K. V., Klotz, S. & Kirichek, O. (2019). J. Neutron Res. 21, 105–116. Google Scholar
Laliena, V., Vicente Álvarez, M. A. & Campo, J. (2019). J. Neutron Res. 21, 95–104. Google Scholar
Lee, W. & Wang, X. L. (2002). Neutron News, 13(4), 30–34. Google Scholar
Lefmann, K. & Nielsen, K. (1999). Neutron News, 10(3), 20–23. CrossRef Google Scholar
Lin, J. Y. Y., Smith, H. L., Granroth, G. E., Abernathy, D. L., Lumsden, M. D., Winn, B., Aczel, A. A., Aivazis, M. & Fultz, B. (2016). Nucl. Instrum. Methods Phys. Res. A, 810, 86–99. CrossRef CAS Google Scholar
Lutterotti, L., Matthies, S. & Wenk, H.-R. (1999). Proceedings of the Twelfth International Conference on Textures of Materials (ICOTOM-12), edited by J. A. Szpunar, Vol. 1, p. 1599. Ottawa: NRC Research Press. Google Scholar
Lutterotti, L., Matthies, S., Wenk, H. -R., Schultz, A. S. & Richardson, J. W. Jr (1997). J. Appl. Phys. 81, 594–600. Google Scholar
Malamud, F., Riffo, A. M., Alvarez, M. V., Vizcaino, P., Li, M., Liu, X., Vogel, S., Law, M., Sumin, V., Luzin, V., Vasin, R. & Santisteban, J. (2018). J. Nucl. Mater. 510, 524–538. Google Scholar
Malamud, F., Santisteban, J. R., Vicente Alvarez, M. A., Bolmaro, R., Kelleher, J., Kabra, S. & Kockelmann, W. (2014). J. Appl. Cryst. 47, 1337–1354. Web of Science CrossRef CAS IUCr Journals Google Scholar
Matthies, S. (1988). Proceedings of the Eighth International Conference on Textures of Materials (ICOTOM8), Santa Fe, New Mexico, 20–25 September 1987, edited by J. S. Kallend & G. Gottstein, pp. 37–48. Warrendale: The Metallurgical Society. Google Scholar
Matthies, S. & Pospiech, J. (1980). Phys. Status Solidi B, 97, 547–556. Google Scholar
Mikula, P., Luká˘s, P. & Vrána, M. (1997). Physica B, 234–236, 1058–1060. Google Scholar
Pawlik, K. (1986). Phys. Status Solidi, 134, 447–483. Google Scholar
Popa, N. C. (1992). J. Appl. Cryst. 25, 611–616. CrossRef Web of Science IUCr Journals Google Scholar
Pospiech, J. & Jura, J. (1974). Z. Metallkd. 65, 324–330. CAS Google Scholar
Pospiech, J., Ruer, D. & Baro, R. (1981). J. Appl. Cryst. 14, 230–233. CrossRef CAS Web of Science IUCr Journals Google Scholar
Rodríguez-Velamazán, J. A., Campo, J., Rodríguez-Carvajal, J. & Noguera, P. (2011). J. Phys. Conf. Ser. 325, 012010. Google Scholar
Rodríguez-Velamazán, J. A. & Noguera, P. (2011). J. Phys. Conf. Ser. 325, 012022. Google Scholar
Roe, R. (1965). J. Appl. Phys. 36, 2024–2031. CrossRef CAS Web of Science Google Scholar
Santisteban, J. R., Edwards, L. & Stelmukh, V. (2006). Physica B, 385–386, 636–638. Web of Science CrossRef CAS Google Scholar
Santisteban, J. R., Vicente-Alvarez, M. A., Vizcaino, P., Banchik, A. D., Vogel, S. C., Tremsin, A. S., Vallerga, J. V., McPhate, J. B., Lehmann, E. & Kockelmann, W. (2012). J. Nucl. Mater. 425, 218–227. Web of Science CrossRef CAS Google Scholar
Šaroun, J., Kornmeier, J. R., Hofmann, M., Mikula, P. & Vrána, M. (2013). J. Appl. Cryst. 46, 628–638. Google Scholar
Šaroun, J. & Kulda, J. (1997). Physica B, 234–236, 1102–1104. Google Scholar
Sato, H., Kamiyama, T. & Kiyanagi, Y. (2011). Mater. Trans. 52, 1294–1302. Web of Science CrossRef CAS Google Scholar
Schaeben, H. (1988). J. Appl. Phys. 64, 2236–2237. CrossRef Web of Science Google Scholar
Schulz, L. G. (1949). J. Appl. Phys. 20, 1030–1033. CrossRef Web of Science Google Scholar
Sears, V. F. (1989). Neutron Optics. New York: Oxford University Press. Google Scholar
Squires, G. L. (1996). Introduction to the Theory of Thermal Neutron Scattering. New York: Dover Publications Inc. Google Scholar
Stoica, A. D., Popovici, M. & Hubbard, C. R. (2001). J. Appl. Cryst. 34, 343–357. Web of Science CrossRef CAS IUCr Journals Google Scholar
Vadon, A. & Heizmann, J. J. (1991). Textures Microstruct. 14, 37–44. CrossRef Web of Science Google Scholar
Vickery, A., Udby, L., Violini, N., Voigt, J., Deen, P. & Lefmann, K. (2013). J. Phys. Soc. Jpn, 82, SA037. Google Scholar
Vogel, S. (1999). PhD thesis, Christian-Albrechts-Universität, Kiel, Germany. Google Scholar
Wenk, H.-R., Lutterotti, L. & Vogel, S. C. (2010). Powder Diffr. 25, 283–296. Web of Science CrossRef CAS Google Scholar
Willendrup, P., Farhi, E., Knudsen, E., Filges, U., Lefmann, K. & Stein, J. (2018). User and Programmers Guide to the Neutron Ray-Tracing Package McStas. Version 2.5. https://www.mcstas.org/documentation/manual/. Google Scholar
Willendrup, P., Farhi, E. & Lefmann, K. (2004). Physica B, 350, E735–E737. Web of Science CrossRef CAS Google Scholar
Willendrup, P. & Lefmann, K. (2019). J. Neutron Res. https://doi.org/10.3233/JNR-190108. Google Scholar
Yvon, K., Jeitschko, W. & Parthé, E. (1977). J. Appl. Cryst. 10, 73–74. CrossRef IUCr Journals Web of Science Google Scholar
Zhao, J. (2011). Nucl. Instrum. Methods Phys. Res. A, 659, 434–441. Google Scholar
Zsigmond, G., Lieutenant, K. & Mezei, F. (2002). Neutron News, 13(4), 11–14. CrossRef Google Scholar

This is an open-access article distributed under the terms of the Creative Commons Attribution (CC-BY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.

JOURNAL OF
APPLIED
CRYSTALLOGRAPHY

ISSN: 1600-5767

Volume 53| Part 2| April 2020| Pages 512-529

https://doi.org/10.1107/S1600576720002290

Open

access