Multiple Bragg reflection by a thick mosaic crystal. II. Simplified transport equation solved on a grid.

The generalized Darwin-Hamilton equations [Wuttke (2014). Acta Cryst. A70, 429-440] describe multiple Bragg reflection from a thick, ideally imperfect crystal. These equations are simplified by making full use of energy conservation, and it is demonstrated that the conventional two-ray Darwin-Hamilton equations are obtained as a first-order approximation. Then an efficient numeric solution method is presented, based on a transfer matrix for discretized directional distribution functions and on spectral collocation in the depth coordinate. Example solutions illustrate the orientational spread of multiply reflected rays and the distortion of rocking curves, especially if the detector only covers a finite solid angle. open access.

Introduction

In a preceding paper, designated as Part I (Wuttke, 2014a ▸), multiple Bragg reflection from a thick, ideally imperfect crystal was studied mainly by analytical means. The planar two-ray transport equations of Darwin (1922 ▸) and Hamilton (1957 ▸) were generalized to account for out-of-plane trajectories. Expanding these equations into a recursive scheme led to some asymptotic results, but did not provide a practicable solution algorithm for the generic case with crystals of finite thickness. Reflection probabilities depend strongly on propagation directions, and with each reflection the next reflection probability can vary by orders of magnitude. This makes the transport equations ill conditioned, and straightforward Monte Carlo simulations inefficient and unreliable.

In this paper, a completely different solution method is presented. Instead of following individual rays through forward and backward reflections, we study reflection-order-independent fluxes (current distributions) I as a function of propagation direction and penetration depth z. They are governed by a system of linear ordinary differential equations in z with separated boundary conditions [equations (1) and (4) below]. We present spectral collocation as a practicable solution method. Solutions are iterated with increasing numbers of collocation points until a required accuracy is reached. Our algorithm is fast enough to be used interactively or/and within complex instrument simulations.

This paper can be read independently of Part I. We re-de­ri­ve most of the theory, making it simpler and more generic. By consequential use of energy conservation, we get rid of one phase-space dimension. By positing translational invariance along the surface of the mosaic plate, two other dimensions are eliminated. This reduction to three argument dimensions is the precondition for an efficient numeric solution on a grid.

Large parts of the theory are now formulated coordinate free. Block normals that fulfill the Bragg condition are parameterized by a polar instead of a Cartesian coordinate; this eliminates an apparent singularity that forced us in Part I to exclude near-backscattering from consideration. Furthermore, the transport equations are simplified and generalized by removal of any reference to two distinct beams.

While all worked-out examples assume a simple geometry with the mean crystallite normal collinear to the mosaic normal, our formalism can be used with any other orientational distribution. One application we have in mind is beam deflection by a rotating stack of tilted mosaic crystals of highly oriented pyrolytic graphite as used in the phase-space transform chopper of third-generation neutron backscattering spectrometers (Meyer et al., 2003 ▸; Frick et al., 2006 ▸; Wuttke et al., 2012 ▸).

In Section 2, we derive the mathematical model to be studied. Discrete grids are chosen in Section 3. In Section 4 our numeric solution method is presented and verified against the two-ray model. Example solutions are shown in Section 5 and conclusions drawn in Section 6. Some derivations, computational details and special cases can be found in Appendices A –D . The supporting information provides the source code and additional documentation of the software MultiBragg developed along with this work.

The mathematical model

Crystal model and current distribution

Following Darwin (1922 ▸), a mosaic crystal is modeled as an assembly of perfectly crystalline blocks that are to some degree orientationally disordered. In an ideally imperfect crystal every block is so thin that it reflects at most a small fraction of the incident beam. Primary extinction and multiple reflections within a block can be neglected. As in Part I, we consider a thick, ideally imperfect crystal, consisting of so many block layers that secondary extinction and multiple reflections are of practical importance.

Since reflections from different blocks add incoherently, the adequate description level is classical transport theory. Our task is to compute the stationary flux (current distribution) . Only elastic diffraction shall be fully accounted for. Inelastic scattering will be dealt with by a loss term. Accordingly, the wavenumber k is a conserved quantity, and can therefore be dropped from the argument list of , leaving over a dependence on the propagation direction .

As in Part I, we concentrate on a mosaic crystal in the form of a slab that can be approximated as an infinite plate (Fig. 1 ▸). Altogether the flux is projected from six-dimensional phase space to the three-dimensional function . This opens the possibility of solving the boundary problem with manageable effort on a grid, thereby overcoming the limitations of the Monte Carlo method used in Part I.

Figure 1

Bragg reflection by a crystalline block within a mosaic plate. Block normals are distributed around , as indicated by the blue cone. The angular width of the distribution is grossly exaggerated; typically, it is a few degrees only. Here, and in all specific examples in this work, we have chosen to be collinear with the real-space depth direction . Much of our theory also holds for .

The price is that we have to neglect the lateral displacement of the beam, which is correlated with the reflection order, which is correlated with the directional spread. At least in the above-mentioned application scenario (graphite deflector in a neutron spectrometer, far from grazing incidence), this is harmless: the mean lateral displacement is at most a low multiple of the crystal thickness, which is a few millimetres, and therefore corresponds to a fraction of a degree at the next optical element, located 2 m downstream, whereas the deflector crystals have a rocking width of several degrees.

Transport equation and boundary conditions

The flux obeys the transport equation (Sears, 1989 ▸, equation 8.1.24),a stationary Boltzmann equation with drift and scattering terms. The linear operator B describes gains by Bragg diffraction, where is the solid-angle differential associated with the integration variable . The kernel μ is reviewed below in Section 2.3. The attenuation operator A is a multiplicative factor, The integral accounts for losses by diffraction. The constant stands for absorption, inelastic scattering, diffuse scattering and diffraction by parasitic reflections (Dorner & Kollmar, 1974 ▸).

To specify boundary conditions, we consider an infinite plate of thickness d, extending from to . The incident flux comes from the half space z < 0. Accordingly, the boundary conditions areOur task is to compute the reflected and the transmitted flux with the indicator function [true] = 1, [false] = 0.

In Part I, we had divided into two functions, , representing the forward and backward traveling beam. Accordingly, the transport equation consisted of two coupled differential equations, generalizing the two-ray Darwin–Hamilton equations of the conventional planar approximation. That notation was useful for describing the reflection-order expansion [Section 3.2 of Part I; see also Grabcev & Stoica (1980 ▸)] as a zigzag walk (Wuttke, 2014b ▸) with a strictly alternating sign of . To distinguish two beams we had to exclude the case of grazing incidence. The notation (1), with just one function I defined for all , is simpler, more generic and more convenient for our present purpose. Only later, when we choose a grid in to approximate I by a histogram, will we take into consideration the effective two-beam geometry.

Reflection kernel

The transport kernel in (2) is a transfer function that gives the probability per unit length for a particle with incident direction to be scattered into an infinitesimal solid angle around the outgoing direction . It is a sum over single-reflection transfer functions, given by an integral over the block transfer function (47) derived in Section A1. The scattering directions depend on the block orientations, and have the statistical distribution .

In certain situations, multiple diffraction by multiple Bragg reflections can be of practical importance (Ohmasa et al., 2016 ▸). Nonetheless, to simplify our exposition, we shall consider only one pair of reflections, hkl and . With the joint distribution we can merge (6) and (7) into the total transfer function An integration, explained in Section A2, reduces the total transfer function (9) toThis simplifies in several ways equation I,25 [denoting equation number (25) in Part I], as discussed in Section A3. We now define the symbols , , and introduced with (10).

The prefactor depends on the unit-cell volume V and structure factor . It has the dimension of an absorption coefficient, i.e. inverse length. The Bragg angle is constant because we consider a fixed reflection hkl and a constant radiation wavenumber k. Both and are independent of the sign of the reflection. The outgoing beam direction is given by the deflection function The parametric curve with contains all possible scattering directions that satisfy the Laue–Bragg condition for an incoming wave direction .

To construct , we choose an orthonormal base for the reciprocal-space vectors and . Note that is not required to coincide with the plate normal (though it does so in our code and our worked-out examples). Choose a rotation matrix so that (for readability, we omit carets in subscripts). The circle of possible scattering directions can then be written It is straightforward to verify that , for all t, satisfies (13).

The condition leaves underdetermined, allowing for an arbitrary rotation around . This is irrelevant because the origin of the polar coordinate t is arbitrary, and only appears under integrals that run from to .

Specializing the distribution of scattering directions

For most mosaic crystals, W is isotropic, i.e. invariant under rotation around the mean block normal . Thereby depends only on .

All the following theoretical developments, including the numeric methodology of Section 4, are independent of what isotropic distribution we choose for . For our numeric examples, however, we need to be more specific. In certain cases (Ohmasa & Chiba, 2018 ▸), W can be a ring-like distribution. Here we concentrate on 00l reflections, where scattering vectors are parallel to the block normals so that is a disc-like distribution. As is customary, we will choose a Mises–Fisher (MF) distribution (a Gaussian on the unit sphere), with the normalization constant . Usually, there is negligible overlap between and so that the sum (8) can be approximated as A mosaic with shall be called normal oriented. Some consequences of the rotational symmetry around are discussed in Appendix C .

In all numeric examples we assume an isotropic, normal oriented mosaic, and we choose . Unless differently stated, the standard deviation is η = 2.5°. The orthographic projection of the circle into the plane is an ellipse. Fig. 2 ▸ shows examples and puts them in relation to .

Figure 2

Crystal orientations that fulfill the Laue–Bragg condition, projected into the plane, form ellipses. The two plots have different Bragg angles . Each plot shows ellipses for three different incident angles , with . The concentric gray discs contain 50% and 90% of all mosaic blocks, assuming a Mises–Fisher distribution that is centered around , with standard deviation η = 2.5°.

Parameterization

In our numeric examples, we will characterize crystals by two dimensionless constants that have a simple intuitive meaning in the two-ray limit (68) for a collimated beam with incident angle . The first of these parameters is the opacity where . The second dimensionless crystal parameter is the relative reflectivity These parameters will enter the following derivations through the products and Note that are not pure material constants but also depend on the wavelength of the used radiation. So, in principle, one could tune to almost arbitrary values by suitable combinations of wavelength and crystal thickness.

Discretization in

Binning

So far, we have assumed that as a function of is a distribution on the unit sphere. We now request that the relevant regions of the unit sphere be partitioned in M bins, and we replace by M histograms with . Each histogram represents a current that is defined as flux integral over the solid angle , Combining this with the definition (2) of the operator B, we getWe assume that μ is a sufficiently smooth function of so that it can be drawn in front of the second integral. We obtain with The attenuation factor, discretized in analogy with (24), is and the transport equation (1) takes the form In this paper, we will not investigate errors introduced by the approximation (23). It is up to practitioners to choose appropriate histogram grids so that both discretization errors and computing time be kept within reasonable bounds.

Grids in zero, one, two dimensions

In Section 4, spectral collocation will be introduced without reference to a particular histogram grid. For our numeric examples, we choose three different grids.

The smallest meaningful grid consists just of bins, representing a forward and a backward traveling beam, with and , respectively. It will be used in Figs. 5 and 6 to illustrate our approach in the simplest possible way, and to allow verification against the known analytical solution. Instead of the indices r = 1, 2, we will use the signs ± to denote the beam direction.

If we are only interested in the total intensity or in the polar distribution of radiation reflected or transmitted by an isotropic, normal oriented mosaic, as in Figs. 9, 10, 12, then we can take an azimuthal average (Section C3), and solve the transport problem on a fine-grained one-dimensional grid in θ.

In all other cases, we need a two-dimensional partition of the unit sphere. Any possible map projection and coordinate system can be chosen to construct the grid. In our examples, we want to preserve the symmetry of the isotropic, normal oriented mosaic and therefore choose a rectangular grid in the spherical coordinates θ and φ.

The one- or two-dimensional grids must not necessarily cover the entire unit sphere. For computational efficiency, we restrict them to two contiguous regions around the transmitted and the reflected beam. These regions can be iteratively adapted, keeping the cut-off error (estimated from the loss channel b, Section B2) under a given tolerance (Section 4.6).

Diffraction matrix

The integral (24) can be carried out at once since the kernel , given by (10), contains a delta function. The result is with the indicator bracket as introduced in Section 2.2.

For given , and sweeping t, the outgoing directions form a one-dimensional manifold on the two-dimensional sphere. This is illustrated in Fig. 3 ▸, which shows these manifolds for three different s. In consequence, for a two-dimensional histogram grid, the matrix B is sparse: most entries are zero, the more so the finer the grid. If each of the two coordinate axes is divided into bins, then B has nonzero entries.

Figure 3

The three bands represent three columns of the reduced reflectivity matrix , with and with three different values of , shown as a function of and . The Bragg angle is = 70°. There are 80 φ bins from 55° to 85°, and 180 θ bins from −180° to 180°. The dimensionless intensity scale applies for and .

Section B1 presents an algorithm for the actual computation of (27). In Section B2, the M directional bins are extended by three loss bins. One of them accounts for absorption; the other two allow us to quantify unphysical losses originating from numeric cut-offs. Therefore we can detect violations of particle conservation, and quantify, and ultimately control, numeric approximation errors.

Spectral collocation in the depth coordinate

Depth rescaling

As we will use Chebyshev polynomials in the depth coordinate, it shall be transformed from to the standard range . We therefore introduce the reduced coordinate and the transformed histogramsThe transport equation (26) becomes with the separated boundary conditions Grids should be constructed such that no bin crosses the equator of the unit sphere.

Equation system

The equation system (30) consists of M coupled first-order linear differential equations in . While a formal solution can easily be written down as a matrix exponential, it is numerically not viable (Moler & Loan, 1978 ▸, 2003 ▸). The method of choice for this kind of problem is spectral collocation; it is based on the expectation that the solution is a smooth function of ζ and therefore can be expanded in Chebyshev polynomials (Gottlieb et al., 1984 ▸; Canuto et al., 1988 ▸; Trefethen, 2000 ▸).

We approximate the functions by polynomials of order N that match in collocation points . We defer the choice of to Section 4.3; until then, we only require . Function values at the collocation points are abbreviated These are unknowns. M of them can immediately be read off from the boundary conditions (31): The others will be obtained from the transport equation (30). To discretize this differential equation, we replace by , by a differentiation matrix D, specified by the requirement[for an introduction to differentiation matrices, see Trefethen (2000 ▸)]. The resulting MN equationsmust hold for all histogram bins () and in all collocation points ().

We collect all linear operators in the matrix with the Kronecker delta so that the transport equation (35) becomes simply These are equations in variables , of which M are known from (33).

The matrix L is sparse because of the Kronecker deltas in its definition (36) and because its component is also sparse when it matters, namely in the computing-intensive case of a two-dimensional grid. In that case, per Section 3.3, of the entries of matrix B, only are nonzero. Overall, of the entries of L, only are nonzero. This is visualized in Fig. 4 ▸.

Figure 4

Visualization of a small part of a small matrix L. Parameters: = 70°, , . Discretization: three collocation points in z; three bins for 55° ≤ θ ≤ 85°; 36 bins for −180° ≤ φ ≤ 180°. Only the entries with 0 ≤ φ ≤ 60° are shown. Since some entries are negative, the figure shows absolute values .

In view of the boundary conditions (4), we now distinguish between histogram bins with forward and backward propagation direction, according to the sign of . We split the sum over s in (37) accordingly, omitting the zero terms in with backward s, and bringing the nonzero terms in with forward s to the right side: This system of inhomogeneous linear equations in MN unknown is overdetermined, due to the loss of information that goes along with the reduction of polynomial order in differentiation.

Overdetermination can in principle be avoided by using a rectangular differentiation matrix (Driscoll & Hale, 2016 ▸; Xu & Hale, 2016 ▸). However, this would be unsuitable for the full multi-ray problem because the necessary ‘downcasting’ interpolation of linear terms would make the matrix L much less sparse. We rather opt for the standard procedure of simply ignoring redundant equations. We choose to delete the M equations with , r forward or , r backward.

Collocation points and differentiation matrix

All of Section 4.2 holds regardless of the collocation points. Their choice, however, is of critical importance for the resulting convergence. We make the standard choice of Chebyshev–Gauss–Lobatto points, which are located at the extrema of the Chebyshev polynomial , They only enter our computation through the corresponding differentiation matrix (34).

This matrix has the outer diagonal entries the interior diagonal entries and the diagonal endpoints For a derivation, see e.g. Gottlieb et al. (1984 ▸) or Trefethen (2000 ▸), but note that our choice of ascending has led to a minus sign on the right-hand side of (42).

Collocation error for the two-ray reference

By the collocation error we understand the error caused by approximating the functions by polynomials . We first consider the two-ray approximation (re-derived in Section D1) for which we can determine the collocation error by comparing with the known analytical results (summarized in Section D2). We choose so that the Bragg operator (68) is simply and the total attenuation .

Fig. 5 ▸ shows currents versus ζ for a moderately thick crystal with realistic attenuation: , . On the linear scale of this plot, the numeric data points and the analytical curves (70) agree perfectly for collocation orders as low as .

Figure 5

Directional currents (29) as a function of depth for , . Lines show the analytical solution (70). Symbols have been computed by spectral collocation with different N.

The rapid convergence of the spectral collocation is further demonstrated in Fig. 6 ▸, which shows the deviation from the analytical result as a function of N. The decrease of the error with increasing N is roughly exponential until some base level is reached.

Figure 6

Accuracy of spectral collocation as a function of the number of collocation points N, for the planar two-ray model, for different crystal opacities. (a), (b) Absolute error of the transmitted and reflected current, determined by comparison with the analytical result (71). (c) Deviation of the total current (transmitted, reflected and absorbed) from the true value 1.

Note that the figure represents the absolute error. This could become a problem in shielding calculations for very thick mosaic crystals where a tiny transmittivity would result in an unacceptable relative error. We exclude this peculiar case from further consideration.

Fig. 6 ▸(c) shows the error of the total current where is the absorption loss channel (Section B2). Analytically, the true value is 1. For odd N, convergence takes about as long as in Figs. 6 ▸(a), 6 ▸(b). In contrast, even for the smallest even N the error is only of the order of machine precision, thanks to pairwise cancelation at collocation points i, . Therefore we generally prefer even N in our computations.

Collocation error for the full model

We now investigate the collocation error for a two-dimensional grid, for representative crystal parameters and for an ideal collimated incoming beam We only consider the reflected flux .

In contrast to Section 4.4, no analytical solution is known. Therefore we estimate the overall collocation error by comparing approximate solutions at successive collocation orders N, In Figs. 7 ▸ and 8 ▸, we increase N from in steps of ; elsewhere, larger values of the hyperparameters and are more convenient (Section 4.6).

Figure 7

Estimates of the collocation error of the reflected current as a function of the collocation order N, for different solver convergence tolerances. Model parameters , , = 70°. Collimated incoming beam (43) with = 70°. Discretization: 45 bins for 58° ≤ θ ≤ 82°; 11 bins for −180° ≤ φ ≤ 180°.

Figure 8

Minimal collocation order N required to keep the overall error (44) of the reflected flux below . In each graph, one of the parameters of the default model of Fig. 7 ▸ is varied. In each graph, a red disc indicates the parameter value that is used in all other graphs.

In Fig. 7 ▸, the error estimate is shown as a function of N for different solver tolerances. The behavior is known from Fig. 6 ▸(b): after a roughly exponential decrease the cross over to a noisy base level that lies safely below the solver tolerance.

The collocation order N is determined by repeating the entire computation for increasing values of N until lies below a given bound (Section 4.6). Fig. 8 ▸ shows how the so-determined minimal collocation order N varies with the grid size and with the model parameters , , η, ρ, τ. All these parameters are found to be uncritical, except the opacity τ, for which Fig. 8 ▸ suggests a possible asymptote . Therefore our computational method should not be applied to extremely thick crystals with . This is of no concern for reflectivity computations: it is always possible to restrict τ to values of about 20 or 30; anything beyond is inconsequential for the reflected current.

Hyperparameters and a posteriori tolerances

The numeric solution is controlled by hyperparameters (so-called in opposition to the regular parameters that describe the crystal model):

(i) Bounds in θ and φ, and numbers of bins, to specify the directional grid.

(ii) Initial number of collocation points and increment for the iterative procedure described in the last paragraph of Section 4.5. Our default choice is and .

(iii) Auxiliary parameters for the numeric evaluation of the reflection kernel (27), as explained in Section B1: bisection accuracy , maximum step size and a kernel cut-off .

(iv) The collocation tolerance (Section 4.5): the sparse equation solver is called with increasing N until the estimated error [equation (44)] falls below . We use .

(v) The solver tolerance required by the sparse equation solver (Section 4.7). It should be smaller than . We use a value of , except when we determine collocation errors as a function of the number of collocation points: Fig. 6 ▸ is computed with and Fig. 7 ▸ with .

Additionally, some a posteriori tolerances are used in checking the numeric integrity of the obtained solution. If any check fails, then we recompute everything with stricter hyperparameters. These tolerances are:

(i) A noise level to disallow currents below , as may be caused by numeric inaccuracies at weak intensities.

(ii) A stricter limit is imposed to the absolute value of the total current (53).

(iii) The cut-off tolerances and are upper limits for the loss channels and that account for unphysical losses due to the finite grid and the deletion of weak matrix elements (Section B2).

All our numeric examples have been checked with and .

Numeric tools

A sparse LU solver is used to invert the remaining MN equations. Our implementation is based on tools from the numerical library Trilinos (Heroux et al., 2003 ▸), namely the sparse compressed row matrix class from package Epetra, the ILU(0) preconditioner from package Ifpack (Sala & Heroux, 2005 ▸) and the GMRES (generalized minimal residual) block solver from package Belos (Bavier et al., 2012 ▸). The mapping of array indices is described in the supporting information.

Results

Open-source code MultiBragg

The software developed along with this work is released under the GNU Public License (GPL v3 or higher), and deposited in the form of a compressed tar.gz archive as supporting information. The code comprises a library MultiBragg for the numeric solution of the transport equation, and application programs that generate the data for the figures in this work. More information on the software is provided in the textual part of the supporting information.

Total reflectivity

Fig. 9 of Part I showed the total reflectivity R for plates of different opacities as a function of the Bragg angle , computed by Monte Carlo integration. The present Fig. 9 ▸ compares old and new results. The perfect accord of both data sets provides strong support for the correctness of both computer codes. While error bars in Part I were of the order , our new results are far more accurate and extend over a wider range. Our new method is also much faster: computing times are typically of the order of some minutes, and only of a few seconds if no azimuthal resolution is needed, as here and in Figs. 10 and 12.

Figure 9

Total reflectivity as in case 1 of Fig. 9 of Part I. Collimated incoming beam with ; mosaicity = 0.025 rad; nominal opacity ; relative reflectivity . In terms of Part I: . Blue symbols with error bars are from the Monte Carlo computation of Part I; red circles are from the present numeric integration.

In the range covered by Part I, corrections from non-planarity amounted to 1% at most. Our numeric method now allows us to compute R up to = 90°. Representative results are shown in Fig. 10 ▸. With approaching 90°, the reflectivity first increases, then decreases rapidly towards 0. These observations are easily understood by looking back to Fig. 2 ▸: close to backscattering, the ellipse of block orientations that fulfill the Bragg condition is almost a circle, and for it is concentric with the disc representing . This makes it plausible that the reflectivity can rise more than 15% above the constant value from planar theory. Even closer to backscattering, however, the ellipse shrinks towards a point, fewer blocks are available for Bragg diffraction, and the reflectivity decreases proportionally.

Figure 10

Total reflectivity near backscattering, for a collimated incoming beam with ; mosaicity η = 2.5°.

Azimuthal distribution

Fig. 11 ▸ shows the directional distribution of the transmitted and reflected radiation for three different incoming beam inclinations . Since the coordinate representation does not account for the different bin sizes , this figure does not show currents (flux integrals per bin) but the directional flux .

Figure 11

Directional distribution of the transmitted and reflected flux for three different inclinations of the incoming collimated beam. Crystal parameters: = 70°, η = 2.5°, , .

In transmission, a bright spot shows the attenuated incoming beam. It is least pronounced at , where the reflectivity is strongest. In the reflected distribution, the parabolic trace comes from one-reflection trajectories, whereas the diffuse cloud represents the sum of all higher reflection orders. While the parabola is known from the approximative treatment of Hennig et al. (2011 ▸), the two-dimensional cloud is only accounted for by the full transport equation solved here.

To visualize the relative importance of this cloud, Fig. 12 ▸ shows the intensity of the direct beam and the single-reflected spray relative to the total transmitted or reflected intensity. Multiple reflections are most important for high opacity τ, high relative reflectivity ρ, and for incident beam directions close to the Bragg condition, . For and , the relative importance of direct transmission goes quickly to 0 for increasing τ. The relative importance of single reflections decreases more slowly, and the limit of 50% is not fully attained within the τ range of the figure.

Figure 12

(a), (b) Relative contribution of the direct beam to the total transmitted intensity; (c), (d) relative contribution of single reflections to the total reflected intensity. All data for = 70°; (a), (c) as a function of τ for different values of ρ, with ; (b), (d) as a function of for different combinations of τ and ρ.

Rocking curves

The standard way to characterize a mosaic crystal experimentally is by measuring a rocking curve (e.g. Schneider, 1974 ▸): the reflected or transmitted intensity is recorded while the crystal is rotated around an axis normal to the scattering plane. In our fixed-crystal frame, this is equivalent to scanning the incident angle while maintaining the detector angle at .

It is well known from theory and experiment (Dorner & Kollmar, 1974 ▸) that rocking curves are generally wider than the underlying crystallite orientation distribution . The width increases with increasing opacity τ. This is illustrated by Fig. 13 ▸ where rocking curves for a very thin () and a very thick () mosaic are shown. In the thin-crystal limit, our solution of the full Darwin–Hamilton equations reproduces Sears’ solution of the planar approximation (71), and both curves coincide almost perfectly with the Gaussian .

Figure 13

Rocking curves for (a) a thin, (b) a thick mosaic crystal with and , respectively. The gray area indicates the Gaussian mosaic distribution with standard variation η = 2.5° as assumed throughout this work. The dashed line is the reflectivity in planar approximation (71). Symbols represent numeric solutions of the full transport equation: small colored symbols show the intensity collected by circular detectors with radius specified as an angle; thick black circles show the entire reflected radiation.

Conversely, for the thick crystal Sears’ solution is much wider than , and our full solution deviates from Sears’ in that it is shifted by about 1° and has a slight asymmetry. This confirms the Monte Carlo result of Fig. I,10, and extends it to Bragg angles further away from backscattering.

So far, we have discussed total reflected intensities. In practice, detectors cover only a finite solid angle. Rocking curves as measured by circular detectors are shown by the colored open symbols in Fig. 13 ▸. Unless the opening angle is considerably larger than η a considerable part of the reflected intensity is indeed lost outside the detector. For the thick crystal, the shape of the rocking curve also varies considerably with the angular coverage.

Conclusions

To summarize, we have simplified the transport equation of Part I (Wuttke, 2014a ▸) by making consequential use of energy conservation and projecting everything to the sphere .

For isotropic, normal oriented mosaics (Section 2.4), azimuthal current distributions are insensitive to the resolution in the polar angle φ; if the polar distribution does not matter, then numeric computations can be accelerated by considering one single φ bin (Sections 3.2, C3). From there, only one more linearization (Section D1) is needed to explain how the original planar Darwin–Hamilton equations got the integral currents essentially right.

Our first numeric result (Section 5.2, Fig. 9 ▸) confirms that off-plane trajectories have very little effect upon the integral currents except near backscattering. The interest of our present work is not in those minor corrections, but in deriving information that is not at all available from the original Darwin–Hamilton equations, namely the directional distribution of the transmitted and reflected radiation.

While Wuttke (2014a ▸) presented some asymptotic results, a formal expansion in reflection order and a Monte Carlo code, we now have derived a numeric scheme that uses spectral collocation in the depth coordinate to compute with high speed and very high accuracy. Fig. 11 ▸ shows an example outcome: a dot represents the direct beam, a parabolic spray comes from single reflections, whereas all higher reflection orders contribute to a diffuse, two-dimensional cloud of propagation directions. Fig. 13 ▸ shows the consequences for rocking-curve measurements.

As mosaic crystals are an important beam optical device, numeric solutions of the transport equation will help to improve instrument and radiation protection simulations. The computer code produced for this work is open source and freely available, and will hopefully find its way into established ray-tracing packages.

Click here for additional data file.

Open-source software MultiBragg as used to generate Figs 3-13. DOI: 10.1107/S2053273320002065/ae5082sup1.gz

Table 1

Correspondence of notations in Sears (1997 ▸) and in the present work

Variable	Sears (1997 ▸)	This work
Crystal thickness	d	d
Depth variable
Rescaled depth variable
Bragg angle	θ
Forward current	I
Backward current
Forward beam polar angle	φ
Backward beam polar angle
Bragg cross section	σ	B
Bragg cross section, reduced, forward beam
Bragg cross section, reduced, backward beam
Loss cross section	μ
Total attenuation cross section
Total attenuation cross section, reduced, forward beam
Total attenuation cross section, reduced, backward beam