Monte Carlo simulation of neutron scattering by a textured polycrystal

A method of simulating the diffraction of thermal neutrons by polycrystalline materials, using the generalized Fourier transform of the orientation distribution function, is developed and implemented in a Monte Carlo code in the McStas software package. Simulations to validate the code and a first application are also presented.


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-tonoise 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 pseudostrains 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 . 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(Dessieux et al., , 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(Santisteban et al., , 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.

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 fx x;ŷ y;ẑ zg 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 where v 0 = 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) 2 SO(3) (the three-dimensional rotation group), so that the crystal orientation at a point r is given by the triad fgðrÞ a; gðrÞ b; gðrÞ cg; ð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 : SOð3Þ ! 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 R where dg is the Haar (invariant) measure on SO(3), normalized so that R SOð3Þ dg ¼ 1: A rotation g can be expressed in terms of the three Euler angles (, , ) as g ¼ g^z z ðÞ g^y y ðÞ g^z z ðÞ; where g^n n ð'Þ denotes the rotation by an angle ' about then n axis. Note that the Euler angles are defined here in terms of rotations about the fixed sample axes fx x;ŷ y;ẑ zg, and and take values in [0, 2] and in [0, ]. In terms of the Euler angles, the invariant measure has the form 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 Jura et al., 1974Jura et al., , 1976Matthies & 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) where D mn l ðgÞ are the Wigner D matrices and C mn The star superscript stands for complex conjugation. For conciseness, here we call C mn l the Fourier coefficients and equation (7) the Fourier series of the ODF, although it is an abuse of language. The relation holds by virtue of the reality of the ODF. In terms of the Euler angles, the Wigner matrices are given by 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(Lutterotti et al., , 1999Wenk et al., 2010).

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 0 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 0 , respectively, and the scattering vector is q = k À k 0 . We deal only with elastic scattering, so that k 0 = 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) 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 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 c f(g)dg, the cross section can be written as The integral over g can be performed before the sum over G, and thus we have to compute 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: Here,q q = q/q is the unit vector along the scattering vector direction,Ĝ G = G/G and with Y m l ðn nÞ being the spherical harmonics evaluated at the pointn 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') ÂðĜ G;q qÞ = 1, since the only non-vanishing Fourier coefficient is C 00 0 = 1, and the well known expression for the scattering by a powder is recovered (Sears, 1989). Note that ÂðĜ G;q 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 ÂðĜ G;q qÞ factor is the modulation inq 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).

