Formulation of perfect-crystal diffraction from Takagi–Taupin equations: numerical implementation in the crystalpy library

Guigay, J.-P.; Sanchez del Rio, M.

doi:10.1107/S160057752400924X

research papers

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 31| Part 6| November 2024| Pages 1469-1480

https://doi.org/10.1107/S160057752400924X

Open

access

Formulation of perfect-crystal diffraction from Takagi–Taupin equations: numerical implementation in the crystalpy library

Jean-Pierre Guigay ^a ^* and Manuel Sanchez del Rio ^a ^*

^aEuropean Synchrotron Radiation Facility, 71 Avenue des Martyrs, F-38000 Grenoble, France
^*Correspondence e-mail: [email protected], [email protected]

Edited by H. Tolentino, Brazilian Synchrotron Light Laboratory, Brazil (Received 21 June 2024; accepted 20 September 2024; online 29 October 2024)

The Takagi–Taupin equations are solved in their simplest form (zero deformation) to obtain the Bragg-diffracted and transmitted complex amplitudes. The case of plane-parallel crystal plates is discussed using a matrix model. The equations are implemented in an open-source Python library crystalpy adapted for numerical applications such as crystal reflectivity calculations and ray tracing.

Keywords: X-ray diffraction; perfect crystal; dynamical theory of diffraction; crystalpy software package; crystal optics.

1. Introduction

Almost every synchrotron radiation beamline operating with hard X-rays makes use of perfect crystals. Most beamlines use a double-crystal monochromator with flat crystals. Multiple reflections are used for high resolution (Ishikawa et al., 2005 ; Shvyd'ko, 2004 ). Curved crystals are used in reflection [polychromators for dispersive X-ray absorption spectroscopy (Tolentino et al., 1988 )] or in transmission [single- (Suortti et al., 1993 ) or double-crystal Laue monochromators (Ren et al., 1999 )]. Plane crystals plates are used to change the polarization state of X-rays (Bouchenoire et al., 2003 ; Detlefs et al., 2012 ). In addition, crystal analyzers are used in most spectroscopy beamlines [see, for example, Rovezzi et al. (2017 )].

Beamline simulation tools used for the design, optimization and commissioning of synchrotron instrumentation implement in software the equations to calculate the reflectivity of perfect crystals. The theory of diffraction [see Authier (2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction. Oxford University Press.]$ ) for a complete reference] is the basis of all numeric implementations.

There are many simulation tools implementing the equations of the dynamical theory in different forms. This variate scenario is even more complex if we consider that the calculation of the crystal structure factor, which is an essential ingredient for calculating diffracted amplitudes and intensities, is obtained from tabulated scattering functions of multiple origins. A wide collection of software methods and tools can be found even in a single application, such as the OASYS suite (Rebuffi & Sanchez del Rio, 2017 ), which provides multiple solutions for calculating diffraction profiles of crystals [e.g. INPRO (https://github.com/oasys-kit/xoppy_external_codes/tree/master/src/INPRO), CRYSTAL (Sanchez del Rio et al., 2015 ), X-RAY Server (Stepanov, 2004 )], as well as beamline simulation tools [based on the ray-tracing code SHADOW (Sanchez del Rio et al., 2011 )] and physical wave-optics simulations with SRW (Chubar & Elleaume, 1998 ; Sutter et al., 2014 ). Most ray-tracing codes used for synchrotron applications incorporate models for crystal diffraction, like SHADOW (Sanchez del Rio et al., 2011), RAY (Schäfers, 2008 ; Baumgärtel et al., 2019 ), XRT (Chernikov & Klementiev, 2017 ; Klementiev & Chernikov, 2023 ). This scenario has inherited decades of advancements and has witnessed the evolution of several generations of synchrotron radiation sources. Our research aims to tackle this challenge by consolidating the resources for crystal diffraction calculations. We have two primary objectives: deducing the equations governing crystal reflectivity from first principles and integrating them into a thoroughly documented open-source software library.

The Takagi–Taupin equations are a powerful tool in Bragg diffraction by deformed crystals for diverse forms of the incident monochromatic wave. They are applied here to the simple particular case of plane parallel crystals and plane incident waves. We derive results found in the conventional dynamical theory described in textbooks (Zachariasen, 1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ; Pinsker, 1978 ; Authier, 2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction. Oxford University Press.]$ ). Therefore, there are no new physical results in the present paper. However, the method presented here is mathematically well defined and simple. It is general in the sense that it deals directly with absorbing crystals. We believe that it represents a valuable and useful shortcut to the conventional method.

In Section 2 we derive the Takagi–Taupin (TT) equations (Takagi, 1962 ; Taupin, 1964 ; Taupin, 1967 ). In Section 3 we solve the TT equations for a plane undeformed-crystal. Given known complex amplitudes at the entrance surface, the complex amplitudes along the incident and diffracted directions at the back surface are calculated via a transfer matrix (Section 3.2). For the Laue case, the transfer matrix is directly used to compute the diffracted and forward-diffracted (or transmitted) complex amplitudes (Section 3.4). For the Bragg case (Section 3.5) the transfer matrix is used to obtain the scattering matrix, which gives the diffracted and transmitted complex amplitudes. Section 4 is dedicated to the software implementation of the library crystalpy. A final summary and conclusions are in Section 5.

2. Takagi–Taupin equations

The scalar time-independent X-ray wave equation in a perfect crystal is the Helmholtz equation,

$[\Delta\Psi+k^{2}\,\big[1+\chi({\bf{r}})\big]\,\Psi({\bf{r}}) = 0, \eqno(1)]$

where $[\Psi({\bf{r}})]$ is the wavefunction, k = 2π/λ, with λ the wavelength, $[\chi({\bf{r}})]$ is the electric susceptibility [refractive index n = (1 + χ)^1/2] that can be expanded in a Fourier series,

$[\chi({\bf{r}}) = \sum_{{\bf{h}}}\chi_{h}\exp(i\,{\bf{h}}.{\bf{r}}), \eqno(2)]$

where the sum goes over all reciprocal lattice vectors h with (hkl) Miller indices. The spacing of the (hkl) reflection is d_hkl = 2π/h, where h = $[|{\bf{h}}|]$ .

Let us consider an incident plane wave $[\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right)]$ in vacuum. Its wavevector $[{\bf{k}}_{0}]$ , with modulus k = $[|{\bf{k}}_{0}|]$ , is close to the Bragg condition for the diffraction vector h. In the `two-beams case' of Bragg diffraction, considered in this paper, we define $[{\bf{k}}_{h}]$ $[\equiv]$ $[{\bf{k}}_{0}+{\bf{h}}]$ , of modulus k_h = $[|{\bf{k}}_{h}|]$ . In general, the direction of $[{\bf{k}}_{h}]$ does not correspond to the Bragg-diffracted wavevector in vacuum, and k_h is slightly different from k. The deviation from the exact Bragg position is expressed by the parameter α ( $[\alpha\ \ll\ 1]$ ), defined as

$[\alpha = {{k^{2}-k_{h}^{2}} \over {k^{2}}} = {{k^{2}-|{\bf{k}}_{0}+{\bf{h}}|^{ 2}} \over {k^{2}}} = -{{{\bf{h}}^{2}+2{\bf{k}}_{0}.{\bf{h}}} \over {k^{2}}}. \eqno(3)]$

The wavefield $[\Psi({\bf{r}})]$ in the crystal is set empirically as the sum of `two modulated plane waves',

$[\Psi({\bf{r}}) = D_{0}({\bf{r}})\exp\left({i\,{\bf{k}}_{0}.{\bf{r}}}\right)+D_{h}({\bf{r}})\exp\left({i\,{\bf{k}}_{h}.{\bf{r}}}\right), \eqno(4)]$

in which the amplitudes $[D_{0,h}({\bf{r}})]$ are considered as `slowly varying functions', thus neglecting their second-order derivatives in $[\Delta\Psi({\bf{r}})]$ ,

$[\eqalign{\Delta\big[D_{0,h}&({\bf{r}})\exp\left(i\,{\bf{k}}_{0,h}.{\bf{r}}\right)\big] = \cr& \big[2i\,{\bf{k}}_{0, h}.\nabla D_{0,h}+(k^{2}-k^{2}_{0,h})D_{0,h}\big] \exp\left(i\,{\bf{k}}_{0,h}.{\bf{r}}\right),}]$

thus giving

$[\eqalignno{\Delta\Psi+k^{2}\Psi = {}&\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right) \big[2i\,{\bf{k}}_{0}.\nabla D_{0}\big] \cr& + \exp\left(i\,{\bf{k}}_{h}.{\bf{r}}\right) \big[2i\,{\bf{k}}_{h}.\nabla D_{ h}+(k^{2}-k_{h}^{2})\,D_{h}\big]. &(5)}]$

In the product $[\chi({\bf{r}})\Psi({\bf{r}})]$ , using the equations (2) and (4), we write separately the terms containing either $[\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right)]$ or $[\exp\left(i\,{\bf{k}}_{h}.{\bf{r}}\right)]$ , and do not consider the other terms,

$[\eqalignno{ \chi\Psi = {}& \big[\chi_{0}D_{0}+\chi_{-h}D_{h}\big] \exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right) \cr& + \big[\chi_{h}D_{0}+\chi_{0}D_{h}\big] \exp\left(i\,{\bf{k}}_{h}.{\bf{r}}\right) + \ldots. &(6) }]$

Inserting equations (5) and (6) into (1), we obtain the TT equations,

$[\eqalignno{ 2i\,{\bf{k}}_{0}.\nabla D_{0}+\chi_{0}k^{2}D_{0}+\chi_{-h}k^{2}D_{h} & = 0, & (7a) \cr 2i\,{\bf{k}}_{h}.\nabla D_{h}+(\alpha+\chi_{0})\,k^{2}D_{h}+\chi_{h}k^{2}D_{0} & = 0. & (7b) }]$

We can use the oblique coordinates (s₀, s_h) in the diffraction plane (the plane containing $[{\bf{k}}_{0}]$ and h, as well as $[{\bf{k}}_{h}]$ ), with origin O on the crystal surface and unit vectors $[\hat{{\bf{s}}}_{0}]$ and $[\hat{{\bf{s}}}_{h}]$ along $[{\bf{k}}_{0}]$ and $[{\bf{k}}_{h}]$ , respectively. A generic spatial position $[{\bf{r}}]$ = (s₀,s_h,s_t) should include a third coordinate s_t along an axis $[\hat{{\bf{s}}}_{t}]$ non-coplanar with $[{\bf{k}}_{0}]$ and $[{\bf{k}}_{h}]$ . We can choose $[\hat{{\bf{s}}}_{t}]$ to lie on the crystal entrance surface and be perpendicular to the intersection line of the diffraction plane and the crystal surface. The chosen direction of $[\hat{{\bf{s}}}_{t}]$ implies $[{\bf{n}}.\hat{{\bf{s}}}_{t}]$ = 0. The equation of the crystal surface is γ₀s₀ + γ_hs_h = 0, with $[\gamma_{0,h}]$ = $[{\bf{n}}.\hat{{\bf{s}}}_{0,h}]$ $[\equiv]$ $[\cos(\theta_{0,h})]$ the director cosines with respect to n, the unit inward normal vector to the entrance plane surface.

The relation ds₀ = $[\nabla s_{0}.d{\bf{r}}]$ = $[\nabla s_{0}.[ds_{0}\hat{{\bf{s}}}_{0}+ds_{h} \hat{{\bf{s}}}_{h}+ds_{t}\hat{{\bf{s}}}_{t}]]$ implies $[\nabla s_{0}.\hat{{\bf{s}}}_{0}]$ = 1 and $[\nabla s_{0}.\hat{{\bf{s}}}_{h,t}]$ = 0. Similarly, $[\nabla s_{h}.\hat{{\bf{s}}}_{h}]$ = 1 and $[\nabla s_{h}.\hat{{\bf{s}}}_{0,t}]$ = 0. Therefore,

$[\eqalignno{ \hat{{\bf{s}}}_{0}.\nabla D = \hat{{\bf{s}}}_{0}.\left[{{ \partial D} \over {\partial s_{0}}}\nabla s_{0}+{{\partial D} \over {\partial s_{h}}} \nabla s_{h}+{{\partial D} \over {\partial s_{t}}}\nabla s_{t}\right] & = {{ \partial D} \over {\partial s_{0}}}, & (8a) \cr \hat{{\bf{s}}}_{h}\nabla D & = {{\partial D} \over {\partial s_{h}}}. & (8b) }]$

Using the approximation¹

$[{\bf{k}}_{h}.\nabla D_{h} = k\left({1-\alpha}\right)^{1/2} \, {{\partial D_{h}}\over{\partial s_{h}}} \simeq k{{\partial D_{h}} \over {\partial s_{h}}}, \eqno(9)]$

we obtain from equation (7),

$[\eqalignno{ {{\partial D_{0}} \over {\partial s_{0}}} & = {{ik}\over{2}} \big[\chi_{0}D _{0}+\chi_{-h}D_{h}\big] = iu_{0}D_{0}+iu_{-h}D_{h}, & (10a) \cr {{\partial D_{h}} \over {\partial s_{h}}} & = {{ik}\over{2}} \big[(\chi_{0}+\alpha)\,D_{h}+\chi_{h}D_{0}\big] = i(u_{0}+\alpha^{\prime})\,D_{h}+iu_{h}D_{0}, & (10b)}]$

where we used the notation

$[\eqalignno{ u_{0,h,-h} & = {{k}\over{2}}\chi_{0,h,-h}, & (11a) \cr \alpha^{\prime} & = {{k}\over{2}}\alpha. & (11b) }]$

Our formulation is written for the case of σ-polarization. For the case of π-polarization, it is sufficient to replace χ_h and χ_−h with Cχ_h and Cχ_−h with C = $[\cos(2\theta_{\rm{B}})]$ . An equivalent form of the TT equations (10) is obtained in Appendix A, using oblique axes along the directions of the geometrical Bragg law, and applying a crystal rotation.

3. Solutions of TT equations for a plane wave incident on a crystal with plane entrance surface

It is interesting to consider first the effects of refraction and absorption without Bragg diffraction. Setting u_−h = 0 in equations (10a), we obtain the following equation for the refracted amplitude,

$[{{\partial D_{0}^{\,\rm{ref}}} \over {\partial s_{0}}} = iu_{0}D_{0}^{\,\rm{ref}}. \eqno(12)]$

Its solution satisfying the boundary conditions $[D_{0}^{\,\rm{ref}}]$ = 1 for γ₀s₀ + γ_hs_h = 0 (equation of the crystal surface) is $[D_{0}^{\,\rm{ref}}]$ = $[\exp\big[iu_{0}(s_{0}+s_{h}\gamma_{h}/\gamma_{0})\big]]$ = exp(iu₀s) where

$[s = s_{0}+{{s_{h}} \over {b}}\semi\qquad b = {{\gamma_{0}} \over {\gamma_{h}}}. \eqno(13)]$

We now consider the solutions of the equations (10) depending on the single variable² s, which means ∂D₀/∂s₀ = $[D_0^{\,\prime}(s)]$ and ∂D_h/∂s_h = $[D_h^{\,\prime}(s)/b]$ . The equations (10) become

$[\eqalignno{ D^{\prime}_{0}(s) & = iu_{0}D_{0}(s)+iu_{-h}D_{h}(s), & (14a) \cr D^{\prime}_{h}(s) & = ib\left(u_{0}+\alpha^{\prime}\right)D_{h}(s)+ibu_{h}D_{0}(s). & (14b) }]$

It is convenient to use the functions B_0,h(s) by setting

$[\eqalignno{ D_{0,h}(s) & = \exp\big[is\,{{u_{0}+b(u_{0}+\alpha^{\prime})}\over{2}}\big] B_{0,h}(s) \cr& = \exp\big[is(u_{0}+\omega)\big]B_{0,h}(s), & (15) }]$

with

$[\omega = {{b(u_{0}+\alpha^{\prime})-u_{0}}\over{2}}. \eqno(16)]$

Equations (14) become

$[\eqalignno{ B^{\prime}_{0}(s) & = -i\omega B_{0}(s)+iu_{-h}B_{h}(s), & (17a) \cr B^{\prime}_{h}(s) & = i\omega B_{h}(s)+ibu_{h}B_{0}(s). & (17b) }]$

They have special solutions of the form³ B₀(s) = exp(ias) and B_h(s) = $[\xi\exp(ias)]$ , which, introduced in equation (17), give ξ = bu_h/(a − ω) = (a + ω)/u_−h and

$[a^{2} = bu_{h}u_{-h}+\omega^{2}. \eqno(18)]$

The general solution of equation (17) is

$[\eqalignno{ B_{0}(s) & = c_{1}\exp(ias)+c_{2}\exp(-ias), & (19a) \cr B_{h}(s) & = c_{1}{{a+\omega}\over{u_{-h}}}\exp(ias) + c_{2}{{\omega-a} \over {u_{-h}}}\exp(-ias). & (19b) }]$

From the case s = 0, we obtain

$[\eqalign{ c_{1} & = B_{0}(0){{a-\omega} \over {2a}}+B_{h}(0){{u_{-h}} \over {2a}}, \cr c_{2} & = B_{0}(0) {{a+\omega}\over {2a}}-B_{h}(0){{u_{-h}} \over {2a}}.} \eqno(20)]$

Consequently

$[\eqalign{ B_{0}(s) = {}& B_{0}(0) {{(a-\omega)\exp(ias)-(a+\omega)\exp(-ias)} \over {2a}} \cr& + B_{h}(0)u_{-h}{{\exp(ias)-\exp(-ias)} \over {2a}},}]$

$[\eqalign{ B_{h}(s) = {}& B_{0}(0)bu_{h}{{\exp(ias)-\exp(-ias)} \over {2a}} \cr& + B_{h}(0){{(a+\omega)\exp(ias)-(\omega-a)\exp(-ias)} \over {2a}},}]$

or, in terms of D_0,h(s) [equation (15)] in a more compact form,

$[\eqalignno{ \exp\big[-is(u_{0}+\omega)\big]D_{0}(s) = {}& \left[\cos(as)-i{{\omega}\over{a}}\sin(as)\right]D_{0}(0) \cr& + i{{u_{-h}} \over {a}}\sin(as)D_{h}(0), & (21a) }]$

$[\eqalignno{ \exp\big[-is(u_{0}+\omega)\big]D_{h}(s) = {}& ib{{u_{h}} \over {a}}\sin(as)D_{0}(0) & (21b) \cr& +\left[\cos(as)+i{{\omega} \over {a}}\sin(as)\right]D_{h}(0). }]$

3.1. Expressions of α, ω and $[{\bi{a}}]$ as functions of the angles

Note that h = $[2k\sin\theta_{\rm{B}}]$ , θ_B being the geometrical Bragg angle, and $[{\bf{k}}_{0}.{\bf{h}}]$ = $[-kh\sin\theta]$ , with θ the glancing angle of $[{\bf{k}}_{0}]$ on the reflecting planes. From equation (3) we obtain

$[\alpha = 4\sin\theta_{\rm{B}}(\sin\theta-\sin\theta_{\rm{B}}) \simeq 2(\theta-\theta_{\rm{B}})\sin(2\theta_{\rm{B}}). \eqno(22)]$

Our definition of α [equation (3)] was made in such a way that α increases when θ increases.⁴ The approximated value of α is not valid far from the Bragg position or when θ_B approaches π/2 (normal incidence); therefore equation (3) is used in the crystalpy software.

α′ [equation (11)] and ω [equation (16)] are

$[\alpha^{\prime} = 2k\sin\theta_{\rm{B}}(\sin\theta-\sin\theta_{\rm{B}}) = h(\sin\theta-\sin\theta_{\rm{B}}), \eqno(23)]$

$[\omega = {{bh} \over {2}} \left(\sin\theta-\sin\theta_{\rm{B}}\right)+{{b-1}\over{2}}u_{0}. \eqno(24)]$

The `corrected Bragg angle' θ_c, that differs from θ_B because of the effect of refraction, is obtained as the θ value such that $[{\rm{Re}}\,\omega]$ = 0, or

$[\sin\theta_{\rm{c}}-\sin\theta_{\rm{B}} = {{1-b} \over {bh}}\,{\rm{Re}}\,(u_{0}), \eqno(25)]$

which, under the usual conditions [ $[\sin\theta_{\rm{c}}-\sin\theta_{\rm{B}}]$ ≃ $[(\theta_{\rm{c}}-\theta_{\rm{B}})\cos\theta_{\rm{B}}]$ ], reduces to

$[\theta_{\rm{c}} \simeq \theta_{\rm{B}}+{{1-b}\over{2b\sin(2\theta_{\rm{B}})}} \,{\rm{Re}}\,(\chi_{0}). \eqno(26)]$

From equations (24) and (25), the value of ω has a simple expression as a function of θ_c and θ,

$[\omega = {{bh} \over {2}}(\sin\theta-\sin\theta_{\rm{c}})+i{{b-1} \over {2}}{\rm{Im}}\, u_{0}. \eqno(27)]$

Note that, in our representation [using waves of the form $[\exp(i\,{\bf{k}}.{\bf{r}})]$ ], we have $[{\rm{Im}}\,u_{0}]$ ≥ 0. Equations (21) are expressed in terms of a, but they depend only on a². Using equation (27) in equation (18) we obtain

$[{\rm{Re}}\,a^{2} = \left[bh{{\sin\theta-\sin\theta_{\rm{c}}} \over {2 }}\right]^{2}-\,\left[{{b-1} \over {2}}{\rm{Im}}\,u_{0}\right]^{2}+\, {\rm{Re}}\,(bu_{h}u_{-h}), \eqno(28a)]$

$[{\rm{Im}}\,a^{2} = {{b(b-1)} \over {2}}\,h\,(\sin\theta-\sin\theta_{\rm{c}}) \, {\rm{Im}}\,u_{0}+{\rm{Im}}\,(bu_{h}u_{-h}). \eqno(28b)]$

Note that a can be expressed in terms of a² as

$[a = \pm\left[{\rm{sgn}}\left({\rm{Im}}\,a^{2}\right) \left({{{|a^{2}|+{\rm{Re}}\,a^{2}}\over{2}}}\right)^{1/2} + \ i\,\left({{{|a^{2}|-{\rm{Re}}\,a^{2}}\over{2}}}\right)^{1/2}\right]. \eqno(29)]$

3.2. The transfer matrix of a parallel crystal slab

The front and back surfaces of a crystal parallel slab correspond to s = 0 and s = t_c/γ₀ = T, respectively, with t_c the `usual' thickness of the crystal. We can express the fields at the back surface (D₀(T), D_h(T)) in terms of those at the front surface (D₀(0), D_h(0)) in matrix form,

$[\left( \matrix{ D_{0}(T) \cr D_{h}(T) } \right) = M \left( \matrix{ D_{0}(0) \cr D_{h}(0) } \right) = \left( \matrix{ m_{11} & m_{12} \cr m_{21} & m_{22} } \right) \left( \matrix{ D_{0}(0) \cr D_{h}(0) } \right). \eqno(30)]$

According to equation (21), the elements of the `transfer matrix' M are

$[m_{11} = \left[\cos(aT)-i\,{{\omega} \over {a}}\sin(aT)\right]\exp[iT(u_{0}+\omega)], \eqno(31a)]$

$[m_{12} = i\,{{u_{-h}}\over{a}}\sin(aT)\exp[iT(u_{0}+\omega)], \eqno(31b)]$

$[m_{21} = ib\,{{u_{h}}\over{a}}\sin(aT)\exp[iT(u_{0}+\omega)], \eqno(31c)]$

$[m_{22} = \left[\cos(aT)+i\,{{\omega} \over {a}}\sin(aT)\right]\exp[iT(u_{0}+\omega)]. \eqno(31d)]$

The determinant of the matrix M is $[{\rm{det}}(M)]$ = $[\exp[2iT(u_{0}+\omega)]]$ . Its modulus |det(M)| ≤ 1. It is 1 for a non-absorbing crystal (u₀ and ω are real). This is in agreement with the expected energy conservation. It can be verified that M(T₁ + T₂) = M(T₂)M(T₁) and M(−T) = [M(T)]⁻¹. Last, but not least, equation (30) is valid for both Bragg and Laue cases (with adequate values of b, a and ω).

3.3. The transfer matrix for the case of a `thick crystal'

Equations (21) and (31) are expressed in terms of a, but they depend only on a². It is possible to write them as

$[\eqalignno{ M = {}& {{\exp\big[i\,(u_{0}+\omega+a)T\big]}\over{2a}} \left( \matrix{ a-\omega & u_{-h} \cr bu_{h} & a+\omega } \right) \cr& + {{\exp\big[i\,(u_{0}+\omega-a)T\big]}\over{2a}} \left( \matrix{ a+\omega & -u_{-h} \cr -bu_{h} & a-\omega } \right), & (32) }]$

where the two terms are interchanged when a is changed in −a. They correspond to the two roots of a². They also correspond to the two branches of the dispersion surface [see, for example, Authier (2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction. Oxford University Press.]$ )]. The real part of the argument of the exponential factors, $[-T[{\rm{Im}}\,u_{0}(b+1)/2\pm{\rm{Im}}\,a]]$ , is related to the absorption. When $[{\rm{Im}}\,a\,\lt\,0]$ , the absorption is less than $[\exp[-T\,{\rm{Im}}\,u_{0}(b+1)/2]]$ for the first term, and more than that for the second one. Similarly, for $[{\rm{Im}}\,a\,\gt\,0]$ the two matrices present opposite behaviour. If $[T|{\rm{Im}}\,a|]$ is large (for example $[T|{\rm{Im}}\,a|\,{\gtrsim}\,5]$ ), we can keep only the largest term in equation (32),

$[M^{\,\rm{thick}} \simeq \left\{ \matrix{ {{\textstyle{\exp[i(u_{0}+\omega+a)T]}}\over{\textstyle{2a}}} \left( \matrix{ a-\omega & u_{-h} \cr bu_{h} & a+\omega } \right),_{\vphantom{\Big|}} \hfill \cr \qquad\qquad\qquad\qquad\qquad {\rm{with\ the\ choice}} \ {\rm{Im}}\, a\,\lt\,0,_{\vphantom{\big|}} \cr {{\textstyle{\exp[i(u_{0}+\omega-a)T]}}\over{\textstyle{2a}}} \left( \matrix{ a+\omega & -u_{-h} \cr -bu_{h} & a-\omega } \right),_{\vphantom{\Big|}} \hfill \cr \qquad\qquad\qquad\qquad\qquad {\rm{with\ the\ choice}} \ {\rm{Im}}\, a\,\gt\,0. } \right. \eqno(33)]$

A beam, that would be absorbed without Bragg diffraction if $[T\,{\rm{Im}}\,u_{0}]$ >> 1 may partially go through the `thick crystal' in the condition of Bragg diffraction. Equation (33) is a clear expression of the Borrmann effect.

3.4. Reflection and transmission amplitudes in the Laue case

In this case, b > 0. The boundary conditions are (D₀(0), D_h(0)) = (1, 0). The reflection and transmission amplitudes, r_L = D_h(T) and t_L = D₀(T), respectively, are directly written from the matrix equation (30),

$[r_{\,\rm{L}} = m_{21} = ibu_{h}{{\sin(aT)}\over{a}} \exp[iT(u_{0}+\omega)], \eqno(34a)]$

$[t_{\,\rm{L}} = m_{11} = \left[\cos(aT)-i\omega{{\sin(aT)}\over{a}}\right] \exp[iT(u_{0}+\omega)]. \eqno(34b)]$

The reflecting power is $[{\cal R}]$ = $[|r_{\,\rm{L}}|^{2}P]$ , where P = |b|⁻¹ is the ratio of the cross sections of the reflected and incident beams (Zachariasen, 1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ). Its plot $[{\cal R}(\theta)]$ as a function of the angle of incidence is the diffraction profile, often referred to as the rocking curve. $[{\cal R}(\theta)]$ and $[{\cal T(\theta)}]$ = $[|t_{\rm{L}}(\theta)|^{2}]$ are hereafter called reflectance and transmittance, respectively. An example is shown in Fig. 1.

Figure 1
Calculated reflectance (a) and transmittance (b) for a 10 µm-thick Laue Si 111 crystal at 8 keV, with 65° of asymmetric angle. The Bragg angle is θ_B = 14.31°.

It is also interesting to consider the case of incidence along the direction of $[{\bf{k}}_{h}]$ (diffraction vector $[-{\bf{h}}]$ ), for which (D₀(0), D_h(0)) = (0, 1). It is directly seen from equation (30) that the transmission and reflection amplitudes are $[\bar{t}_{\,\rm{L}}]$ = D_h(T) = m₂₂ and $[\bar{r}_{\,\rm{L}}]$ = D₀(T) = m₁₂ (note that the reflection power factor is P = |b| in this case). These results can be written as

$[\left( \matrix{ t_{\,\rm{L}} & \bar{r}_{\,\rm{L}} \cr r_{\,\rm{L}} & \bar{t}_{\,\rm{L}} } \right) = \left( \matrix{ m_{11} & m_{12} \cr m_{21} & m_{22} } \right) = M. \eqno(35)]$

This means that, for the Laue case, the matrix M can be considered not only as the `transfer-matrix' of the crystal slab but also as the `scattering-matrix' (S-matrix) which relates the vacuum waves leaving the crystal to the vacuum waves entering it, in analogy with the S-matrix used in general scattering theory.

The exponential factor in equation (34) gives a damping factor, which is {using u₀ + ω = [(b + 1)u₀ + bα′]/2 from equation (16), and noting that α′ is real}

$[\exp\big\{{\rm{Re}}\,[iT(u_{0}+\omega)]\big\} = \exp\left[-T{{b+1}\over{2}}\,{\rm{Im}}\,(u_{0})\right]. \eqno(36)]$

The Pendellösung effect is due to the oscillations of $[|\sin(aT)\,|^{2}]$ = $[\sin^{2}({\rm{Re}}\,aT)+\sinh^{2}({\rm{Im}}\,aT)]$ . The Pendellösung distance (depending on θ) along s₀ is thus equal to $[\pi/|{\rm{Re}}\,a|]$ = $[\lambda/|{\rm{Re}}\,\left({b\chi_{h}\chi_{-h}+w^{2}}\right)^{1/2}|]$ , where w = (λ/π)ω. At θ = θ_c, $[b\chi_{h}\chi_{-h}+w^{2}]$ = $[b\chi_{h}\chi_{-h}-[({\rm{Im}}\,\chi_{0}(b-1)/2]^{2}]$ . In the symmetric Laue case (b = 1) we obtain the well known formula of the Pendellösung distance along the direction normal to the crystal surface [see, for example, equation (3.48) of Pinsker (1978)],

$[\Lambda = {{ \lambda\cos\theta_{\rm{B}} }\over{ \left|{\rm{Re}}\,\left({\chi_{h}\chi_{-h}}\right)^{1/2}\right| }}. \eqno(37)]$

3.5. Reflection and transmission amplitudes in the Bragg case – the S-matrix

In this case, b < 0. We set D₀(0) = 1 and D_h(T) = 0 (which means no beam entering the crystal slab on the back surface). The reflection and transmission amplitudes are r_B = D_h(0) and t_B = D₀(T), respectively. Equation (30) is

$[\left( \matrix{ t_{\rm{B}} \cr 0 } \right) = \left( \matrix{ m_{11} & m_{12} \cr m_{21} & m_{22} } \right) \left( \matrix{ 1 \cr r_{\rm{B}} } \right), \eqno(38)]$

from which we obtain

$[r_{\rm{B}} = -{{m_{21}}\over{m_{22}}} = {{-ibu_{h}\sin(aT)}\over{a\cos(aT)+i\omega\sin(aT)}}, \eqno(39a)]$

$[t_{\rm{B}} = m_{11}+m_{12}r_{\rm{B}} = {{a\exp[iT(u_{0}+\omega)]}\over{a\cos(aT)+i\omega\sin(aT)}}. \eqno(39b)]$

These solutions, as well as those for Laue in equations (34), can also be obtained by direct integration of the TT equations (17) using Laplace transforms (see Appendix B). Similarly, in the case of incidence on the crystal back side along the direction $[{\bf{k}}_{h}]$ (diffraction vector $[-{\bf{h}}]$ ), we set D_h(T) = 1, D₀(T) = $[\bar{r}_{\rm{B}}]$ , D₀(0) = 0 and D_h(0) = $[\bar{t}_{\rm{B}}]$ . Therefore, equation (30) gives

$[\left( \matrix{ \bar{r}_{\rm{B}} \cr 0 } \right) = \left( \matrix{ m_{11} & m_{12} \cr m_{21} & m_{22} } \right) \left( \matrix{ 0 \cr \bar{t}_{\rm{B}} } \right), \eqno(40)]$

from which we obtain

$[\bar{t}_{\rm{B}} = {{1}\over{m_{22}}} = {{a\exp[-iT(u_{0}+\omega)]}\over{a\cos(aT)+i\omega\sin(aT)}}, \eqno(41a)]$

$[\bar{r}_{\rm{B}} = m_{12}\bar{t}_{\rm{B}} = {{iu_{-h}\sin(aT)}\over{a\cos(aT)+i\omega\sin(aT)}}. \eqno(41b)]$

Consequently, the S-matrix for the Bragg case, defined as

$[\left( \matrix{ D_{0}(T) \cr D_{h}(0) } \right) = S\left( \matrix{ D_{0}(0) \cr D_{h}(T) } \right), \eqno(42)]$

$[S = \left( \matrix{ t_{\rm{B}} & \bar{r}_{\rm{B}} \cr r_{\rm{B}} & \bar{t}_{\rm{B}} } \right) = \left( \matrix{ m_{11}-{{m_{12}m_{21}}\over{m_{22}}} & {{m_{12}}\over{m_{22}}} \cr -{{m_{21}}\over{m_{22}}} & {{1}\over{m_{22}}} } \right). \eqno(43)]$

The diffraction profile (reflectance) is $[{\cal R}(\theta)]$ = $[|r_{\rm{B}}(\theta)|^{2}P]$ , with P = 1/|b|, and the transmittance is $[{\cal T}(\theta)]$ = $[|t_{\rm{B}}(\theta)|^{2}]$ . An example is shown in Fig. 2.

Figure 2
Calculated reflectance (a) and transmittance (b) of a symmetrical Bragg Si 111 crystal with 10 µm thickness at 8 keV.

The field inside the crystal, i.e. D₀(s) and D_h(s), can be calculated using equation (19), with D₀(0) = 1 and D_h(0) = r_B from equation (39a) we obtain

$[\eqalignno{ D_{h}(s) & = {{ibu_{h}\sin(as-aT)}\over{a\cos(aT)+i\omega\sin(aT)}} \exp[is(u_{0}+\omega)] \cr& = r_{\rm{B}}{{\sin(aT-as)}\over{\sin(aT)}} \exp[is(u_{0}+\omega)], & (44a) }]$

$[D_{0}(s) = {{a\cos(aT-as)+i\omega\sin(aT-as)}\over{a\cos(aT)+i\omega\sin(aT)}} \exp[is(u_{0}+\omega)]. \eqno(44b)]$

An example of simulation of the field inside the crystal using equation (44) is shown in Fig. 3. For the Laue case, also shown in this figure, we observe that the field at coordinate s is simply calculated by the equations (34) replacing T by s.

Figure 3
Calculations for a symmetric Si 111 at 8 keV with thickness t_c = 50 µm. The graphs show the electric field intensity inside the crystal as a function of the deviation angle θ − θ_B and penetration ratio −s/T (equivalent to a depth ratio − t/t_c), for (a) Bragg |D_h|², (b) Bragg |D₀|², (c) Laue |D_h|², (d) Laue |D₀|².

Fig. 3(b) shows that the penetration of the incident wave inside the crystal is small in a limited interval around θ_c. In Fig. 4 we fitted the intensity profile of |D₀(s)|² versus depth for each value of θ − θ_B. The fact that for a thick crystal in Bragg geometry |D₀(s)|² has significant values only in the vicinity of the crystal surface, in the central region, can be explained from equation (44). The moduli of the functions $[\sin(aT-as)]$ and $[\cos(aT-as)]$ are approximately proportional to $[\exp[(T-s)|{\rm{Im}}\,a|]]$ if the argument of this exponential function is sufficiently large. Consequently, |D₀(s)|² is nearly proportional to $[\exp\{-s[2|{\rm{Im}}\,a|+(b+1)\,{\rm{Im}}\,u_{0}]\}]$ . Writing $[|D_{0}(s)|^{2}]$ = $[\exp[-s/s_{\rm{ext}}]]$ , with $[s_{\rm{ext}}]$ = $[[(2|{\rm{Im}}\,a|+(b+1)\,{\rm{Im}}\,u_{0})]^{-1}]$ the extinction length (measured along the s₀ axis), we obtain for the Bragg symmetric (b = −1) case,

$[s_{\rm{ext}} = {{1}\over{2|{\rm{Im}}\,a|}}. \eqno(45)]$

Figure 4
Fit of |D₀|² versus t of Fig. 3

(b) with a function $[\exp(-t/t_{\rm{ext}})]$ to obtain the extinction depth t_ext. This result is compared with the calculated value from equation (45)

that gives $[t_{\rm{ext}}]$ = $[\sin\theta_{\rm{B}}s_{\rm{ext}}]$ = $[\sin\theta_{\rm{B}}/(2|{\rm{Im}}\,a|)]$ .

3.5.1. Reflection amplitude for a thick absorbing crystal in the Bragg case

In the case of thick (or semi-infinite) Bragg crystal, the reflection amplitude, given by equation (39), takes the form [using M from equation (33)]

$[r_{\rm{B}}^{\rm{thick}} = -{{m_{21}}\over{m_{22}}} = \left\{ \matrix{ -{{bu_{h}}\over{a+\omega}}={{\omega-a}\over{u_{-h}}}, & \!{\rm{with\ the\ choice\ Im}}\,a\,\lt\,0, \cr {{bu_{h}}\over{a-\omega}}={{\omega+a}\over{u_{-h}}}, & \!{\rm{with\ the\ choice\ Im}}\,a\,\gt\,0. } \right. \eqno(46)]$

Both equations are equivalent assuming that the sign in a = ±(a²)^1/2 is correctly selected. The condition on $[{\rm{Im}}\,(a)]$ in equation (46) is in accordance with the physical condition that |r_B| goes to zero when $[|\sin\theta-\sin\theta_{\rm{c}}|]$ is large. The condition $[{\rm{Im}}\,a\,\lt\,0]$ [ $[{\rm{Im}}\,a\,\gt\,0]$ ] is equivalent⁵ to $[{\rm{sgn}}({\rm{Re}}\,a)]$ = $[{\rm{sgn}}({\rm{Re}}\,\omega)]$ [ $[{\rm{sgn}}({\rm{Re}}\,a)]$ = $[-{\rm{sgn}}({\rm{Re}}\,\omega)]$ ] for large values of $[|\sin\theta-\sin\theta_{\rm{c}}|]$ .

Equation (46) is a useful result, as most crystal monochromators used in synchrotron radiation are thick crystals in Bragg (reflection) mode.

3.5.2. Reflection amplitude for non-absorbing crystals in the Bragg case

In this case $[u_{-h}]$ = u^*_h; ω [see equation (16)] and a² = $[\omega^{2}-|b|u_{h}u_{h}^{*}]$ are real. We can distinguish two cases.

If a² ≤ 0, or $[\omega^{2}]$ $[\leq]$ $[|b|u_{h}u_{h}^{*}]$ , then a = $[i\left({|b|u_{u}u_{h}^{*}-\omega^{2}}\right)^{1/2}]$ ; therefore, according to equation (46),

$[r_{\rm{B}}^{\rm{thick,transparent}} = {{1}\over{u_{h}^{*}}} \left[ \omega+i\left({|b|u_{h}u_{h}^{*}-\omega^{2}}\right)^{1/2} \right]. \eqno(47)]$

If a² > 0, or $[\omega^{2}]$ > $[|b|u_{h}u_{h}^{*}]$ ,

$[r_{\rm{B}}^{\rm{thick,transparent}} = {{1}\over{u_{h}^{*}}} \left[\omega-{\rm{sgn}}(\omega) \left({\omega^{2}-|b|u_{h}u_{h}^{*}}\right)^{1/2} \right]. \eqno(48)]$

Equation (48) represents the tails of the reflection profile. As discussed previously, the sign selection is such that |r_B| tends to zero for large values of |ω|. Equation (47) corresponds to the zone of total reflection. The reflectance is

$[{\cal R}_{\rm{B}}^{\rm{transparent,thick}} = \left\{ \matrix{ 1,\hfill & {\rm{if}}\ |y|\,\leq\,1,\hfill \cr \big[\,y-\left({y^{2}-1}\right)^{1/2}\big]^{2},\hfill & {\rm{if}}\ |y|\,\gt\,1,\hfill } \right. \eqno(49)]$

with y = $[\omega/\left({|b|u_{h}u_{h}^{*}}\right)^{1/2}.]$

3.6. Calculation of reflection and transmission amplitudes using the transfer matrix

The matrix method permits the complex reflection and transmission amplitudes of a crystal made by layers of different crystals (or the same crystal with different orientations) to be obtained. For that, (i) construct the transfer matrix of the total crystal by multiplication⁶ of the transfer matrices of the different layers [each one calculated using equation (31)]; (ii) if the geometry is Laue, obtain reflection and transmission amplitudes using the coefficients m₂₁ and m₁₁, respectively, of this matrix [equation (34)]; otherwise (Bragg geometry), compute the scattering matrix using equation (43) and the reflection and transmission amplitudes are given in the matrix terms s₂₁ and s₁₁ [equation (39)], respectively.

A first example shows how simple is the application of this recipe of multiplication of transfer matrices to obtain the reflectance of a simple two-layer crystal. Consider a bilayer of two identical crystal layers of thickness T and transfer matrix M for each one. Using matrix analysis, the transfer matrix of the bilayer is [M(T)]² = M(2T) from which it is easy to compute the reflectivity in Bragg geometry [equation (39)]. Otherwise, if this result would be obtained via the reflectivities (r and $[\bar{r}]$ ) and transmittivities (t and $[\bar{t}]$ ) of the single layer (S-matrix), the reflectivity of the bilayer results from an infinite series as shown in Fig. 5.

Figure 5
Example of calculation of the reflection amplitude r₂ of a Si 111 crystal of 2 µm thickness from the amplitudes of the half-layer (1 µm). The reflectivity of the bilayer r₂ can be obtained as an infinite sum r₂ = $[r+rt\bar{t}]$ + $[rt\bar{t}(r\bar{r})]$ + $[rt\bar{t}(r\bar{r})^{2}]$ + … = $[r[1+t\bar{t}\sum_{n\,=\,0}^{\infty}(r\bar{r})^{n}]]$ = $[r\{1+[t\bar{t}/(1-r\bar{r})]\}]$ . Calculations done with crystalpy for a photon energy of 8 keV.

A second example is the Bragg reflection of a crystal layer on a thick substrate. The transfer matrix is calculated as

$[M^{\,\prime} = M^{\rm{thick}}\times M, \eqno(50)]$

with M^thick the transfer matrix of the substrate and M the transfer matrix of the thin layer. We are interested in the Bragg reflectivity or

$[r_{\rm{B}} = -{{m^{\prime}_{21}}\over{m^{\prime}_{22}}} = -{{m_{21}^{\rm{thick}}m_{ 11}+m_{22}^{\rm{thick}}m_{21}} \over {m_{21}^{\rm{thick}}m_{12}+m_{22}^{\rm{thick}}m_{22}}}. \eqno(51)]$

This can be expressed as a function of the substrate reflectivity $[r_{\rm{S}}]$ = $[-m_{21}^{\rm{thick}}/m_{22}^{\rm{thick}}]$ giving

$[r_{\rm{B}} = {{r_{\rm{S}}m_{11}-m_{12}}\over{m_{22}-r_{\rm{S}}m_{11}}}. \eqno(52)]$

The method of transfer matrix multiplication can also be used for analysing distorted and bent crystals and will be explored in a future work.

3.7. The direction of the diffracted wave in vacuum

In some applications, as in ray tracing, it is essential to know the diffracted wavevector $[{\bf{k}}_{\rm{out}}]$ that exits from the crystal. As mentioned before, the choice of $[{\bf{k}}_{0}]$ in equation (4) is somewhat arbitrary. In our choice, $[{\bf{k}}_{0}]$ corresponds exactly to the wavevector of the incident plane wave. The vector $[{\bf{k}}_{h}]$ (see Section 2), defined as $[{\bf{k}}_{h}]$ = $[{\bf{k}}_{0}+{\bf{h}}]$ , does not correspond in general to the wavevector of the outgoing ray or wave outside the crystal. $[{\bf{k}}_{\rm {out}}]$ has the form

$[{\bf{k}}_{\rm{out}} = {\bf{k}}_{h}+\beta{\bf{n}},]$

with n the unit vector along the inward normal to the crystal exit surface. The (real) coefficient β is obtained by writing that the modulus of $[{\bf{k}}_{\rm{out}}]$ is equal to k,

$[|{\bf{k}}_{\rm{out}}|^{2} = |{\bf{k}}_{h}+\beta{\bf{n}}|^{2} = k^{2}.]$

Note that $[|{\bf{k}}_{h}|^{2}]$ = $[k^{2}(1-\alpha)]$ and $[{\bf{k}}_{h}.{\bf{n}}]$ = $[\gamma_{h}k\left({1-\alpha} \right)^{1/2}]$ , from which we obtain the equation

$[k^{2} = k^{2}(1-\alpha)+2\beta k\gamma_{h}\left({1-\alpha}\right)^{1/2}+\beta^{2}.]$

Its solutions are

$[\beta = -k\gamma_{h}\left({1-\alpha}\right)^{1/2} \,\pm\, k\left[{\alpha+\gamma_{h}^{2}(1-\alpha)}\right]^{1/2},]$

where the ± sign is chosen in such a way that β = 0 when α = 0 (i.e. positive for γ_h > 0 and negative for γ_h < 0).

4. The crystalpy library

Crystalpy is a Python library that performs calculations on diffraction from perfect crystals using the formalism introduced in the previous chapters.

The motivation of crystalpy was to create a modern, extensible, well documented and friendly library to overcome the difficulties of integrating ancient software tools based on the dynamical diffraction theory. It is specifically designed for two objectives: support for new versions of the crystal diffraction codes delivered in platforms like OASYS (Rebuffi & Sanchez del Rio, 2017), and to provide a core for ray-tracing simulations with crystals. The crystalpy library is written in the Python language and uses standard libraries (NumPy and SciPy). It makes use of vector calculus and stack operations to accelerate the calculations. Therefore, it is adapted for being used in ray-tracing tools, such as the future SHADOW (Sanchez del Rio et al., 2011) versions.

To simulate a diffraction experiment using a perfect crystal, crystalpy offers functions that implement the theory described previously. Two input objects must be prepared: (i) the incident wave(s) or photon ray(s), and (ii) the information on the crystal (diffraction setup). The objects representing these two entities are described here.

The Photon class is a minimum class for a photon, containing the energy (in eV) and a unit direction vector, implemented in the Vector class. It deals with the storage and operations (addition, scalar product, cross product, normalization, rotation around an axis, etc.) of a 3D vector. A superclass of Photon is ComplexAmplitudePhoton, that contains the scalar complex amplitude for σ and π polarizations). These classes ( Vector, Photon and ComplexAmplitudePhoton) can hold stacks (the internal storage is done with arrays to speed-up vector operations). The ComplexAmplitudePhoton class has a corresponding ComplexAmplitudePhotonBunch superclass, decorated with methods to deal with multiple waves or beams (bunches or sets of photons).

The information on the crystal itself (e.g. particular crystal material and crystal structure), its preparation (crystal cut) and related physical parameters (like the structure factor) are managed by the DiffractionSetup classes. crystalpy allows multiple options to retrieve the crystal structure and the scattering functions needed to calculate the structure factors. The DiffractionSetupAbstract class defines the methods to access the basic information of the crystal (defined as a string, e.g. `Si') such as angleBragg, dSpacing and unitCellVolume, and to compute the structure factors: F0, FH, FH_ bar. These parameters can be obtained from several libraries external to crystalpy. We implemented three options: (i) DiffractionSetupXraylib using the $[{\it{xraylib}}]$ library (Schoonjans et al., 2011 ), (ii) DiffractionSetupDabax that uses the dabax library (Sanchez del Rio, 2021 ), and (iii) using an ad hoc generated data file. This modular structure permits disconnecting the calculation part from the access to optical and physical constants. Indeed, when using ad hoc data files we do not have to import $[{\it{xraylib}}]$ or dabax packages. We implemented this for the crystal material files of the SHADOW (Sanchez del Rio et al., 2011) code in the traditional version ( DiffractionSetupShadowPreprocessorV1), and in a version supporting d-spacing crystals ( DiffractionSetupShadowPreprocessorV2). The DiffractionSetup classes handle the information about the crystal setup and collect all the parameters needed to fully define the physical system we are modelling: geometry_ type (among BraggDiffraction, BraggTransmission, LaueDiffraction and LaueTransmission), crystal_ name (a string, e.g. Si, Ge), thickness (crystal thickness in SI units [m]), miller _ h, miller _ k, miller_ l (the Miller indices) and asymmetry _ angle [angle in degrees between the crystal surface and the planes hkl as definedby Sanchez del Rio et al. (2015)].

The determination of crystal structure factors (necessary to compute χ_h, χ_−h and thus u_h and u_−h) is not trivial, and requires the list of the crystallographic parameters (basically the cell parameters and a list of the atoms of the unit cell, with their occupation and coordinates). Both $[{\it{dabax}}]$ and $[{\it{xraylib}}]$ libraries use similar methods that are detailed by Yu et al. (2022 ). This implementation allows any possible crystal structure. Complex crystals such as alpha-quartz (Sutter et al., 2022 ; Sutter et al., 2023 ), or YB₆₆ (Yu et al., 2022) are considered. However, some particularities regarding chirality, strong anisotropy or temperature dependence may not be included accurately and are sometimes modelled by phenomenological parameters.

To perform the main calculations (reflectivities, transfer matrices, diffracted photons, etc.) several methods in the Diffraction class are used, getting the crystal setup and the photon bunch as inputs. For the moment, only flat perfect crystals are coded (in the PerfectCrystalDiffraction class) which directly implements the formulation and theory in Section 3. For completeness, crystalpy also includes the equations of Zachariasen (Zachariasen, 1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ) and can be used instead of the formalism described in this paper. Typical angle or photon scans, as shown in Fig. 2, are calculated defining a ComplexAmplitudePhoton entity for each point, grouping them in a ComplexAmplitudePhotonBunch and then calculating the diffraction by the crystal using calculateDiffractedComplexAmplitudes.

A user-friendly application has been written in the OASYS environment to compute diffraction profiles using crystalpy (Fig. 6). The applications automatically generates a script that can be used for further batch calculations.

Figure 6
Interactive application for computing the perfect crystal diffraction profiles using crystalpy and available in OASYS.

In the ray-tracing SHADOW4, all calculations related to crystal optics are delegated to crystalpy. Ray-tracing permits simulations of beamline optics including a realistic description of the source. It also allows the simulation of curved crystals, under the assumption that the local reflectivity of the curved crystal is the same as for the flat crystal. This assumption is not always granted and has to be verified before performing ray-tracing simulations including curved crystals.

5. Summary and conclusions

We have presented a theoretical and numerical description of dynamical diffraction in perfect crystals. In the first part of this paper we presented a new perspective of the well known dynamical theory of diffraction applied to undeformed perfect crystals. We deduced the equations of diffraction amplitudes (as well as intensity reflectance and transmittance) starting from basic principles via the solution of the Takagi–Taupin equations. We calculated the transfer matrix, useful to compute the diffraction of stacked crystal layers, and also the scattering matrix, of interest for the Bragg case. For completeness, our results are compared with those presented in the well known textbook by Zachariasen (1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ) (see Appendix C). In the second part we presented crystalpy, a software library completely written in Python that implements the theory previously discussed. This open source tool can be used to predict the diffraction properties of any crystal structure, like Si, Ge or diamond typically used in synchrotron beamlines, but also for any other crystal provided its crystalline structure is known. This library is intended to replace multiple scattered pieces of software in packages like OASYS (Rebuffi & Sanchez del Rio, 2017) and is designed to be the kernel of the crystal calculations in version 4 of the SHADOW (Sanchez del Rio et al., 2011) ray-tracing code. The crystalpy library and its documentation are available from https://github.com/oasys-kit/crystalpy.

APPENDIX A

Derivation of the TT equations for a rotating perfect crystal

In the `rotating crystal mode', the crystal is rotated around an axis perpendicular to the `diffraction plane' which contains the diffraction vector h and the wavevector $[{\bf{k}}_{0}]$ of the fixed incident plane-wave in vacuum. The crystal rotation from the exact geometrical Bragg position may be viewed as a special kind of crystal deformation. We propose to use the Takagi–Taupin approach for the deformed crystal to derive the basic results of the dynamic theory for perfect crystal diffraction. The X-ray wavefield inside the crystal is set as

$[\eqalignno{ \Psi({\bf{r}}) & = \exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right) \left[A_{0}({\bf{r}})+\exp\left(i\,{\bf{h}}_{\rm{B}}.{\bf{r}}\right) A_{h}({\bf{r}})\right] \cr& = A_{0}({\bf{r}})\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right) + A_{h}({\bf{r}})\exp\left(i\,{\bf{k}}_{h{\rm{B}}}.{\bf{r}}\right). & (53) }]$

$[{\bf{h}}_{\rm{B}}]$ is the position of the diffraction vector when the crystal is in the exact Bragg condition for the fixed incident wavevector $[{\bf{k}}_{0}]$ . The vector $[{\bf{k}}_{h{\rm{B}}}]$ = $[{\bf{k}}_{0}]$ + $[{\bf{h}}_{\rm{B}}]$ is therefore such that $[|{\bf{k}}_{0}|]$ = $[|{\bf{k}}_{h{\rm{B}}}|]$ = k = $[2\pi/\lambda]$ . The Fourier coefficients χ_h of the perfect crystal susceptibility are replaced by the function $[\chi_{h}\exp[i\phi({\bf{r}})]]$ , in which $[\phi({\bf{r}})]$ = $[-{\bf{h}}_{\rm{B}}.{\bf{u}}({\bf{r}})]$ , with $[{\bf{u}}({\bf{r}})]$ the displacement field of the deformed crystal. In such conditions, the following form of the TT equations,

$[2i\,{\bf{k}}_{0}.\nabla A_{0}+\chi_{0}k^{2}A_{0}+\chi_{-h}k^{2} \exp(-i\phi)A_{h} = 0,]$

$[2i\,{\bf{k}}_{h{\rm{B}}}.\nabla A_{h}+\chi_{0}k^{2}A_{h}+\chi_{h}k^{2} \exp(+i\phi)A_{0} = 0,]$

is obtained by inserting equation (53) into (1), with the following approximations: the second-order derivatives of A_0,h, supposed to be slowly varying amplitudes, are neglected and only the terms containing $[\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right)]$ in the product χΨ are considered. Introducing oblique coordinates (s₀, s_h) in the diffraction plane, along the directions of $[{\bf{k}}_{0}]$ and $[{\bf{k}}_{h{\rm{B}}}]$ , so that $[{\bf{k}}_{0}.\nabla A_{0}]$ = $[k{{\partial A_{0}}/{\partial s_{0}}}]$ and $[{\bf{k}}_{h{\rm{B}}}.\nabla A_{h}]$ = $[k{{\partial A_{h}}/{\partial s_{h}}}]$ , the TT equations are

$[{{\partial A_{0}} \over {\partial s_{0}}} = {{ik} \over {2}}\left[\chi_{0}A _{0}+\chi_{-h}\exp(-i\phi)A_{h}\right],]$

$[{{\partial A_{h}} \over {\partial s_{h}}} = {{ik} \over {2}}\left[\chi_{0}A _{h}+\chi_{h}\exp(+i\phi)A_{0}\right].]$

Performing the transformation A₀ = D₀ and $[\exp(i\phi)A_{h}]$ = D_h, we obtain

$[{{\partial D_{0}} \over {\partial s_{0}}} = {{ik} \over {2}}\left[\chi_{0}D _{0}+\chi_{-h}D_{h}\right], \eqno(54a)]$

$[\displaystyle{{\partial D_{h}} \over {\partial s_{h}}} = {{ik} \over {2}}\left[\left(\chi_{0} +{{2} \over {k}}{{\partial\phi} \over {\partial s_{h}}}\right)D_{h}+\chi_{h}D_{0}\right]. \eqno(54b)]$

This is identical to equation (10) considering that $[\alpha]$ = $[({{2}/{k}})({{\partial\phi}/{\partial s_{h}}})]$ , for which a demonstration follows.

A1. Demonstration of α = (2/k)(∂ϕ/∂s_h)

Let $[{\bf{i}}_{0,h}]$ be unit vectors along the fixed directions of $[{\bf{k}}_{0}]$ and $[{\bf{k}}_{h{\rm{B}}}]$ ; the crystal rotation transforms $[{\bf{i}}_{0,h}]$ into $[{\bf{j}}_{0,h}]$ . A position vector $[s_{0}{\bf{i}}_{0}+s_{h}{\bf{i}}_{h}]$ is transformed into $[s_{0}\,{\bf{j}}_{0}+s_{h}\,{\bf{j}}_{h}]$ . The displacement field is $[{\bf{u}}(s_{0},s_{h})]$ = $[s_{0}(\,{\bf{j}}_{0}-{\bf{i}}_{0})]$ + $[s_{h}(\,{\bf{j}}_{h}-{\bf{i}}_{h})]$ ; $[{\bf{h}}_{\rm{B}}]$ = $[k({\bf{i}}_{h}-{\bf{i}}_{0})]$ . Hence

$[\phi(s_{0},s_{h}) = {\bf{h}}_{\rm{B}}.{\bf{u}} = ks_{0}({\bf{i}}_{h}-{\bf{i}}_{ 0}).(\,{\bf{j}}_{0}-{\bf{i}}_{0})+ks_{h}({\bf{i}}_{h}-{\bf{i}}_{0}).(\,{\bf{j}}_{h}-{\bf{i}}_{h}).]$

We note that

$[\eqalign{ {\bf{i}}_{0}.{\bf{j}}_{0} & = {\bf{i}}_{h}.{\bf{j}}_{h} = \cos \Delta\theta_{\rm{B}}, \cr \displaystyle{\bf{i}}_{0}.{\bf{j}}_{h} & = \cos(2\theta_{\rm{B}}+\Delta\theta), \cr \displaystyle{\bf{i}}_{h}.{\bf{j}}_{0} & = \cos(2\theta_{\rm{B}}-\Delta\theta), \cr \displaystyle{\bf{i}}_{0}.{\bf{i}}_{h} & = \cos 2\theta_{\rm{B}}, }]$

and therefore

$[\eqalign{ {{2} \over {k}}{{\partial\phi} \over {\partial s_{h}}} & = 2({\bf{i}}_{h}- {\bf{i}}_{0}).(\,{\bf{j}}_{h}-{\bf{i}}_{h}) \cr& = 2(\cos\Delta\theta-\cos(2\theta_{\rm{B}}+\Delta\theta)-1+\cos 2\theta _{\rm{B}}) \cr& = 2[(\cos\Delta\theta-1)(1-\cos 2\theta_{\rm{B}})+\sin 2\theta_{\rm{B}}\sin \Delta\theta] \cr& = 4\sin\theta_{\rm{B}}[\sin\theta_{\rm{B}}(\cos\Delta\theta-1)+\cos\theta_{B }\sin\Delta\theta] \cr& = 4\sin\theta_{\rm{B}}[\sin(\theta_{\rm{B}}+\Delta\theta)-\sin\theta_{\rm{B}}] = \alpha. }]$

APPENDIX B

Solutions of TT equations (17) using the Laplace transform

B1. Laue solution based on Laplace transform

Let us denote $[\bar{F}(p)]$ the Laplace transform of a function F(s)

$[\bar{F}(p) = \int\limits_{0}^{\infty} {\rm{d}}s\,\exp(-ps)\,F(s).]$

Applying the Laplace transform to equations (17) we get

$[\eqalign{ (\,p+i\omega)\bar{B}_{0}(\,p)-iu_{-h}\bar{B}_{h}(\,p) & = 1, \cr (\,p-i\omega)\bar{B}_{h}(\,p)-ibu_{h}\bar{B}_{0}(\,p) & = 0. }]$

The solutions are

$[\eqalign{ \bar{B}_{0}(\,p) & = {{(\,p-i\omega)} \over {p^{2}+a^{2}}}, \cr \bar{B}_{h}(\,p) & = {{ibu_{h}} \over {p^{2}+a^{2}}}, }]$

with, as previously defined, a² = ω² + bu_hu_−h, a = (ω² + bu_hu_−h)^1/2, hence one retrieves the same results of equations (34) using the fact that (p² + a²)⁻¹ and p(p² + a²)⁻¹ are the Laplace transforms of sin(as)/a and cos(as), respectively.

B2. Bragg solution based on Laplace transform

By Laplace transform of equation (17), and calling r′ = B_h(0), we obtain

$[\eqalign{ (\,p+i\omega)\bar{B}_{0}(\,p)-iu_{-h}\bar{B}_{h}(\,p) & = 1, \cr (\,p-i\omega)\bar{B}_{h}(\,p)-iu_{h}\bar{B}_{0}(\,p) & = r, }]$

$[\eqalign{ \bar{B}_{0}(\,p) & = {{p-i\omega+iru_{-h}} \over {p^{2}+a^{2}}}, \cr \bar{B}_{h}(\,p) & = {{r\,(\,p+i\omega)+ibu_{h}} \over {p^{2}+a^{2}}}, }]$

with (the same as before) a² = ω² + bu_hu_−h. Hence,

$[B_{0}(s) = \cos(as)+i\,(ru_{-h}-\omega){{\sin(as)}\over{a}},]$

$[B_{h}(s) = r\left[\cos(as)+i\omega{{\sin(as)}\over{a}}\right] + ibu_{h}{{\sin (as)} \over {a}}.]$

r and then the reflected and transmitted amplitudes are obtained using the condition D_h(T) = B_h(T) = 0. With some calculation, we obtain the results of equations (39).

APPENDIX C

Equivalence of amplitudes in equation (39) and (34) with Zachariasen's formulation

The diffracted and transmitted intensities (not the amplitudes) are derived in the book of Zachariasen (Zachariasen, 1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ) (first edition published in 1944). It is shown hereafter that the amplitudes can be easily derived from Zachariasen's formalism. For that purpose, we use Zachariasen's notation and equations.

In the Laue case, using equations [3.127] and [3.128] of Zachariasen (1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ), we obtain

$[t_{\rm{L}} = c_{1}D^{\,\prime}_{0}+c_{2}D^{\,\prime\prime}_{0} = {{c_{1}x _{2}-c_{2}x_{1}} \over {x_{2}-x_{1}}}, \eqno(55a)]$

$[r_{\rm{L}} = x_{1}c_{1}D^{\,\prime}_{0}+x_{2}c_{2}D^{\,\prime\prime}_{0} = x_ {1}x_{2}{{c_{1}-c_{2}} \over {x_{2}-x_{1}}}. \eqno(55b)]$

Similarly, in the Bragg case, using Zachariasen's equations [3.127] and [3.135], we obtain⁷

$[t_{\rm{B}} = c_{1}D^{\,\prime}_{0}+c_{2}D^{\,\prime\prime}_{0} = c_{1}c_{2} {{x_{2}-x_{1}} \over {c_{2}x_{2}-c_{1}x_{1}}}, \eqno(56a)]$

$[r_{\rm{B}} = x_{1}D^{\,\prime}_{0}+x_{2}D^{\,\prime\prime}_{0} = x_{1}x_{2} {{c_{2}-c_{1}} \over {c_{2}x_{2}-c_{1}x_{1}}}. \eqno(56b)]$

The symbols c and x are

$[\eqalignno{ c_{1} & = \exp(-2\pi ik_{0}\delta^{\,\prime}_{0}t/\gamma_{0}), \cr c_{2} & = \exp(-2\pi ik_{0}\delta^{\,\prime\prime}_{0}t/\gamma_{0}). &(57) }]$

γ₀ (γ_h) is the direction cosine of the incident (diffracted) wave and the other quantities are defined as

$[\left(\matrix{\delta_{0}^{\,\prime}\cr \delta_{0}^{\,\prime\prime}}\right) = {{1} \over {2}}\left(\Psi_{0}-z\pm X\right),]$

$[x_{1,2} = {{-z\pm X} \over {\Psi_{\bar{H}}}},]$

$[z = {{1-b} \over {2}}\Psi_{0}+{{b} \over {2}}\alpha_{Z},]$

$[\alpha_{Z} = {{1} \over {|{\bf{k}}_{0}|^{2}}}\left[|{\bf{B}}_{H}|^{2} + 2{\bf{k}}_{0}\cdot{\bf{B}}_{H}\right],]$

with X = (q + z²)^1/2, q = $[b\Psi_{H}\Psi_{\bar{H}}]$ , Ψ_H is the Fourier component of the electrical susceptibility Ψ₀ and b = γ₀/γ_h is the asymmetry factor.

It is easy to see that $[x_{2}-x_{1}]$ = $[2X/\Psi_{\bar{H}}]$ , x₁x₂ = $[-b\Psi_{H}/\Psi_{\bar{H}}]$ . Introducing the variables c = $[\exp[-2\pi ik_{0}(\Psi_{0}-z)t/(2\gamma_{0})]]$ and m = −2πk₀Xt/(2γ₀), we have

$[c_{1}-c_{2} = c\,\big[\exp(im)-\exp(-im)\big] = 2ic\sin(m),]$

and

$[\eqalign{ x_{2}c_{1}-x_{1}c_{2} & = {{c}\over{\Psi_{\bar{H}}}} \left[-(X+z)\exp(im)-(X-z)\exp(-im)\right] \cr & = {{2c} \over {\Psi_{\bar{H}}}}\big[-X\cos(m)-iz\sin(m)\big]. }]$

Replacing in equation (56) the terms obtained here we finally get

$[\eqalign{ r^{\,Z}_{\rm{L}} & = icb\Psi_{H}{{\sin(m)} \over {X}}, \cr t^{\,Z}_{\rm{L}} & = c\left[\cos(m)+i{{z} \over {X}}\sin(m)\right]. }]$

For the Bragg case, we pre-calculate

$[\eqalign{ x_{2}c_{2}-x_{1}c_{1} & = {{c} \over {\Psi_{\bar{H}}}}\big[\!-\!(X+z)\exp(-im)-(X-z)\exp(im)\big] \cr& = {{2c} \over {\Psi_{\bar{H}}}}\big[\!-\!X\cos(m)+iz\sin(m)\big], }]$

that we introduced in equation (56) and we finally get

$[\eqalignno{ r^{\,Z}_{\rm{B}} & = ib\Psi_{H}{{\sin(m)} \over {-X\cos(m)+iz\sin(m)}}, & (58a) \cr t^{\,Z}_{\rm{B}} & = {{-cX} \over {-X\cos(m)+iz\sin(m)}}. & (58b) }]$

Considering the equivalence of notations between this work and Zachariasen' book (see Table 1), we can verify that equations (55) are identical to (34) and, similarly, equations (56) are identical to (39).

Table 1
Correspondences of notation in this work and Zachariasen (1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ )

Zachariasen	This work
$[\exp(-2\pi i\,{\bf{k}}_{0}.{\bf{r}})]$	$[\exp\left(i\,{\bf{k}}_{0}.{\bf{r}}\right)]$
k₀ = 1/λ	k = 2π/λ
α_Z	−α
Ψ₀	χ₀
Ψ_H	χ_h
z	−(λ/π)ω
X	(λ/π)a
t₀	t_c = T/γ₀
m	aT
c	$[\exp[iT(u_{0}+\omega)]]$

Footnotes

¹If the approximation in equation (9) were not used, the coefficients of D_0,h in (10b) should be multiplied by (1 − α)^−1/2 ≃ 1 + α/2 +…; since χ₀, χ_h and α are much smaller than 1, we neglect the high-order terms.

²For any point $[{\bf{r}}]$ = (s₀,s_h,s_t), s is the path length inside the crystal along $[\hat{{\bf{s}}}_{0}]$ : the ray through r enters the crystal at the point of coordinates (s′₀, s_h, s_t) such that γ₀s′₀ + γ_hs_h = 0, so that s = s₀ − s′₀ = s₀ + s_h/b.

³These solutions are the Ewald wavefields discussed in detail by Authier (2003 $[Authier, A. (2003). Dynamical Theory of X-ray Diffraction. Oxford University Press.]$ ) using the dispersion surface.

⁴Our α has the opposite sign to the α defined in equation [3.114b] of Zachariasen (1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ).

⁵To see this, suppose $[\sin\theta-\sin\theta_{\rm{c}}]$ > 0; this implies [equation (27)] $[{\rm{Re}}\,\omega]$ < 0 (note that b < 0 for the Bragg case) and [equation (28a)] $[{\rm{Re}}\,a^{2}]$ > 0, therefore [equation (29)] $[{\rm{Re}}\,a]$ > 0 if $[{\rm{Im}}\,a^{2}]$ > 0. Similarly, supposing $[\sin\theta-\sin\theta_{\rm{c}}]$ < 0, we obtain $[{\rm{Re}}\,\omega]$ > 0 and $[{\rm{Re}}\,a^{2}]$ < 0, therefore $[{\rm{Re}}\,a]$ < 0 if $[{\rm{Im}}\,a]$ > 0.

⁶The multiplication should be done from bottom to top, i.e. M = M_nM_n−1…M₂M₁.

⁷Note that in equation (56b) we write (c₂ − c₁) rather than (c₁ − c₂) in Zachariasen's equation [3.137]. This does not affect the result shown by Zachariasen as the amplitudes are squared to give intensities. However, for calculating the amplitudes, the correct sign (as shown here) is important.

Acknowledgements

We recognize the work of Edoardo Cappelli and Mark Glass who developed the first version of crystalpy implementing the theory in Zachariasen (1994 $[Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications.]$ ). We acknowledge Ali Khounsary and Yujia Ding (CSRII, Illinois Institute of Technology, Chicago, IL 60616, USA) for helpful discussions on stacked crystals.

References

Authier, A. (2003). Dynamical Theory of X-ray Diffraction. Oxford University Press. Google Scholar
Baumgärtel, P., Grundmann, P., Zeschke, T., Erko, A., Viefhaus, J., Schäfers, F. & Schirmacher, H. (2019). AIP Conf. Proc. 2054, 060034. Google Scholar
Bouchenoire, L., Brown, S. D., Thompson, P., Duffy, J. A., Taylor, J. W. & Cooper, M. J. (2003). J. Synchrotron Rad. 10, 172–176. Web of Science CrossRef CAS IUCr Journals Google Scholar
Chernikov, R. & Klementiev, K. (2017). Proc. SPIE, 10388, 1038806. Google Scholar
Chubar, O. & Elleaume, P. (1998). Proceedings of the 6th European Particle Accelerator Conference (EPAC1998), pp. 1177–1179. THP01G. Google Scholar
Detlefs, C., Sanchez del Rio, M. & Mazzoli, C. (2012). Eur. Phys. J. Spec. Top. 208, 359–371. Web of Science CrossRef CAS Google Scholar
Ishikawa, T., Tamasaku, K. & Yabashi, M. (2005). Nucl. Instrum. Methods Phys. Res. A, 547, 42–49. Web of Science CrossRef CAS Google Scholar
Klementiev, K. & Chernikov, R. (2023). Synchrotron Radiat. News, 36(5), 23–27. CrossRef Google Scholar
Pinsker, Z. G. (1978). Dynamical Scattering of X-rays in Crystals. Springer-Verlag. Google Scholar
Rebuffi, L. & Sanchez del Rio, M. (2017). Proc. SPIE, 10388, 103880S. Google Scholar
Ren, B., Dilmanian, F., Chapman, L., Ivanov, I., Wu, X., Zhong, Z. & Huang, X. (1999). Nucl. Instrum. Methods Phys. Res. A, 428, 528–550. Web of Science CrossRef CAS Google Scholar
Rovezzi, M., Lapras, C., Manceau, A., Glatzel, P. & Verbeni, R. (2017). Rev. Sci. Instrum. 88, 013108. Web of Science CrossRef PubMed Google Scholar
Sanchez del Rio, M. (2021). DABAX Python library, https://github.com/oasys-kit/dabax. Google Scholar
Sanchez del Rio, M., Canestrari, N., Jiang, F. & Cerrina, F. (2011). J. Synchrotron Rad. 18, 708–716. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sanchez del Rio, M., Perez-Bocanegra, N., Shi, X., Honkimäki, V. & Zhang, L. (2015). J. Appl. Cryst. 48, 477–491. Web of Science CrossRef CAS IUCr Journals Google Scholar
Schäfers, F. (2008). The BESSY Raytrace Program RAY, pp. 9–41. Berlin, Heidelberg: Springer. Google Scholar
Schoonjans, T., Brunetti, A., Golosio, B., Sanchez del Rio, M., Solé, V. A., Ferrero, C. & Vincze, L. (2011). At. Spectrosc. 66, 776–784. Web of Science CrossRef CAS Google Scholar
Shvyd'ko, Y. (2004). X-ray optics: high-energy-resolution applications, Vol. 98 of Springer Series in Optical Sciences. Springer Science & Business Media. Google Scholar
Stepanov, S. A. (2004). Proc. SPIE, 5536, 16–26. CrossRef Google Scholar
Suortti, P., Thomlinson, W., Chapman, D., Gmür, N., Siddons, D. & Schulze, C. (1993). Nucl. Instrum. Methods Phys. Res. A, 336, 304–309. CrossRef Web of Science Google Scholar
Sutter, J. P., Chubar, O. & Suvorov, A. (2014). Proc. SPIE, 9209, 92090L. Google Scholar
Sutter, J. P., Pittard, J., Filik, J. & Baron, A. Q. R. (2022). J. Appl. Cryst. 55, 1011–1028. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sutter, J. P., Pittard, J., Filik, J. & Baron, A. Q. R. (2023). Proc. SPIE, 12697, 126970A. Google Scholar
Takagi, S. (1962). Acta Cryst. 15, 1311–1312. CrossRef CAS IUCr Journals Web of Science Google Scholar
Taupin, D. (1964). Bull. Soc. Fr. Miner. Crystallogr. 87, 469–511. CAS Google Scholar
Taupin, D. (1967). Acta Cryst. 23, 25–35. CrossRef CAS IUCr Journals Web of Science Google Scholar
Tolentino, H., Dartyge, E., Fontaine, A. & Tourillon, G. (1988). J. Appl. Cryst. 21, 15–22. CrossRef CAS Web of Science IUCr Journals Google Scholar
Yu, X. J., Chi, X., Smulders, T., Wee, A. T. S., Rusydi, A., Sanchez del Rio, M. & Breese, M. B. H. (2022). J. Synchrotron Rad. 29, 1157–1166. Web of Science CrossRef CAS IUCr Journals Google Scholar
Zachariasen, W. H. (1994). Theory of X-ray Diffraction in Crystals. Dover Publications. Google Scholar

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

JOURNAL OF
SYNCHROTRON
RADIATION

ISSN: 1600-5775

Volume 31| Part 6| November 2024| Pages 1469-1480

https://doi.org/10.1107/S160057752400924X

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}}\)

Formulation of perfect-crystal diffraction from Takagi–Taupin equations: numerical implementation in the crystalpy library

1. Introduction

2. Takagi–Taupin equations

3. Solutions of TT equations for a plane wave incident on a crystal with plane entrance surface

3.1. Expressions of α, ω and as functions of the angles

3.2. The transfer matrix of a parallel crystal slab

3.3. The transfer matrix for the case of a `thick crystal'

3.4. Reflection and transmission amplitudes in the Laue case

3.5. Reflection and transmission amplitudes in the Bragg case – the S-matrix

3.5.1. Reflection amplitude for a thick absorbing crystal in the Bragg case

3.5.2. Reflection amplitude for non-absorbing crystals in the Bragg case

3.6. Calculation of reflection and transmission amplitudes using the transfer matrix

3.7. The direction of the diffracted wave in vacuum

4. The crystalpy library

5. Summary and conclusions

APPENDIX A

Derivation of the TT equations for a rotating perfect crystal

A1. Demonstration of α = (2/k)(∂ϕ/∂sh)

APPENDIX B

Solutions of TT equations (17) using the Laplace transform

B1. Laue solution based on Laplace transform

B2. Bragg solution based on Laplace transform

APPENDIX C

Equivalence of amplitudes in equation (39) and (34) with Zachariasen's formulation

Footnotes

Acknowledgements

References

research papers

3.1. Expressions of α, ω and $[{\bi{a}}]$ as functions of the angles

A1. Demonstration of α = (2/k)(∂ϕ/∂s_h)