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

Bornemann, F.; Li, Y.Y.; Wuttke, J.

doi:10.1107/S2053273320002065

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 76| Part 3| May 2020| Pages 376-389

https://doi.org/10.1107/S2053273320002065

Open

access

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

Folkmar Bornemann,^a Yun Yvonna Li ^b and Joachim Wuttke ^b ^*

^aZentrum Mathematik M3, Technische Universität München, 85747 Garching, Germany, and ^bJülich Centre for Neutron Science at MLZ, Forschungszentrum Jülich GmbH, 85747 Garching, Germany
^*Correspondence e-mail: j.wuttke@fz-juelich.de

Edited by A. Altomare, Institute of Crystallography - CNR, Bari, Italy (Received 28 November 2019; accepted 13 February 2020; online 16 April 2020)

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.

Keywords: mosaic crystals; multiple scattering; Darwin–Hamilton equations; spectral collocation.

1. 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 $[\hat {\bf k}]$ 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-derive 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 $[\hat {\bf G}]$ 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 $[\hat {\bf k}]$ 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.

2. The mathematical model

2.1. 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) $[I({\bf k},{\bf r})]$ . 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 $[I({\bf k},{\bf r})]$ , leaving over a dependence on the propagation direction $[\hat {\bf k}]$ .

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 $[I(\hat{\bf k},z)]$ . 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 $[\hat {\bf G}]$ are distributed around $[\hat{\bf H}]$ , 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 $[\hat{\bf H}]$ to be collinear with the real-space depth direction $[\hat{\bf z}]$ . Much of our theory also holds for $[\hat{\bf H}\not\parallel\hat{\bf z}]$ .

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.

2.2. Transport equation and boundary conditions

The flux obeys the transport equation (Sears, 1989 , equation 8.1.24),

$[(\hat {k}_z \partial_z + A)I(\hat {\bf k},z) - B\circ I(\hat {\bf k},z) = 0, \eqno(1)]$

a stationary Boltzmann equation with drift and scattering terms. The linear operator B describes gains by Bragg diffraction,

$[B\circ I(\hat {\bf k}): = \textstyle\int\,{\rm d}^2\hat {\bf k}' \mu(\hat {\bf k},\hat {\bf k}') I(\hat {\bf k}'),\eqno(2)]$

where $[{\rm d}^2\hat {\bf k}']$ is the solid-angle differential associated with the integration variable $[\hat {\bf k}']$ . The kernel μ is reviewed below in Section 2.3. The attenuation operator A is a multiplicative factor,

$[A(\hat {\bf k}): = \textstyle\int\,{\rm d}^2\hat {\bf k}'\mu(\hat {\bf k}',\hat {\bf k}) + \mu_{\rm a}.\eqno(3)]$

The integral accounts for losses by diffraction. The constant $[\mu_{\rm a}]$ 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 z = 0 to z = d. The incident flux I_in comes from the half space z < 0. Accordingly, the boundary conditions are

$[\matrix{ I(\hat {\bf k},0) = I_{\rm in}(\hat{\bf k}) &{\rm for}\,\,\hat {k}_z\,\gt\,0,\cr I(\hat {\bf k},d) = 0\hfill &{\rm for}\,\,\hat {k}_z\,\lt\,0.\cr }\eqno(4)]$

Our task is to compute the reflected and the transmitted flux

$[\matrix{ I_{\rm refl} (\hat{\bf k}) &: = [\hat {k}_z\,\lt\,0]I(\hat{\bf k},0),\cr I_{\rm trans}(\hat{\bf k}) &: = [\hat {k}_z\,\gt\,0]I(\hat{\bf k},d)\cr }\eqno(5)]$

with the indicator function [true] = 1, [false] = 0.

In Part I, we had divided $[I({\bf k},{\bf r})]$ into two functions, $[I_\pm]$ , 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 k_z. To distinguish two beams we had to exclude the case of grazing incidence. The notation (1), with just one function I defined for all $[\hat {\bf k}]$ , is simpler, more generic and more convenient for our present purpose. Only later, when we choose a grid in $[\hat {\bf k}]$ to approximate I by a histogram, will we take into consideration the effective two-beam geometry.

2.3. Reflection kernel

The transport kernel in (2) is a transfer function that gives the probability per unit length $[\mu(\hat {\bf k}',\hat {\bf k}){\rm d}^2\hat {\bf k}']$ for a particle with incident direction $[\hat {\bf k}]$ to be scattered into an infinitesimal solid angle $[{\rm d}^2\hat {\bf k}']$ around the outgoing direction $[\hat {\bf k}']$ . It is a sum

$[\mu(\hat {\bf k}',\hat {\bf k}) = \textstyle\sum\limits_{hkl}\mu_{hkl}(\hat {\bf k}',\hat {\bf k})\eqno(6)]$

over single-reflection transfer functions, given by an integral

$[\mu_{hkl}(\hat {\bf k}',\hat {\bf k}) = \textstyle\int\,{\rm d}^2\hat {\bf G}W_{hkl}(\hat {\bf G})\mu^{\rm block}_{hkl}(\hat {\bf k}',\hat {\bf k}\semi \hat {\bf G})\eqno(7)]$

over the block transfer function $[\mu^{\rm block}_{hkl}]$ (47) derived in Section A1. The scattering directions $[\hat {\bf G}]$ depend on the block orientations, and have the statistical distribution W_hkl.

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 $[{\overline h}\, {\overline k}\, {\overline l}]$ . With the joint distribution

$[W(\hat {\bf G}): = W_{hkl}(\hat {\bf G}) + W_{{\overline h}\, {\overline k}\, {\overline l}} (\hat {\bf G})\eqno(8)]$

we can merge (6) and (7) into the total transfer function

$[\mu(\hat {\bf k}',\hat {\bf k}) = \textstyle\int\,{\rm d}^2\hat {\bf G}W(\hat {\bf G})\mu_{\rm block}(\hat {\bf k}',\hat {\bf k}\semi \hat {\bf G}).\eqno(9)]$

An integration, explained in Section A2, reduces the total transfer function (9) to

$[\mu(\hat {\bf k}',\hat {\bf k}) = \overline\mu \cos\theta_{\rm B} \textstyle\int\limits_{0}^{2\pi}\,{\rm d}tW[\hat {\boldGamma}(\hat {\bf k},t)]\delta^2[\hat {\bf k}'-\hat {\boldkappa}(\hat {\bf k},t)].\eqno(10)]$

This simplifies in several ways equation I,25 [denoting equation number (25) in Part I], as discussed in Section A3. We now define the symbols $[\overline\mu]$ , $[\theta_{\rm B}]$ , $[\hat {\boldGamma}]$ and $[\hat {\boldkappa}]$ introduced with (10).

The prefactor

$[\overline\mu: = {{(2\pi)^3|F_{hkl}|^2} \over {2V^2 k^3 \cos\theta_{\rm B}\sin\theta_{\rm B}}}\eqno(11)]$

depends on the unit-cell volume V and structure factor F_hkl. It has the dimension of an absorption coefficient, i.e. inverse length. The Bragg angle $[\theta_{\rm B}]$ is constant because we consider a fixed reflection hkl and a constant radiation wavenumber k. Both $[\theta_{\rm B}]$ and $[\overline\mu]$ are independent of the sign of the reflection. The outgoing beam direction $[\hat {\bf k}']$ is given by the deflection function

$[\hat {\boldkappa}(\hat {\bf k},t): = \hat {\bf k} - 2\sin\theta_{\rm B} \hat {\boldGamma}(\hat {\bf k},t).\eqno(12)]$

The parametric curve $[\hat {\boldGamma}(\hat {\bf k},t)]$ with $[t\in[0\semi 2\pi]]$ contains all possible scattering directions $[\hat {\bf G}]$ that satisfy the Laue–Bragg condition

$[\hat {\bf G}\hat {\bf k} = \sin\theta_{\rm B}\eqno(13)]$

for an incoming wave direction $[\hat {\bf k}]$ .

To construct $[\hat {\boldGamma}]$ , we choose an orthonormal base $[\{\hat{\bf e}_x,\hat{\bf e}_y,\hat{\bf e}_z\}]$ for the reciprocal-space vectors $[\hat {\bf k}]$ and $[\hat {\bf G}]$ . Note that $[\hat{\bf e}_z]$ is not required to coincide with the plate normal $[\hat{\bf z}]$ (though it does so in our code and our worked-out examples). Choose a rotation matrix $[R_{\bf k}]$ so that $[\hat {\bf k} = R_{\bf k}\hat{\bf e}_z]$ (for readability, we omit carets in subscripts). The circle of possible scattering directions $[\hat {\bf G}]$ can then be written

$[\hat {\boldGamma}(\hat {\bf k},t): = R_{\bf k} \left(\matrix{ \cos\theta_{\rm B}\cos t\cr \cos\theta_{\rm B}\sin t\cr \sin\theta_{\rm B} }\right) = R_{\bf k}\hat {\boldGamma}(\hat{\bf e}_z,t).\eqno(14)]$

It is straightforward to verify that $[\hat {\bf G} = \hat {\boldGamma}(\hat {\bf k},t)]$ , for all t, satisfies (13).

The condition $[\hat {\bf k} = R_{\bf k}\hat{\bf e}_z]$ leaves $[R_{\bf k}]$ underdetermined, allowing for an arbitrary rotation around $[\hat{\bf e}_z]$ . This is irrelevant because the origin of the polar coordinate t is arbitrary, and $[\hat {\boldGamma}(\hat {\bf k},t)]$ only appears under integrals that run from t = 0 to $[2\pi]$ .

2.4. Specializing the distribution of scattering directions

For most mosaic crystals, W is isotropic, i.e. invariant under rotation around the mean block normal $[\hat{\bf H}]$ . Thereby $[W(\hat {\bf G})]$ depends only on $[\hat {\bf G}\hat{\bf H}]$ .

All the following theoretical developments, including the numeric methodology of Section 4, are independent of what isotropic distribution we choose for $[W(\hat {\bf G})]$ . 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 $[W(\hat {\bf G})]$ is a disc-like distribution. As is customary, we will choose a Mises–Fisher (MF) distribution (a Gaussian on the unit sphere),

$[W_{\pm(hkl)}(\hat {\bf G}) = {{1} \over {2\pi\eta^2 C_{\rm MF}(\eta)}}\exp\left[-{{(\hat {\bf G}\mp\hat{\bf H})^2} \over {2\eta^2}}\right]\eqno(15)]$

with the normalization constant $[C_{\rm MF}(\eta) = 1-\exp(-2/\eta^2)]$ . Usually, there is negligible overlap between W_hkl and $[W_{{\overline h}\, {\overline k}\, {\overline l}}]$ so that the sum (8) can be approximated as

$[W(\hat {\bf G}) \simeq {{1} \over {2\pi\eta^2 C_{\rm MF}(\eta)}}\exp\left({{-1+|\hat {\bf G}\hat{\bf H}|} \over {\eta^2}}\right).\eqno(16)]$

A mosaic with $[\hat{\bf H} = \hat{\bf z}]$ shall be called normal oriented. Some consequences of the rotational symmetry around $[\hat{\bf z}]$ are discussed in Appendix C.

In all numeric examples we assume an isotropic, normal oriented mosaic, and we choose $[\hat{\bf e}_z: = \hat{\bf z}]$ . Unless differently stated, the standard deviation is η = 2.5°. The orthographic projection of the circle $[\hat {\boldGamma}(\hat {\bf k},t)]$ into the $[\hat{\bf e}_x,\hat{\bf e}_y]$ plane is an ellipse. Fig. 2 shows examples and puts them in relation to $[W(\hat {\bf G})]$ .

Figure 2
Crystal orientations $[\hat {\bf G}]$ that fulfill the Laue–Bragg condition, projected into the $[\hat{\bf e}_x,\hat{\bf e}_y]$ plane, form ellipses. The two plots have different Bragg angles $[\theta_{\rm B}]$ . Each plot shows ellipses for three different incident angles $[\theta = \arctan(k_z/k_x)]$ , with k_y = 0. The concentric gray discs contain 50% and 90% of all mosaic blocks, assuming a Mises–Fisher distribution $[W(\hat {\bf G})]$ that is centered around $[\hat{\bf H} = \hat{\bf e}_z]$ , with standard deviation η = 2.5°.

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 $[\theta_{\rm in} = \theta_{\rm B}]$ . The first of these parameters is the opacity

$[\tau: = {{d} \over {\sin\theta_{\rm B}}}\left[\overline\mu {\cal G}(0) + \mu_{\rm a}\right],\eqno(17)]$

where $[{\cal G}(0) = 1/[(2\pi)^{1/2}\eta]]$ . The second dimensionless crystal parameter is the relative reflectivity

$[\rho: = {{\overline\mu {\cal G}(0)} \over {\overline\mu {\cal G}(0) + \mu_{\rm a}}}.\eqno(18)]$

These parameters will enter the following derivations through the products

$[d{\overline\mu} = \sin\theta_{\rm B} {{\rho\tau} \over {{\cal G}(0)}}\eqno(19)]$

and

$[d\mu_{\rm a} = \sin\theta_{\rm B} (1-\rho)\tau.\eqno(20)]$

Note that $[\overline\mu,\mu_{\rm a}]$ are not pure material constants but also depend on the wavelength of the used radiation. So, in principle, one could tune $[\tau,\rho]$ to almost arbitrary values by suitable combinations of wavelength and crystal thickness.

3. Discretization in $[{\hat{\bf k}}]$

3.1. Binning

So far, we have assumed that $[I(\hat {\bf k},z)]$ as a function of $[\hat {\bf k}]$ 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 $[I(\hat {\bf k},z)]$ by M histograms I_r(z) with $[r = 1,\ldots,M]$ . Each histogram represents a current that is defined as flux integral over the solid angle $[\Delta\Omega_r]$ ,

$[I_{r}(z): = \textstyle\int\limits_{\Delta\Omega_r}\,{\rm d}^2\hat {\bf k} \ I(\hat {\bf k},z).\eqno(21)]$

Combining this with the definition (2) of the operator B, we get

$[\eqalignno{ (B\circ I)_r & = \textstyle\int\limits_{\Delta\Omega_r}\,{\rm d}^2\hat {\bf k} \ B\circ I(\hat {\bf k},z)&\cr & = \textstyle\int\limits_{\Delta\Omega_r}\,{\rm d}^2\hat {\bf k} \sum\limits_{s = 1}^M\int\limits_{\Delta\Omega_s}\,{\rm d}^2\hat {\bf k}'\mu(\hat {\bf k},\hat {\bf k}')I(\hat {\bf k}',z). &(22)}]$

We assume that μ is a sufficiently smooth function of $[\hat {\bf k}']$ so that it can be drawn in front of the second integral. We obtain

$[(B\circ I)_r \simeq \textstyle\sum\limits_{s = 1}^M B_{rs}I_s(z)\eqno(23)]$

with

$[B_{rs}: = \textstyle\int\limits_{\Delta\Omega_r}\,{\rm d}^2\hat {\bf k} \ \mu(\hat {\bf k},\hat {\bf k}_s).\eqno(24)]$

The attenuation factor, discretized in analogy with (24), is

$[A_r: = \textstyle\sum\limits_{s = 1}^M B_{sr} + \mu_{\rm a}\eqno(25)]$

and the transport equation (1) takes the form

$[(\hat{k}_{r,z}\partial_z + A_r)I_{r}(z) - \textstyle\sum\limits_{s = 1}^M B_{rs}I_{s}(z) = 0.\eqno(26)]$

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.

3.2. 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 M = 2 bins, representing a forward and a backward traveling beam, with $[\hat {k}_z\,\gt\,0]$ and $[\hat {k}_z\,\lt\,0]$ , 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 $[\epsilon_{\rm b}]$ (Section 4.6).

3.3. Diffraction matrix

The integral (24) can be carried out at once since the kernel $[\mu(\hat {\bf k}',\hat {\bf k})]$ , given by (10), contains a delta function. The result is

$[B_{rs} = \overline\mu \cos\theta_{\rm B}\textstyle\int\limits_{0}^{2\pi}\,{\rm d}t W[\hat {\boldGamma}(\hat {\bf k}_s,t)] \left[\hat {\boldkappa}(\hat {\bf k}_s,t)\in\Delta\Omega_r\right]\eqno(27)]$

with the indicator bracket as introduced in Section 2.2.

For given $[\hat {\bf k}_s]$ , and sweeping t, the outgoing directions $[\hat {\boldkappa}(\hat {\bf k}_s,t)]$ 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 $[{\cal O}(M^{1/2})]$ bins, then B has $[{\cal O}(M^{3/2})]$ nonzero entries.

Figure 3
The three bands represent three columns of the reduced reflectivity matrix B_rs, with $[\varphi_s = 0]$ and with three different values of $[\theta_s]$ , shown as a function of $[\varphi_r]$ and $[\theta_r]$ . The Bragg angle is $[\theta_{\rm B}]$ = 70°. There are 80 φ bins from 55° to 85°, and 180 θ bins from −180° to 180°. The dimensionless intensity scale applies for $[\rho = 1]$ and $[\tau = 1]$ .

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.

4. Spectral collocation in the depth coordinate

4.1. Depth rescaling

As we will use Chebyshev polynomials in the depth coordinate, it shall be transformed from [0,d] to the standard range [-1,1]. We therefore introduce the reduced coordinate

$[\zeta: = -1+2z/d \eqno(28)]$

and the transformed histograms

$[J_r(\zeta): = I_r\left(d{{1+\zeta} \over {2}}\right).\eqno(29)]$

The transport equation (26) becomes

$[(2\hat{k}_{r,z}\partial_\zeta + dA_r)J_{r}(\zeta) - \textstyle\sum\limits_{s = 1}^M dB_{rs}J_{s}(\zeta) = 0\eqno(30)]$

with the separated boundary conditions

$[\matrix{ J_r(-1) = I_{{\rm in},r} &{\rm for}\,\,\hat{k}_{r,z}\,\gt\,0,\cr J_r(+1) = 0\hfill &{\rm for}\,\,\hat{k}_{r,z}\,\lt\,0. }\eqno(31)]$

Grids should be constructed such that no bin crosses the equator of the unit sphere.

4.2. Equation system

The equation system (30) consists of M coupled first-order linear differential equations in $[J_r(\zeta)]$ . 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 J_r by polynomials P_r of order N that match J_r in N+1 collocation points $[\zeta_i]$ . We defer the choice of $[\{\zeta_i\}]$ to Section 4.3; until then, we only require $[-1 = \zeta_0\,\lt\,\zeta_1\,\lt\,\ldots\,\lt\,\zeta_{N-1}\,\lt\,\zeta_N = 1]$ . Function values at the collocation points are abbreviated

$[p_{ir}: = P_r(\zeta_i).\eqno(32)]$

These are M(N+1) unknowns. M of them can immediately be read off from the boundary conditions (31):

$[\matrix{ p_{0r} = I_{{\rm in},r} &{\rm for}\,\,\hat{k}_{r,z}\,\gt\,0,\cr p_{Nr} = 0\hfill &{\rm for}\,\,\hat{k}_{r,z}\,\lt\,0. }\eqno(33)]$

The others will be obtained from the transport equation (30). To discretize this differential equation, we replace J_r by P_r, $[\partial_\zeta]$ by a differentiation matrix D, specified by the requirement

$[\textstyle\sum\limits_{j = 0}^N D_{ij} p_{jr} = P'_r(\zeta_i)\eqno(34)]$

[for an introduction to differentiation matrices, see Trefethen (2000)]. The resulting MN equations

$[\textstyle\sum\limits_{j = 0}^N 2 \hat{k}_{r,z} D_{ij} p_{jr} + d A_r p_{ir} - d \sum\limits_{s = 0}^M B_{rs} p_{is} = 0\eqno(35)]$

must hold for all histogram bins ( $[r = 1,\ldots,M]$ ) and in all collocation points ( $[i = 0,\ldots,N]$ ).

We collect all linear operators in the matrix

$[L_{irjs}: = (2 \hat{k}_{r,z} D_{ij}+dA_r\delta_{ij})\delta_{rs} - dB_{rs}\delta_{ij}\eqno(36)]$

with the Kronecker delta $[\delta_{ij}: = [i = j]]$ so that the transport equation (35) becomes simply

$[\forall_{i = 0}^N \forall_{r = 1}^M: \textstyle\sum\limits_{j = 0}^N \sum\limits_{s = 1}^M L_{irjs} p_{js} = 0.\eqno(37)]$

These are M(N+1) equations in M(N+1) variables p_js, 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 B_rs 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 M² entries of matrix B, only $[{\cal O}(M^{3/2})]$ are nonzero. Overall, of the N²M² entries of L, only $[{\cal O}(N^2M + N M^{3/2})]$ are nonzero. This is visualized in Fig. 4.

Figure 4
Visualization of a small part of a small matrix L. Parameters: $[\theta_{\rm B}]$ = 70°, $[\rho = 0.9]$ , $[\tau = 5]$ . Discretization: three collocation points in z; three bins for 55° ≤ θ ≤ 85°; 36 bins for −180° ≤ φ ≤ 180°. Only the 126² entries with 0 ≤ φ ≤ 60° are shown. Since some entries are negative, the figure shows absolute values $[|L_{irjs}|]$ .

In view of the boundary conditions (4), we now distinguish between histogram bins with forward and backward propagation direction, according to the sign of $[\hat{k}_{r,z}]$ . We split the sum over s in (37) accordingly, omitting the zero terms in p_Ns with backward s, and bringing the nonzero terms in p_0s with forward s to the right side:

$[\left(\textstyle\sum\limits_{j = 1}^N\sum\limits_{s}^{\rm forward} +\sum\limits_{j = 0}^{N-1}\sum\limits_{s}^{\rm backward} \right) L_{irjs} p_{js} = - \textstyle\sum\limits_{s}^{\rm forward} L_{ir0s}J_{{\rm in},s}.\eqno(38)]$

This system of M(N+1) inhomogeneous linear equations in MN unknown p_js 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 i = 0, r $[\in]$ forward or i = N, r $[\in]$ backward.

4.3. 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 T_N,

$[\zeta_i: = -\cos{{\pi i} \over {N}}.\eqno(39)]$

They only enter our computation through the corresponding differentiation matrix (34).

This matrix has the outer diagonal entries

$[D_{ij} = {{1+\delta_{i0}+\delta_{iN}} \over {1+\delta_{j0}+\delta_{jN}}} {{(-1)^{i+j}} \over {\zeta_i-\zeta_j}} \quad(i\ne j),\eqno(40)]$

the interior diagonal entries

$[D_{jj} = - {{\zeta_j} \over {2(1-\zeta_j^2)}}\quad(0\,\lt\,j\,\lt\,N)\eqno(41)]$

and the diagonal endpoints

$[D_{00} = -D_{NN} = -{{2N^2+1} \over {6}}.\eqno(42)]$

For a derivation, see e.g. Gottlieb et al. (1984) or Trefethen (2000), but note that our choice of ascending $[\zeta_i]$ has led to a minus sign on the right-hand side of (42).

4.4. Collocation error for the two-ray reference

By the collocation error we understand the error caused by approximating the functions J_r by polynomials P_r. 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 $[\theta_{\rm in} = \theta_{\rm B}]$ so that the Bragg operator (68) is simply $[B = \rho\tau]$ and the total attenuation $[A = \tau]$ .

Fig. 5 shows currents $[J_\pm(\zeta)]$ versus ζ for a moderately thick crystal with realistic attenuation: $[\tau = 5]$ , $[\rho = 0.9]$ . On the linear scale of this plot, the numeric data points and the analytical curves (70) agree perfectly for collocation orders as low as $[N\ge4]$ .

Figure 5
Directional currents (29)

as a function of depth for $[\tau = 5]$ , $[\rho = 0.9]$ . 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 J_tot: = J₊(1)+J_-(-1)+J_a(-1) where J_a 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, N-i. Therefore we generally prefer even N in our computations.

4.5. Collocation error for the full model

We now investigate the collocation error for a two-dimensional $[\hat {\bf k}]$ grid, for representative crystal parameters and for an ideal collimated incoming beam

$[I_{\rm in}(\hat {\bf k}) = \delta^2(\hat{\bf k},\hat{\bf k}_{\rm in}).\eqno(43)]$

We only consider the reflected flux $[I_{\rm refl}(\hat {\bf k}_r) = P_{r}(-1)]$ .

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,

$[\epsilon^{(N)}: = \left\{\textstyle\sum\limits_{r}^{\rm backward} [P^{(N)}_r-P^{(N-\Delta N)}_r]^2\right\}^{1/2}.\eqno(44)]$

In Figs. 7 and 8, we increase N from N_ini = 2 in steps of $[\Delta N = 1]$ ; elsewhere, larger values of the hyperparameters N_ini and $[\Delta N]$ are more convenient (Section 4.6).

Figure 7
Estimates of the collocation error $[\epsilon^{(N)}]$ of the reflected current $[I_{\rm refl}^{(N)}(\hat {\bf k})]$ as a function of the collocation order N, for different solver convergence tolerances. Model parameters $[\tau = 25]$ , $[\rho = 0.9]$ , $[\theta_{\rm B}]$ = 70°. Collimated incoming beam (43)

with $[\theta_{\rm in}]$ = 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 $[\epsilon_{\rm C} = 10^{-5}]$ . 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 $[\epsilon^{(N)}]$ 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 $[\epsilon^{(N)}]$ 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 $[\epsilon^{(N)}]$ lies below a given bound $[\epsilon_{\rm C}]$ (Section 4.6). Fig. 8 shows how the so-determined minimal collocation order N varies with the grid size $[M = M_\theta M_\varphi]$ and with the model parameters $[\theta_{\rm B}]$ , $[\theta_{\rm in}-\theta_{\rm B}]$ , η, ρ, τ. All these parameters are found to be uncritical, except the opacity τ, for which Fig. 8 suggests a possible asymptote $[N\sim\tau^{1/2}]$ . Therefore our computational method should not be applied to extremely thick crystals with $[\tau\gg 100]$ . 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.

4.6. 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 N_ini = 8 and $[\Delta N = 4]$ .

(iii) Auxiliary parameters for the numeric evaluation of the reflection kernel (27), as explained in Section B1: bisection accuracy $[\delta t = 10^{-15}]$ , maximum step size $[\Delta t/\eta = 0.2]$ and a kernel cut-off $[\epsilon_{\rm K} = 10^{-5}]$ .

(iv) The collocation tolerance $[\epsilon_{\rm C}]$ (Section 4.5): the sparse equation solver is called with increasing N until the estimated error $[\epsilon ^{(N)}]$ [equation (44)] falls below $[\epsilon_{\rm C}]$ . We use $[\epsilon_{\rm C} = 10^{-5}]$ .

(v) The solver tolerance $[\epsilon_{\rm S}]$ required by the sparse equation solver (Section 4.7). It should be smaller than $[\epsilon_{\rm C}]$ . We use a value of 10^-7, except when we determine collocation errors as a function of the number of collocation points: Fig. 6 is computed with $[\epsilon_{\rm S} = 3\times 10^{-15}]$ and Fig. 7 with $[\epsilon_{\rm S} = 1\times 10^{-9}]$ .

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 $[\epsilon_{\rm N}]$ to disallow currents below $[-\epsilon_{\rm N}]$ , as may be caused by numeric inaccuracies at weak intensities.

(ii) A stricter limit $[\epsilon_{\rm F}]$ is imposed to the absolute value of the total current (53).

(iii) The cut-off tolerances $[\epsilon_{\rm b}]$ and $[\epsilon_{\rm c}]$ are upper limits for the loss channels I_b and I_c that account for unphysical losses due to the finite $[\hat {\bf k}]$ grid and the deletion of weak matrix elements (Section B2).

All our numeric examples have been checked with $[\epsilon_{\rm N} = 10^{-5}]$ and $[\epsilon_{\rm F} = \epsilon_{\rm b} = \epsilon_{\rm c} = 10^{-4}]$ .

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

5. Results

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

5.2. Total reflectivity

Fig. 9 of Part I showed the total reflectivity R for plates of different opacities as a function of the Bragg angle $[\theta_{\rm B}]$ , 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 $[\Delta R\,\lesssim\, 0.001]$ , our new results are far more accurate and extend over a wider $[\theta_{\rm B}]$ 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 $[\theta_{\rm in} = \theta_{\rm B}]$ ; mosaicity $[\eta]$ = 0.025 rad; nominal opacity $[\tau = 12.5]$ ; relative reflectivity $[\rho = 0.8]$ . In terms of Part I: $[\tau = (\tilde\mu_0 +\tilde\sigma_0)d,]$ $[\rho = \tilde\mu_0/(\tilde\mu_0+\tilde\sigma_0)]$ . Blue symbols with error bars are from the Monte Carlo computation of Part I; red circles are from the present numeric integration.

In the $[\theta_{\rm B}]$ 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 $[\theta_{\rm B}]$ = 90°. Representative results are shown in Fig. 10. With $[\theta_{\rm B}]$ 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 $[\alpha,\beta]$ that fulfill the Bragg condition is almost a circle, and for $[\theta_{\rm in} = \theta_{\rm B}]$ it is concentric with the disc representing $[W(\alpha,\beta)]$ . 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 $[\theta_{\rm in} = \theta_{\rm B}]$ ; mosaicity η = 2.5°.

5.3. Azimuthal distribution

Fig. 11 shows the directional distribution of the transmitted and reflected radiation for three different incoming beam inclinations $[\theta_{\rm in}]$ . Since the coordinate representation does not account for the different bin sizes $[\Delta\Omega_r]$ , this figure does not show currents (flux integrals per bin) but the directional flux $[I(\hat {\bf k}_r,z): = I_r(z)/\Delta\Omega_r]$ .

Figure 11
Directional distribution of the transmitted and reflected flux for three different inclinations $[\theta_{\rm in}]$ of the incoming collimated beam. Crystal parameters: $[\theta_{\rm B}]$ = 70°, η = 2.5°, $[\tau = 5]$ , $[\rho = 0.9]$ .

In transmission, a bright spot shows the attenuated incoming beam. It is least pronounced at $[\theta_{\rm in}= \theta_{\rm B}]$ , 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, $[\theta_{\rm in}\simeq\theta_{\rm B}]$ . For $[\rho = 1]$ and $[\theta_{\rm in} = \theta_{\rm B}]$ , 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 $[\theta_{\rm B}]$ = 70°; (a), (c) as a function of τ for different values of ρ, with $[\theta_{\rm in}= \theta_{\rm B}]$ ; (b), (d) as a function of $[\theta_{\rm in}]$ for different combinations of τ and ρ.

5.4. 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 $[\theta_{\rm in}]$ while maintaining the detector angle at $[2\theta_{\rm B}-\theta_{\rm in}]$ .

It is well known from theory and experiment (Dorner & Kollmar, 1974) that rocking curves are generally wider than the underlying crystallite orientation distribution $[G_\theta]$ . The width increases with increasing opacity τ. This is illustrated by Fig. 13 where rocking curves for a very thin ( $[\tau = 0.1]$ ) and a very thick ( $[\tau = 10]$ ) 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 $[G_\theta]$ .

Figure 13
Rocking curves for (a) a thin, (b) a thick mosaic crystal with $[\tau = 0.1]$ and $[\tau = 10]$ , 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 $[G_\theta]$ , 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.

6. 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 $[|{\bf k}| = k]$ .

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 $[I(\theta,\varphi)]$ 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.

APPENDIX A

Transfer function

A1. Block transfer function

In Part I, the transfer function of a single-crystalline block was derived in three-dimensional wavevector space (equation I,12),

$[\eqalignno{ \mu_{\rm block}({\bf k}',{\bf k}\semi \hat {\bf G})&= \overline\mu k \cos\theta_{\rm B}&\cr &\quad\times\delta({\bf k}\hat {\bf G}-k\sin\theta_{\rm B})&\cr &\quad\times\delta^3({\bf k'}-{\bf k}+2k\sin\theta_{\rm B}\hat {\bf G}). &(45)}]$

The prefactor $[\overline\mu]$ , introduced in (11), agrees with equations I,7 and I,26, with equation 53 in Sears (1997), and equation A7 in Grabcev & Stoica (1980). The three-dimensional delta function in (45) can be decomposed as

$[\delta^3({\bf k}'-{\bf k}) = k^{-2} \delta(k'-k) \delta^2(\hat {\bf k}'-\hat {\bf k})\eqno(46)]$

where $[\delta^2]$ is the delta function on the unit sphere. The scalar delta function is implicit in our redefinition of distribution functions, $[I^{\rm Part\ I}({\bf k}) = \delta(k-k_0)I^{\rm here}(\hat {\bf k})]$ . By comparing the three- and two-dimensional variants of the transport equation we see that it cancels. The factor k^-2 cancels when casting the three-dimensional integral (equation I,38) to our two-dimensional definition (2) of the Bragg operator. A factor 1/k can be drawn out of the first delta function in (45). Altogether we obtain

$[\eqalignno{ \mu_{\rm block}(\hat {\bf k}',\hat {\bf k}\semi \hat {\bf G})& = \overline\mu \cos\theta_{\rm B}&\cr &\quad\times\delta(\hat {\bf k}\hat {\bf G}-\sin\theta_{\rm B})&\cr &\quad\times\delta^2(\hat {\bf k}'-\hat {\bf k}+2\sin\theta_{\rm B}\hat {\bf G}). &(47)}]$

The first delta function enforces the Laue–Bragg condition (13). The second delta function in (47) ensures that the diffracted wave propagates in a direction $[\hat {\bf k}']$ given by the deflection function $[\hat {\boldkappa}(\hat {\bf k},t)]$ , defined in (12).

A2. Total transfer function

To carry out the integral (9), we consider how an integral in $[\hat {\bf G}]$ acts on the first delta function of (47). We parameterize $[\hat {\bf G}]$ in spherical coordinates $[\theta,t]$ with respect to the axis $[\hat {\bf k}]$ , and introduce a test function f to find

$[\eqalignno{ &\textstyle\int\,{\rm d}^2\hat {\bf G} \ \delta(\hat {\bf k}\hat {\bf G}-\sin\theta_{\rm B})f(\hat {\bf G})&\cr & = \textstyle\int\limits_{-\pi/2}^{+\pi/2}\,{\rm d}\theta\cos\theta\int\limits_{0}^{2\pi}\,{\rm d}t \delta(\sin\theta-\sin\theta_{\rm B})f[\hat {\bf G}(\theta,t)]&\cr & = \textstyle\int\limits_{0}^{2\pi}\,{\rm d}t \ f[\hat {\bf G}(\theta_{\rm B},t)]. &(48)}]$

The remaining integral involves a full circle in $[\hat {\bf G}]$ . This circle arises as the intersection of the unit sphere with the plane defined by the Bragg condition (13). For given $[\hat {\bf k}]$ , we will write this circle as $[\hat {\boldGamma}(\hat {\bf k},t): = \hat {\bf G}(\theta_{\rm B},t)]$ . Because it only appears under integrals running from t = 0 to $[2\pi]$ , we can leave the origin in t (the orientation of the spherical coordinates around $[\hat {\bf k}]$ ) unspecified.

A3. Relation to Part I

The above is a simplification in three ways over Part I (equations 25,27): thanks to the parameterization in t, there is no more need to sum over two half $[\hat {\bf G}]$ circles (or ellipses, when referring to the orthographic parameterization of Part I). The singularity $[1/(\cos^2\theta_{\rm B}-\beta^2)^{1/2}]$ in the correction factor h has canceled under the substitution $[{\rm d}\beta = {\rm d}t(\cos^2\theta_{\rm B}-\beta^2)^{1/2}]$ , so that there is no longer a singularity near backscattering. And the β dependence in the second fraction in equation I,27 is gone for good as we no longer use the approximation $[{\rm d}^2{\hat {\bf u}}\simeq{\rm d}\hat u_x{\rm d}\hat u_y]$ implicit in equation I,13.

APPENDIX B

Diffraction matrix

B1. Numeric computation

To compute the matrix element B_rs for given s, we divide the integration domain in (27) in finite intervals. We use bisection to determine points $[t_0 = 0 \,\lt\, t_1 \,\lt\, \ldots \,\lt\, t_{m-1} \,\lt\, t_m = 2\pi]$ such that for $[t\in[t_n,t_{n+1}]]$ all $[\hat {\boldkappa}(\hat {\bf k}_s,t)]$ lie in one and the same bin $[\rho_{sn}]$ . If necessary, intervals are further divided to ensure $[t_{n+1}-t_n\,\lt\,\eta/5]$ . This simplifies (27) to

$[B_{rs} = \overline\mu \cos\theta_{\rm B}\textstyle\sum\limits_{n = 0}^{m-1} [\rho_{sn} = r] \int\limits_{t_n}^{t_{n+1}}\,{\rm d}t W[\hat {\boldGamma}(\hat {\bf k}_s,t)].\eqno(49)]$

Integrals now extend over ranges so small that Gauss–Legendre three-point quadrature is good enough; the accuracy has been ascertained by comparison with higher-order rules. We compute $[\hat {\boldkappa}(\hat {\bf k}_s,t)]$ for the midpoint t = (t_n+1-t_n)/2. From this, we easily deduce $[\rho_{sn}]$ and increment the corresponding $[B_{\rho_{sn}s}]$ . Quite some computational effort can be saved in the special case of an isotropic, normal oriented mosaic, as discussed in Section C2.

To make the matrix B sparser, and thereby speed up the solution of the discretized transport equation, we delete matrix entries that are too tiny to have any consequence for the question under study. Specifically, we set B_rs to zero if it is smaller than $[\epsilon_{\rm K}\max(B_{rs})/M]$ . Our default choice for the cut-off hyperparameter is $[\epsilon_{\rm K} = 10^{-5}]$ (Section 4.6).

B2. Loss channels

If $[\mu_{\rm a} = 0]$ , then the total current in the $[\hat{\bf z}]$ direction is constant,

$[\partial_z\textstyle\sum\limits_{r = 1}^M\hat{k}_{r,z}I_r(z) = 0.\eqno(50)]$

This identity can be used to check the accuracy of a numeric solution.

To maintain (50) even in the presence of non-diffractive losses, we increment M by 1 to allow for a loss channel I_a(z). Numeric experimentation shows that the preferred propagation direction is backward, say $[\hat {\bf k}_{\rm a} = -\hat{\bf z}]$ , which implies that the proper boundary condition is I_a(d) = 0. The diffraction matrix acquires the additional entries

$[\matrix{ B_{{\rm a}s} = \mu_{\rm a}&{\rm for }\ s\ne {\rm a},\cr B_{r{\rm a}} = 0\hfill &{\rm for\ all }\ r, }\eqno(51)]$

and the term $[\mu_{\rm a}]$ appears no longer explicitly in the attenuation factor (25),

$[A_r = \textstyle\sum\limits_{s = 1}^M B_{sr}.\eqno(52)]$

This approach comes to fruition with two more loss channels, I_b(z) and I_c(z), which account for errors introduced by two approximations that reduce the number of nonzero matrix entries B_rs:

One approximation consists of choosing a grid that does not cover the full unit sphere (Section 3.2). It makes (22) inexact because for some s there is a finite probability B_bs of diffraction towards some $[\hat {\bf k}']$ that is not in the grid.

The other approximation is the zeroing of matrix entries B_rs that are too tiny to have any practical consequence (Section B1). To assess the overall error made by this approximation, we set B_cs to the sum of all deleted entries B_rs for given s.

As per (51), there is no scattering out of these channels: B_rb = B_rc = 0. Starting from I_b(d) = I_c(d) = 0, intensity accumulates with decreasing z. The total approximation losses can be read off from I_b(0) and I_c(0). If they are below tolerances $[\epsilon_{\rm b}, \epsilon_{\rm c}]$ (Section 4.6), then one can be sure that approximations made by restricting the grid and by zeroing some B_rs are unproblematic.

Altogether, the total current, including all loss channels, is

$[\textstyle\sum\limits_{r}^{\rm {forward}} \hat{k}_{r,z}[I^{\rm in}_r(0)-I_r(d)] -\sum\limits_{r}^{\rm {backward}} \hat{k}_{r,z}I_r(0) -\sum\limits_{l}^{\rm {a,b,c}} I_l = 0.\eqno(53)]$

Allowing for some numeric inaccuracy, the absolute value of the left-hand side is requested to stay below a tolerance $[\epsilon_{\rm F}]$ (Section 4.6).

APPENDIX C

Isotropic, normal oriented mosaic

C1. Spherical coordinates

Isotropic, normal oriented mosaics (defined in Section 2.4) have a rotational symmetry around $[\hat{\bf z}]$ that allows us to simplify and accelerate some computations. For given $[\hat{\bf e}_z = \hat{\bf z} = \hat {\bf H}]$ , we choose orthonormal vectors $[\hat{\bf e}_x,\hat{\bf e}_y]$ . For a given reciprocal-space direction $[\hat {\bf k}]$ , we define spherical coordinates $[\theta_{\bf k},\varphi_{\bf k}]$ with respect to the base $[\hat{\bf e}_x,\hat{\bf e}_y,\hat{\bf e}_z]$ :

$[\hat {\bf k} =: \left(\matrix{ \cos\theta_{\bf k}\cos\varphi_{\bf k}\cr \cos\theta_{\bf k}\sin\varphi_{\bf k}\cr \sin\theta_{\bf k} }\right).\eqno(54)]$

We choose the rotation matrix in equation (14) as $[R_{\bf k} = R_{\bf z}(\varphi_{\bf k})R_{\bf y}(\pi/2-\theta_{\bf k})]$ , where $[R_{\bf z},R_{\bf y}]$ are rotations around the z,y axes. It is easily verified that $[\hat {\bf k} = R_{\bf k}\hat{\bf z}]$ . The deflection function (12) is

$[\eqalignno{ \hat {\boldkappa}(\hat {\bf k},t) & = \hat {\bf k} - 2\sin\theta_{\rm B}\hat {\boldGamma}(\hat {\bf k},t)&\cr & = R_{\bf k}\left[ \hat{\bf z} - 2\sin\theta_{\rm B} \hat {\boldGamma}(\hat{\bf z},t)\right]&\cr & = R_{\bf z}(\varphi_{\bf k})\left[R_{\bf y}(\pi/2-\theta_{\bf k})\hat {\boldkappa}(\hat{\bf z},t)\right]. &(55)}]$

With the last line, we obtained a factorization of the $[\varphi_{\bf k}]$ and $[\theta_{\bf k}]$ dependence of $[\hat {\boldkappa}]$ .

C2. Computing the diffraction matrix

The algorithm to compute the diffraction matrix B_rs described in Section B1 has an outer loop that runs over the bin index $[s = 1,\ldots,M]$ . If for an isotropic, normal oriented mosaic bins are chosen on a rectangular grid in the spherical coordinates $[\theta_{\bf k},\varphi_{\bf k}]$ , then that loop must only be executed for $[s = 1,\ldots,M_\theta]$ at one fixed $[\varphi_{\bf k}]$ . In a second step, results are then transcribed to all other $[\varphi_{\bf k}]$ ; the outgoing bin indices r are translated accordingly. This is particularly easy if the φ grid extends over $[2\pi]$ so that a periodic boundary condition applies. A further factor of 2 in computational effort can be saved by using the symmetry of $[\hat {\boldkappa}]$ and W under a change of sign of t.

Given the circulancy of B in φ, one could also decouple equations by Fourier transform. However, this would come at an expense in sparsity, and it would not help for future applications with non-normal oriented mosaics. Therefore, we have not explored this idea any further.

C3. Azimuthal average

If we are not interested in the azimuthal distribution of the radiation reflected or transmitted from an isotropic, normal oriented mosaic then the rotational symmetry allows us to integrate out the coordinate φ so that only the θ dependence of the flux

$[I(\theta,z): = \textstyle\int\limits_0^{2\pi}\,{\rm d}\varphi I[\hat {\bf k}(\theta,\varphi),z]\eqno(56)]$

is studied further. To integrate over the transfer function (10), we note that $[W[\hat {\boldGamma}(\hat {\bf k},t)]]$ is independent of $[\varphi_{\bf k}]$ . We find the deflection operator (2)

$[B\circ I(\theta): = \textstyle\int\limits_{0}^{2\pi}\,{\rm d}\varphi B\circ I[\hat {\bf k}(\theta,\varphi)] = \int\limits_{-\pi/2}^{\pi/2}\,{\rm d}\theta'\mu(\theta,\theta') I(\theta')\eqno(57)]$

with the transfer function

$[\mu(\theta',\theta): = \overline\mu \cos\theta_{\rm B} \textstyle\int\limits_{0}^{2\pi}\,{\rm d}t W[\hat {\boldGamma}(\theta,t)] \delta[\theta'-\Theta(\theta,t)]\eqno(58)]$

and the polar deflection function

$[\Theta(\theta,t): = \arcsin\left[\hat{\bf z} R_{\bf y}(\pi/2-\theta_{\bf k})\hat {\boldkappa}(\hat{\bf z},t) \right].\eqno(59)]$

We integrate (58) to obtain the attenuation operator

$[A(\hat {\bf k}) = \overline\mu w(\hat {\bf k}) + \mu_{\rm a}\eqno(60)]$

with the dimensionless total deflection probability

$[w(\hat {\bf k}): = \cos\theta_{\rm B}\textstyle\int\limits_{0}^{2\pi}\,{\rm d}tW[\hat {\boldGamma}(\hat {\bf k},t)].\eqno(61)]$

APPENDIX D

The two-ray model

D1. Derivation

If the reflected polar angle (59) is expanded in $[\theta-\theta_{\rm B}]$ and t, then the lowest order is just [in accord with equation I,64, but note that the next order has been corrected in Wuttke (2020 )]

$[\Theta(\theta,t) \,\doteq\, 2\theta_{\rm B} - \theta.\eqno(62)]$

Accordingly, (58) and (61) are simplified to

$[\mu(\theta',\theta) \,\doteq\, \overline\mu w(\theta) \delta(\theta'+\theta-2\theta_{\rm B})\eqno(63)]$

and

$[w(\theta) = \cos\theta_{\rm B} \textstyle\int\limits_{0}^{2\pi}\,{\rm d}t W[\hat {\boldGamma}(\theta,t)].\eqno(64)]$

An incoming collimated beam, through arbitrarily many reflections, will propagate forward along $[\theta_{\rm in}]$ and backward along $[2\theta_{\rm B}-\theta_{\rm in}]$ , as described by the two-ray Darwin–Hamilton equations.

Realistic mosaic distributions are narrow; in radians, $[\eta\ll 1]$ . This justifies the above truncation, and motivates a series expansion of $[\hat {\boldGamma}]$ in the small arguments $[|\theta|-\theta_{\rm B}]$ and t, for use in (63). In lowest order, we find

$[(\hat {\bf G}\pm\hat{\bf H})^2 \,\doteq\, {(|\theta|-\theta_{\rm B})}^2 + \cos^2\theta_{\rm B}t^2.\eqno(65)]$

Specifically, with the Mises–Fisher distribution (16), we can carry out the integral (64) to find

$[w(\theta) \,\doteq\, {\cal G}(|\theta|-\theta_{\rm B})\eqno(66)]$

with the normalized univariate Gaussian

$[{\cal G}(\delta): = {{1} \over {(2\pi )^{1/2}\eta}}\exp\left(-{{\delta^2} \over {2\eta^2}}\right).\eqno(67)]$

In the two-ray notation of Part I, the Bragg operator (57) becomes

$[B\circ I_\pm = \overline\mu {\cal G}(|\theta|-\theta_{\rm B}) I_\mp,\eqno(68)]$

in agreement with the standard treatment of the planar Darwin–Hamilton equations (Zachariasen, 1945 $[Zachariasen, W. H. (1945). Theory of X-ray Diffraction in Crystals. New York: Wiley.]$ , equation 4.19; Sears, 1989, equation 5.2.70).

D2. Analytical solution

The two-ray boundary problem with plate geometry has been solved in full generality by Sears (1997). Since our notation deviates from Sears' in potentially confusing ways, translations are given in Table 1. According to (68), the Bragg cross sections are the same for both rays: $[B: = B(\theta_+) = B(\theta_-)]$ . The abbreviations from Sears' equation 12 are in our notation

$[\eqalignno{ p&: = {{dA} \over {2}} \left({{1} \over {\sin\theta_+}} + {{1} \over {\sin\theta_-}} \right),&\cr q&: = {{dA} \over {2}} \left({{1} \over {\sin\theta_+}} - {{1} \over {\sin\theta_-}} \right),&\cr r&: = {\left(p^2 - {{dB} \over {\sin\theta_+}}{{dB} \over {\sin\theta_-}} \right)}^{1/2}. &(69)}]$

The resulting directional currents $[I_\pm(z)]$ are then given by Sears' equation 15.

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

Variable	Sears (1997)	This work
Crystal thickness	d	d
Depth variable	$[z = 0\ldots d]$	$[z = 0\ldots d]$
Rescaled depth variable		$[\zeta = -1\ldots+1]$
Bragg angle	θ	$[\theta_{\rm B}]$
Forward current	I	I₊
Backward current		I_-
Forward beam polar angle	φ	$[\theta_+\equiv\theta_{\rm in}]$
Backward beam polar angle	$[\varphi']$	$[\theta_-\equiv 2\theta_{\rm B}-\theta_{\rm in}]$
Bragg cross section	σ	B
Bragg cross section, reduced, forward beam	$[b\equiv d\sigma/\sin\varphi]$	$[dB/\sin\theta_+]$
Bragg cross section, reduced, backward beam	$[b\zeta\equiv d\sigma/\sin\varphi']$	$[dB/\sin\theta_-]$
Loss cross section	μ	$[\mu_{\rm a}]$
Total attenuation cross section	$[\sigma+\mu]$	$[A = B+\mu_{\rm a}]$
Total attenuation cross section, reduced, forward beam	a+b	$[dA/\sin\theta_+]$
Total attenuation cross section, reduced, backward beam	$[(a+b)\zeta]$	$[dA/\sin\theta_-]$

We now specialize to the symmetric case $[\theta_{\rm in}= \theta_{\rm B}]$ , hence $[\theta_+ = \theta_-]$ , as addressed in Fig. 5. With our abbreviations τ and ρ from Section 2.5 and with the further abbreviation $[\omega: = ({1-\rho^2} )^{1/2}]$ , (69) reduces to $[p = \tau]$ , q = 0 and $[r = \tau\omega]$ . Sears' equation 15 yields

$[\eqalignno{ J_+(\zeta) &= {{\omega\cosh[\tau\omega(1-\zeta)/2]+\sinh[\tau\omega(1-\zeta)/2]}\over{\omega\cosh(\tau\omega) + \sinh(\tau\omega)}}I_{\rm in},&\cr J_-(\zeta) &= {{\rho\sinh[\tau\omega(1-\zeta)/2]}\over {\omega\cosh(\tau\omega) + \sinh(\tau\omega)}}I_{\rm in}.&(70)} ]$

To discuss rocking curves, we need the transmittivity and reflectivity for arbitrary $[\theta_{\rm in}]$ (Sears' equation 16)]:

$[\eqalignno{ T &: = {{J_+(+1)} \over {I_{\rm in}}} = \exp(-q){{r} \over {r\cosh r+ p \sinh r}},&\cr R &: = {{J_-(-1)} \over {I_{\rm in}}} = {{dB} \over {\sin\theta_+}}{{\sinh r} \over {r\cosh r+ p \sinh r}}. &(71)}]$

For $[\theta_{\rm in} = \theta_{\rm B}]$ , we obtain by specializing either (70) or (71)

$[\eqalignno{ T & = {{\omega} \over {\omega\cosh(\tau\omega)+ \sinh(\tau\omega)}},&\cr R & = {{\rho\sinh(\tau\omega)} \over {\omega\cosh(\tau\omega)+ \sinh(\tau\omega)}}. &(72)}]$

Supporting information

Open-source software MultiBragg as used to generate Figs 3-13. DOI: https://doi.org/10.1107/S2053273320002065/ae5082sup1.tgz

Acknowledgements

We thank an anonymous reviewer for detailed and very helpful suggestions.

References

Bavier, E., Hoemmen, M., Rajamanickam, S. & Thornquist, H. (2012). Sci. Program. 20, 241–255. Google Scholar
Canuto, C., Quarteroni, A., Hussaini, M. Y. & Zang, T. (1988). Spectral Methods in Fluid Mechanics. New York: Springer. Google Scholar
Darwin, C. G. (1922). London Edinb. Dubl. Philos. Mag. J. Sci. 43, 800–829. CrossRef CAS Google Scholar
Dorner, B. & Kollmar, A. (1974). J. Appl. Cryst. 7, 38–41. CrossRef CAS IUCr Journals Web of Science Google Scholar
Driscoll, T. A. & Hale, N. (2016). IMA J. Numer. Anal. 36, 108–132. Google Scholar
Frick, B., Bordallo, H. N., Seydel, T., Barthélémy, J.-F., Thomas, M., Bazzoli, D. & Schober, H. (2006). Physica B, 385–386, 1101–1103. Web of Science CrossRef CAS Google Scholar
Gottlieb, D., Hussaini, M. Y. & Orszag, S. A. (1984). In Spectral Methods for Partial Differential Equations, edited by R. G. Voigt, D. Gottlieb & M. Y. Hussaini. Philadelphia: SIAM. Google Scholar
Grabcev, B. & Stoica, A. D. (1980). Acta Cryst. A36, 510–519. CrossRef CAS IUCr Journals Web of Science Google Scholar
Hamilton, W. C. (1957). Acta Cryst. 10, 629–634. CrossRef CAS IUCr Journals Web of Science Google Scholar
Hennig, M., Frick, B. & Seydel, T. (2011). J. Appl. Cryst. 44, 467–472. Web of Science CrossRef CAS IUCr Journals Google Scholar
Heroux, M., Bartlett, R., Howle, V., Hoekstra, R., Hu, J., Kolda, T., Lehoucq, R., Long, K., Pawlowski, R., Phipps, E., Salinger, A., Thornquist, H., Tuminaro, R., Willenbring, J. & Williams, A. (2003). An Overview of Trilinos. Report SAND2003-2927. Albuquerque: Sandia National Laboratories. Google Scholar
Meyer, A., Dimeo, R. M., Gehring, P. M. & Neumann, D. A. (2003). Rev. Sci. Instrum. 74, 2759–2777. Web of Science CrossRef CAS Google Scholar
Moler, C. & Van Loan, C. (1978). SIAM Rev. 20, 801–836. CrossRef Web of Science Google Scholar
Moler, C. & Van Loan, C. (2003). SIAM Rev. 45, 3–49. Web of Science CrossRef Google Scholar
Ohmasa, Y. & Chiba, A. (2018). Acta Cryst. A74, 681–698. Web of Science CrossRef IUCr Journals Google Scholar
Ohmasa, Y., Shimomura, S. & Chiba, A. (2016). J. Appl. Cryst. 49, 835–844. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sala, M. & Heroux, M. (2005). Robust Algebraic Preconditioners with IFPACK 3.0. Report SAND-0662. Albuquerque: Sandia National Laboratories. Google Scholar
Schneider, J. R. (1974). J. Appl. Cryst. 7, 541–546. CrossRef IUCr Journals Web of Science Google Scholar
Sears, V. F. (1989). Neutron Optics. Oxford University Press. Google Scholar
Sears, V. F. (1997). Acta Cryst. A53, 35–45. CrossRef CAS Web of Science IUCr Journals Google Scholar
Trefethen, L. N. (2000). Spectral Methods in MATLAB. Philadelphia: SIAM. Google Scholar
Wuttke, J. (2014a). Acta Cryst. A70, 429–440. Web of Science CrossRef IUCr Journals Google Scholar
Wuttke, J. (2014b). J. Phys. A Math. Theor. 47, 215203. Web of Science CrossRef Google Scholar
Wuttke, J. (2020). Acta Cryst. A76, 215. Web of Science CrossRef IUCr Journals Google Scholar
Wuttke, J., Budwig, A., Drochner, M., Kämmerling, H., Kayser, F.-J., Kleines, H., Ossovyi, V., Pardo, L. C., Prager, M., Richter, D., Schneider, G. J., Schneider, H. & Staringer, S. (2012). Rev. Sci. Instrum. 83, 075109. Web of Science CrossRef PubMed Google Scholar
Xu, K. & Hale, N. (2016). IMA J. Numer. Anal. 36, 618–632. Web of Science CrossRef Google Scholar
Zachariasen, W. H. (1945). Theory of X-ray Diffraction in Crystals. New York: Wiley. Google Scholar

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

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 76| Part 3| May 2020| Pages 376-389

https://doi.org/10.1107/S2053273320002065

Open

access

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Format		BIBTeX
		EndNote
		RefMan
		Refer
		Medline
		CIF
		SGML
		Plain Text
		Text

Search IUCr Journals		doi		Advanced search
Author		volume	page

research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

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

1. Introduction

2. The mathematical model

2.1. Crystal model and current distribution

2.2. Transport equation and boundary conditions

2.3. Reflection kernel

2.4. Specializing the distribution of scattering directions

2.5. Parameterization

3. Discretization in

3.1. Binning

3.2. Grids in zero, one, two dimensions

3.3. Diffraction matrix

4. Spectral collocation in the depth coordinate

4.1. Depth rescaling

4.2. Equation system

4.3. Collocation points and differentiation matrix

4.4. Collocation error for the two-ray reference

4.5. Collocation error for the full model

4.6. Hyperparameters and a posteriori tolerances

4.7. Numeric tools

5. Results

5.1. Open-source code MultiBragg

5.2. Total reflectivity

5.3. Azimuthal distribution

5.4. Rocking curves

6. Conclusions

APPENDIX A

Transfer function

A1. Block transfer function

A2. Total transfer function

A3. Relation to Part I

APPENDIX B

Diffraction matrix

B1. Numeric computation

B2. Loss channels

APPENDIX C

Isotropic, normal oriented mosaic

C1. Spherical coordinates

C2. Computing the diffraction matrix

C3. Azimuthal average

APPENDIX D

The two-ray model

D1. Derivation

D2. Analytical solution

Supporting information

Acknowledgements

References

research papers

3. Discretization in $[{\hat{\bf k}}]$