Total cross section
The total elastic coherent cross section is obtained by integrating the differential cross section over the direction of the scattered neutron,k k 0 : For a textured polycrystal the differential cross section is given by (12). Hence, equation (16) actually involves a double integration over dg and d^k k 0 . To get an explicit expression in terms of the Fourier coefficients of the ODF we find it convenient to perform the integral over d^k k 0 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 where H(x) is the Heaviside step function, which is 0 for x < 0 and 1 for x > 0, and 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 | 2 , 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 mn l . 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 mn l 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 mn l with low l.
The total cross section computed from pole figures has been used to analyse transmission experiment data (Santisteban et research papers 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.

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., 2004Willendrup et al., , 2019 (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 0 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 nfp ðÞ = expðÀÞ, where is the linear attenuation coefficient (or macroscopic cross section), given by 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 0 , where p 2 [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), k 0 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 k 0 . 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 k 0 .
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) 2 , 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 ofk k andq 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 Here, in general ^n n and '^n n denote the polar coordinates of the unit vectorn n in the sample reference system and P P where P n l ðxÞ are the associated Legendre functions defined in Appendix A. It is convenient to work with the polar components of V n l , defined by where R n l ðĜ GÞ and n l ðĜ GÞ are the modulus and the argument of the complex number V n l ðĜ GÞ, respectively. Then, equations (15) and (18) can be expressed as where 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½. . . replaced by sin½. . ., has to vanish. The computation of Â and Ç is too expensive if l max is large, even using expressions (23) and (24), with precomputed R m l ðGÞ and m l ðGÞ. As mentioned above, to speed up the computations U l ðG;n nÞ is precomputed, for each G and l, on a dense two-dimensional grid of cos ^n n 2 ½À1; 1 and '^n n 2 ½0; 2. 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, ÂðĜ G;q qÞ is precomputed, for each G, on a dense two-dimensional grid of cos ^q q and '^q 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 c Þ coh el , 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, k 0 . 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 where is the normalization factor. The values of P S ð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 P S ð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 Thus, as is well known,k k 0 lies on the surface of a cone whose axis is given byk k and whose angle is s (the Debye-Scherrer cone). It only remains to sample the azimuthal angle, ' 0 2 ½0; 2, around the cone axis. Introducing two unit vectorŝ t t 1 andt t 2 so that ft t 1 ;t t 2 ;k kg forms a right-handed orthonormal triad, we havê and Let us denote by P G ðk; ' 0 Þ the probability density of ' 0 . 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 forq q into ÂðĜ G;q qÞ, which is computed numerically from equation (23). Hence, given G and k, the first step is to obtain cos s , sin s , and the two vectorst t 1 andt t 2 . Then, ' 0 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 P G ðk; ' 0 Þ 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 ' 0 is sampled. The simple rejection method is thus straightforward and works as follows. An angle ' 0 is uniformly selected in [0, 2], andq q is computed according to equation (30). Then P G ðk; ' 0 Þ 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 ' 0 is accepted,k k 0 is computed from equation (29) and k 0 = kk k 0 . Otherwise, another value of ' 0 is selected uniformly in [0, 2] and the process is repeated.

research papers
Clearly, the sampling of ' 0 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 ÂðĜ G;q qÞ with reasonably efficiency (by interpolating precomputed values), this is not enough. For an importance sampling method, like rejection, we need the maximum in ' 0 of ÂðĜ G;q qÞ, withq q given by equation (30), which depends on the wavevector of the incident neutron, k. To estimate the maximum, ÂðĜ G;q qÞ has to be evaluated many times. Alternatively, we might consider using the weight factor transformation, 1 frequently used in McStas (Willendrup et al., 2019(Willendrup et al., , 2018. A similar problem arises in this case, however, since the weight is given by the normalized probability density, P G ðk; ' 0 Þ, which is proportional to ÂðĜ G;q qÞ, the proportionality factor being the inverse of the integral of ÂðĜ G;q qÞ over ' 0 . To compute the normalization factor, which again depends on k, one has to evaluate ÂðĜ G;q qÞ many times. A simple brute-force possibility would be to precompute either the maximum in ' 0 of ÂðĜ G;q qÞ or the normalization factor of P G ðk; ' 0 Þ 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 P G ðk; ' 0 Þ on the ' 0 grid each time that ' 0 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(Willendrup et al., , 2018. Using the SPLIT n keyword, each incoming neutron is reused n times. The values of coh and P G ðk; ' 0 Þ 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, fâ a;b b;ĉ cg, 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 ^n n ; '^n n Þ space where ÂðĜ G;q qÞ and U l ðĜ G;k kÞ are precomputed, n ct Â n ' ; and the size of the grid used to sample ' 0 , n ' 0 .
(c) One real number, the packing factor p.
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.

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 = P mn jC mn l j=ð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 ! 1.
The dependence on l max of the pole figures P G ðr 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) wherer 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 The dependence of the cross sections on l max is also interesting. Fig. 4  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 = P m;n jC l mn j=ð2l þ 1Þ. 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 ! 1 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. Fig. 5 displays the probability density, P G ðk; ' 0 Þ, of ' 0 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.  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.
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. 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.

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 ðĜ G;n nÞ and ÂðĜ G;q qÞ on a 201 Â 601 uniform grid in ðcos ; 'Þ 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ŷ y L direction in the laboratory reference frame, given by the orthonormal triad fx x L ;ŷ y L ;ẑ z L g. The sample is initially positioned so thatx x =x x L ,ŷ y = Àẑ z L andẑ z = y y L . The vectors of the sample reference frame, fx x;ŷ y;ẑ zg, 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 aboutx x and subsequently by an angle ' aboutẑ z, and a Monte Carlo simulation is performed. This process is repeated in steps of 5 in both and ', starting from = 0 and ' = 0. 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 ! 1 limit, they are normalized to the value at the largest l max .

Figure 5
The probability density, P G ðk; ' 0 Þ, of ' 0 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.

Figure 6
The same as Fig. 5 but for the Zr-2.5Nb pressure tube.     The same as Fig. 8 but for the Zr-2.5Nb pressure tube with l max = 30.
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.
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 k 0 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 À4 , 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 0 () and I 1 () 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 1 =I 0 Þ=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.
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 2 Ca 3 Al 2 F 14 (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.   (31) with the corresponding cut-off.

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 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 AE 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  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.

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). 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.
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.
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. The texture has an important influence on the SNR, as expected. The left-hand panel of Fig. 15 displays the ratio SNR 1 /SNR 2 , where SNR 1 and SNR 2 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 pointto-point ratio of Figs. 14(b) and 14(a).

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 expres-sion 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  (a) The ratio of SNR 1 /SNR 2 of the SNR of a pressure cell with Zr walls with the texture of the Zr-2.5Nb pressure tubes (SNR 1 ) and with a uniform texture (SNR 2 ). (b) A cut of panel (a) along the equatorial plane of the detector system, height = 0. 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 The irreducible representations of SO(3) are characterized by the non-negative integers l associated with the eigenvalues, l(l + 1), of L 2 = L 2 x þ L 2 y þ L 2 z , and have dimension 2l + 1. An orthonormal basis of the representation space is given by the eigenvectors jlmi of L 2 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 where L AE ¼ L x AE iL y . For each g 2 SO(3), the Wigner D matrix is defined by 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 Þ ¼ P l m 0 ¼Àl D l mm 0 ðg 1 Þ D l m 0 n ðg 2 Þ: In addition, the unitarity of the representation implies that the Wigner matrices are unitary and Below we will use the following relation (Galindo & Pascual, 1990): where Y m l ð; 'Þ are the spherical harmonics, defined by with 0 and 0 ' < 2, and P m l ðxÞ are the Legendre associated functions, given by P m l ðxÞ ¼ ð1 À x 2 Þ m=2 ðÀ1Þ m l!2 l d lþm dx lþm ðx 2 À 1Þ l : 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,n n, given by the polar angle and the azimuthal angle ', so that we can use the notation Y m l ðn 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: where S ðx x;ŷ yÞ is the Dirac delta function on the sphere, which in polar coordinates and ' reads S ðx x;ŷ yÞ ¼ ðcos ^y y À cos ^x x Þ P ð'^y y À '^x x Þ; ð45Þ and P (') is the periodic Dirac delta function. Then, introducing the Fourier representation of the ODF, f(g), we have research papers where I l mn ¼ R SOð3Þ dg D l mn ðgÞ S ðq q; gĜ GÞ; withq q ¼ ðk À k 0 Þ=jk À k 0 j. Let g^n n denote a rotation that brings a unit vectorn n toẑ z, so that g^n nn n ¼ẑ z. Since the Dirac delta function on the sphere is rotationally invariant, we have S ðq q; gĜ GÞ ¼ S ðẑ z; g^q q gg À1 G Gẑ zÞ: Using the invariance of the measure we have The integral can be readily evaluated in terms of the Euler angles, (, , ), that characterize g, since gẑ z ¼ sin cos x x þ sin sin ŷ y þ cos ẑ z; so that S ðẑ z; gẑ zÞ ¼ ðcos À 1Þ P ðÞ: Taking into account that d l mn ð1Þ ¼ mn , we obtain Z SOð3Þ dg D l m 0 n 0 ðgÞ S ðẑ z; gẑ zÞ ¼ 1 4 m 0 0 n 0 0 ; ð53Þ and therefore The Euler angles corresponding to a rotation g^n n that bringsn n toẑ z are = 0, ¼ ^n n and ¼ À '^n n , where ^n n and '^n n are the polar coordinates ofn n, and therefore, using relation (41), we have D l 0m ðg^n n Þ ¼ Using the above expression, equation (54) becomes 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 d k 0 the differential cross section given by equation (14). However, it is easier to integrate equation (12) over d k 0 before performing the integral over dg. Since the scattering is elastic, we have k 0 ¼ kk k 0 , so that ðk À k 0 À gGÞ ¼ ðk À gG À kk k 0 Þ ¼ 1 k 2 ðjk À gGj À kÞ Sk k À 1 k gG;k k 0 : The integral over d^k k 0 of the above expression is 1, since the first Dirac delta function forcesk k À gG=k to lie on the unit sphere. Therefore we have Let us remember that we defined g^n n as the rotation that brings the unit vectorn n toẑ z, so that g^n nn n ¼ẑ z. Then we have Using the invariance of the measure we can write R dg f ðgÞ kðjk À gGjÀ kÞ ¼ R dg f ðg À1 k gg G Þ ðjkẑ z À Ggẑ zjÀ kÞ; and, using the Fourier expansion of the ODF and the group properties of the Wigner matrices, we have R dg f ðgÞ kðjk À gGj À kÞ ¼ P 1 where S l m 0 n 0 ðG; kÞ ¼ R dg D l m 0 n 0 ðgÞ kðjkẑ z À Ggẑ zj À kÞ: Using equation (51) we have jkẑ z À Ggẑ zj ¼ k 2 þ G 2 À 2kG cos À Á 1=2 ; which is independent of and . Therefore, the integrals over these angles give 4 2 m 0 0 n 0 0 . Making the change of variable z = cos and taking into account that d l 00 ðzÞ = P l ðzÞ, we have S l m 0 n 0 ðG; kÞ ¼ m 0 0 n 0 0 1 2 Z 1 À1 dz P l ðzÞ k k 2 þ G 2 À 2kG cos À Á 1=2 Àk h i : The argument of the delta function vanishes for z = G/2k, and its derivative at this point is ÀG. Then we have S l m 0 n 0 ðG; kÞ ¼ m 0 0 n 0 0 k 2G P l G 2k Inserting the above result into (61) and using equation (55) we get Z dg f ðgÞ k ðjk À gGj À kÞ ¼ k 2G ÇðG; kÞ: