A finite difference scheme for integrating the Takagi–Taupin equations on an arbitrary orthogonal grid

Carlsen, M.; Simons, H.

doi:10.1107/S2053273322004934

research papers

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 78| Part 5| September 2022| Pages 395-401

https://doi.org/10.1107/S2053273322004934

Open

access

A finite difference scheme for integrating the Takagi–Taupin equations on an arbitrary orthogonal grid

Mads Carlsen ^a ^* and Hugh Simons ^a

^aDepartment of Physics, Technical University of Denmark (DTU), Fysikvej, Building 311, 2800 Kgs. Lyngby, Denmark
^*Correspondence e-mail: madsac@dtu.dk

Edited by I. A. Vartaniants, Deutsches Electronen-Synchrotron, Germany (Received 1 March 2022; accepted 9 May 2022; online 8 July 2022)

Calculating dynamical diffraction patterns for X-ray diffraction imaging techniques requires numerical integration of the Takagi–Taupin equations. This is usually performed with a simple, second-order finite difference scheme on a sheared computational grid in which two of the axes are aligned with the wavevectors of the incident and scattered beams. This dictates, especially at low scattering angles, an oblique grid of uneven step sizes. Here a finite difference scheme is presented that carries out this integration in slab-shaped samples on an arbitrary orthogonal grid by implicitly utilizing Fourier interpolation. The scheme achieves the expected second-order convergence and a similar error to the traditional approach for similarly dense grids.

Keywords: dynamical diffraction; X-ray topography.

1. Introduction

Simulations based on the propagation of coherent wavefronts are becoming an increasingly common tool for the development of X-ray diffraction imaging techniques, where they are regularly used to evaluate the viability of new methods (Pedersen et al., 2018 ; Holstad et al., 2022 ) and to investigate the effect of experimental errors (Shabalin et al., 2017 ; Carnis et al., 2019 ). Furthermore, by accurately predicting coherent interference, such simulation methods are particularly relevant given the recent arrival of highly coherent X-ray sources, such as fourth-generation synchrotrons and free-electron lasers.

When simulating the propagation of coherent wavefronts through large and near-perfect single crystals, multiple scattering effects (i.e. dynamical diffraction) become important. One often tries to avoid these dynamical effects [even if occasionally they are the subject of interest (Rodriguez-Fernandez et al., 2021 )] by using highly deformed samples, small grains or relying on the `weak beam approximation', i.e. measuring at the tails of the rocking curve (Shabalin et al., 2017; Holstad et al., 2022). However, in many cases dynamical effects are unavoidable and must be accounted for in the simulation framework by solving the Takagi–Taupin equations (TTEs): a set of coupled, first-order PDEs (partial differential equations) that, in general, must be integrated numerically (Takagi, 1962 ; Taupin, 1967 ). When numerically integrating the TTEs (in the two-beam case), it is natural to choose a computational grid with two of its axes aligned with the wavevector of the incident and scattered waves (Taupin, 1967; Authier et al., 1968 ). With this approach, the TTEs have previously been solved via finite difference integration on a structured grid of constant (Authier et al., 1968) or varying (Epelboin, 1981 ) step sizes, or by iterative approaches (Bremer, 1984 ; Yan & Li, 2014 ). However, the use of sheared grids may complicate matters by requiring the use of a connecting interpolation step when the scattering calculation is combined with other numerical methods, for example to generate input for the incident wavefront, to generate the input for the deformed crystal structure, or when the calculated diffraction patterns are further input into simulations of the downstream optics. Approaches using different coordinate systems have typically involved the use of symmetric scattering geometries in which the sheared coordinate system coincides with a rectangular one (Kolosov & Punegov, 2005 ; Osterhoff, 2012 ), or a finite element approach that can solve the TTEs on an unstructured grid (Honkanen et al., 2018 ). These approaches, however, are either geometrically restrictive or require third-party software.

Ideally, one should be able to straightforwardly integrate the TTEs on an orthogonal grid that requires no intermediate interpolation step. In the kinematical case (i.e. not incorporating multiple scattering), Li et al. (2020 ) described a method for carrying out scattering simulations on an orthogonal grid that implicitly uses Fourier interpolation to avoid making cumulative interpolation errors that would otherwise cause such a calculation to fail. Inspired by this approach, we here describe how implicit Fourier interpolation can also be utilized to integrate the dynamical TTEs on an orthogonal grid. Our approach is based on exponential Rosenbrock-type methods (Hochbruck & Ostermann, 2010 ) for the numerical integration, and yields a a result similar to the mixed real-space/reciprocal-space methods called `multistep methods' regularly used to model a wide range of optical scenarios (Li et al., 2017 ; Hare & Morrison, 1994 ). The method we present is applicable for slab-shaped samples (two parallel surfaces and infinite extent in the orthogonal directions) in Laue geometry, making it ideal for simulating X-ray diffraction images from lightly deformed (i.e. strained) materials.

2. The Takagi–Taupin equations

The most general framework for treating dynamical diffraction from strained crystals involves the TTE (Takagi, 1962, 1969 ; Taupin, 1967) which, for the two-beam case, are

$[\eqalignno{ 2i({\bf k}_{0}\cdot\nabla)E_{0}({\bf r})& = k^{2}\left[\chi_{0}E_{0}({\bf r})+\chi_{\overline{h}}^{ \prime}({\bf r})E_{h}({\bf r})\right]&\cr 2i({\bf k}_{h}\cdot\nabla)E_{h}({\bf r})& = k^{ 2}\left[(\chi_{0}+\beta)E_{h}({\bf r})+\chi_{h}^{\prime}({\bf r})E_{0}({\bf r})\right],&(1)}]$

where E₀ and E_h are the complex envelopes of the monochromatic fields of the incident and scattered beams, respectively, $[{\bf k}_{0}]$ and $[{\bf k}_{h} = {\bf k}_{0}+{\bf Q}]$ are the vacuum wavevectors of the incident and scattered beams, respectively,¹ and $[{\bf Q}]$ is the scattering vector. $[\chi_{0}]$ is the average electric susceptibility of the crystal, while $[\chi_{h}^{\prime}]$ and $[\chi_{\overline{h}}^{\prime}]$ are the spatially varying Fourier components of the electric susceptibility corresponding to the scattering vectors $[{\bf Q}]$ and $[{\bf -Q}]$ , respectively, given by

$[\eqalignno{\chi_{h}^{\prime}({\bf r})& = \exp\left[i{\bf Q}\cdot{\bf u}({\bf r})\right]\chi_{h},&\cr \chi_{\overline{h}}^{\prime}({\bf r})& = \exp\left[-i{\bf Q}\cdot{\bf u}({\bf r})\right]\chi_{\overline{h}},&(2)}]$

where $[{\bf u}({\bf r})]$ is the displacement field of the crystal. These susceptibility terms are related to the form factors F_h and $[F_{\overline{h}}]$ through

$[\eqalignno{\chi_{h} &= -\left({{4\pi r_{0}} \over {k^{2}V_{{\rm u. c.}}}}\right)F_{h},&\cr \chi_{\overline{h}} &= -\left({{4\pi r_{0}} \over {k^{2}V_{{\rm u.c. }}}}\right)F_{\overline{h}},&(3)}]$

where r₀ is the classical electron radius and $[V_{{\rm u.c.}}]$ is the volume of the crystal unit cell. Finally, $[\beta = (|{\bf k}_{h}|^{2}-k^{2})/k^{2}]$ is a measure of the deviation away from the exact Bragg condition, where k is the magnitude of $[{\bf k}_{0}]$ . If we consider only a right-handed `rocking' rotation around an axis parallel to $[{\bf Q}\times{\bf k}_{0}]$ (` ×' denotes the cross-product) by an angle ϕ, we can write $[\beta = 2\sin(2\theta)\phi]$ , where $[\phi = 0]$ corresponds to the exact satisfaction of the Bragg condition. Rotation of the crystal can also be simulated by the addition of a rotational component to the displacement field (e.g. Shabalin et al., 2017). This alternative approach is also compatible with the method presented here.

If we ignore the scattering terms, the equations (1) are a pair of convection equations whose solution involves the interpolation of the initial condition through the integration volume. Since direct application of a finite difference scheme in a Cartesian coordinate system is equivalent to linear interpolation and would accumulate errors at each step, the traditional approach is to solve the equations in an oblique coordinate system with the axes aligned with the incident and scattered wavevectors. In this way, the interpolation from slice to slice in the computational grid becomes a simple shifting of array elements. Though simple and elegant, the use of such oblique coordinate systems requires additional interpolation steps when the scattering simulation interfaces with material models or other optical simulation steps that require orthogonal grids. The approach we describe here utilizes an orthogonal grid, but avoids cumulative interpolation errors by re-stating the problem in reciprocal space and utilizing Fourier interpolation implicitly.

3. Rewriting the TTEs in terms of transverse Fourier transforms

We wish to arrive at a mixed real-/reciprocal-space solution to the TTEs. To achieve this, we first state the TTEs in an orthogonal coordinate system, then Fourier transform the resulting differential equations along the transverse coordinates.

We define an orthogonal grid with the three orthonormal unit vectors $[\hat{{\bf x}}]$ , $[\hat{{\bf y}}]$ and $[\hat{{\bf z}}]$ . The only restriction on the choice of coordinate system is that

$[{\bf k}_{0}\cdot\hat{{\bf z}}\,\gt\,0\quad{\rm and}\quad{\bf k}_{h}\cdot \hat{{\bf z}}\,\gt\,0,\eqno(4)]$

such that the z axis takes the role of a quasi-optical axis and we can treat z as the dynamical variable and x and y as transverse variables. To this end, we decompose the vectors $[{\bf k}_{0}]$ and $[{\bf k}_{h}]$ into their z components and their projection onto the x–y plane, i.e. $[k_{0} = k_{0,z}\hat{{\bf z}}+{\bf k}_{0,\perp}]$ and $[k_{h} = k_{h,z}\hat{{\bf z}}+{\bf k}_{h,\perp}]$ , where k_0,z,k_h,z refer to the z components of the respective vectors. We can now rewrite equations (1) as

$[\eqalignno{ 2k_{0,z}{{\partial} \over {\partial z}}E_{0}({\bf r})& = -ik^{2}\chi_{0}E_{0}({\bf r})-2({\bf k}_{0,\perp}\cdot \nabla_{\perp})E_{0}({\bf r})&\cr &\quad -ik^{2}\chi_{\overline{h}}^{\prime}({\bf r })E_{h}({\bf r}),&\cr 2k_{h,z}{{\partial} \over {\partial z}}E_{h}({\bf r})& = -ik^{2}(\chi_{0}+\beta)E_{h}({\bf r})-2({\bf k}_{h,\perp} \cdot\nabla_{\perp})E_{h}({\bf r})&\cr &\quad -ik^{2}\chi_{h}^{\prime}({\bf r})E_{0} ({\bf r}),&(5)}]$

where $[\nabla_{\perp} = [\partial/\partial x,\partial/\partial y,0]]$ . We now introduce the transverse Fourier transform:

$[{\cal F}_{\perp}\{E(x,y,z)\} = \textstyle\int\int E(x,y,z)\exp[-i2\pi(xq_{x}+yq_{y})]\, {\rm d}x\,{\rm d}y,\eqno(6)]$

which is inverted by the inverse Fourier transform, appropriately defined as

$[{\cal F}^{-1}_{\perp}\{\tilde{E}(q_{x},q_{y},z)\} = \textstyle\int\int\tilde{E}(q_{x},q _{y},z)\exp[i2\pi(xq_{x}+yq_{y})]\,{\rm d}x\,{\rm d}y.\eqno(7)]$

Tildes are used to denote the transforms of functions, namely $[\tilde{E}_{0}(q_{x},q_{y},z)]$ = $[{\cal F}_{\perp}\{\tilde{E}_{0}(x,y,z)\}]$ and $[\tilde{E}_{h}(q_{x},q_{y},z)]$ = $[{\cal F}_{\perp}\{\tilde{E}_{h}(x,y,z)\}]$ .

With this definition, we Fourier transform the equations (5):

$[\eqalignno{{{\partial} \over {\partial z}}\tilde{E}_{0}(q_{x},q_{y}, z)& = \left[{{-ik^{2}} \over {2k_{0,z}}}\chi_{0}-{{i2\pi} \over {k_{0,z }}}{\bf q}\cdot{\bf k}_{0,\perp}\right]\tilde{E}_{0}(q_{x},q_{y},z)&\cr &\quad-{{ik^{2}} \over {2k_{0,z}}}{\cal F}_{\perp}\lbrace\chi_{ \overline{h}}^{\prime}(x,y,z)E_{h}(x,y,z)\rbrace,&\cr {{\partial} \over {\partial z}}\tilde{E}_{h}(q_{x},q_{y},z)& = \left[{{-ik^{2}} \over {2k_{h,z}}}(\chi_{0}+\beta)-{{i2\pi} \over {k_ {h,z}}}{\bf q}\cdot{\bf k}_{0,\perp}\right]\tilde{E}_{h}(q_{x},q_{y},z) &\cr &\quad-{{ik^{2}} \over {2k_{h,z}}}{\cal F}_{\perp}\lbrace\chi_{h}^{ \prime}(x,y,z)E_{0}(x,y,z)\rbrace,&(8)}]$

where $[{\bf q} = [q_{x},q_{y},0]]$ . Here we have assumed that $[\chi_{0}]$ is constant throughout the simulated volume, which limits this approach to slab-shaped crystals. We introduce the angles $[\alpha_{0}]$ and $[\alpha_{h}]$ given by $[{\bf k}_{0}\cdot\hat{{\bf z}}]$ = $[|{\bf k}_{0}|\cos(\alpha_{0})]$ and $[{\bf k}_{h}\cdot\hat{{\bf z}}]$ = $[|{\bf k}_{h}|\cos(\alpha_{h})]$ to give

$[\eqalignno{{{\partial} \over {\partial z}}\tilde{E}_{0}(q_{x},q_{y}, z)& = \left[{{-ik} \over {2\cos(\alpha_{0})}}\chi_{0}-{{i2\pi} \over { \cos(\alpha_{0})}}q_{0,\perp}\right]\tilde{E}_{0}(q_{x},q_{y},z)&\cr &\quad-{{ik} \over {2\cos(\alpha_{0})}}{\cal F}_{ \perp}\lbrace\chi_{\overline{h}}^{\prime}(x,y,z)E_{h}(x,y,z)\rbrace,&\cr {{\partial} \over {\partial z}}\tilde{E}_{h}(q_{x},q_{y},z)& = \left[{{-ik} \over {2\cos(\alpha_{h})}}\chi_{0}-{{i2\pi} \over {\cos(\alpha_{h})}}q_{h,\perp}\right]\tilde{E}_{h}(q_{x},q_{y},z)&\cr &\quad-{{ik} \over {2\cos(\alpha_{h})}}{\cal F}_{ \perp}\lbrace\chi_{h}^{\prime}(x,y,z)E_{0}(x,y,z)\rbrace &(9)}]$

where $[q_{0,\perp} = {\bf q}\cdot{\bf k}_{0,\perp}/k]$ and $[q_{h,\perp} = {\bf q}\cdot{\bf k}_{h,\perp}/k]$ .

In cases where $[\chi_{h}^{\prime}]$ is constant or depends only on z, the equations can be solved analytically with Green's function methods. In the general case where $[\chi_{h}^{\prime}]$ varies as a function of all coordinates, the scattering term $[{\cal F}_{\perp}\{\chi_{\overline{h}}^{\prime}(x,y,z)E_{h}(x,y,z)\}]$ cannot be simplified and finite difference methods must be used.

We note that in cases when both $[{\bf k}_{0}]$ and $[{\bf k}_{h}]$ lie within the x–z plane, the 2D Fourier transforms may be replaced by 1D Fourier transforms along the x direction.

4. Discretization and numerical evaluation

We introduce a computational grid with axes parallel to those in the coordinate system defined above. It has step sizes d_x, d_y and d_z, and number of points N_x, N_y and N_z in each dimension. The thickness of the slab-shaped crystal is $[t = d_{z}(N_{z}-1)]$ and the simulated volume has side lengths $[L_{x} = d_{x}(N_{x}-1)]$ and $[L_{y} = d_{y}(N_{y}-1)]$ in the transverse directions. A point on the grid P = (i_xd_x,i_yd_y,i_zd_z) is then indexed by the numbers i_x, i_y, i_z where i_z = 0, 1, 2… $[N_{z}-1]$ .

In order to utilize discrete Fourier transform methods when solving these equations on a finite grid, we impose zero Dirichlet boundary conditions in the two transverse dimensions, x and y:

$[\eqalignno{ E_{0}(0,y,z)& = E_{0}(L_{x},y,z) = E_{h}(0,y,z) = E_{h}(L_ {x},y,z) = 0,&\cr E_{0}(x,0,z)& = E_{0}(x,L_{y},z) = E_{h}(x,0,z) = E_{h}(x,L_{y},z) = 0.&\cr &&(10)}]$

These boundary conditions require that the sample grid be large enough to fit the Borrmann triangle extending from every point where the initial condition is non-zero. If the initial condition is only non-zero on a domain Ω on the surface (i.e. at z = 0), the direct projection of this domain along the directions of $[{\bf k}_{0}]$ and $[{\bf k}_{h}]$ must lie within the sample grid (see Fig 1). This condition is fulfilled if the domain Ω is fully contained in the rectangle defined by

$[\eqalignno{&\max(0,tk_{0x}/k,tk_{hx}/k)\,\lt\,x&\cr &\lt\,\min(L_{x},L_{x}+tk_{0x}/k,L_{x}+tk _{hx}/k),&(11)}]$

$[\max(0,tk_{0x}/k,tk_{hx}/k)\,\lt\,y\,\lt\,\min(L_{y},L_{y}+tk_{0y}/k,L_{y}+tk _{hy}/k)\eqno(12)]$

which can always be made true for a finitely bounded initial condition if the computational grid is sufficiently large.

Figure 1
Scattering geometry inside the sample volume and the finite support of the initial condition.

As we are using discrete Fourier transforms we also specify a grid in ( q_x,q_y) space, which is related to the (x–y) grid in real space. In full period frequency, this corresponds to the points

$[\eqalignno{ q_{x}& = \left[-{{1} \over {2d_{x}}},-{{ 1} \over {2d_{x}}}+{{1} \over {L_{x}}},-{{1} \over {2d_{x}}}+{{2} \over {L_{x}}}\ldots,0,\ldots,{{1} \over {2d_{x}}}-{{1} \over {L_{x}}}\right],&\cr q_{y}& = \left[-{{1} \over {2d_{y}}},-{{1} \over {2d_{y}}}+ {{1} \over {L_{y}}},-{{1} \over {2d_{y}}}+{{2} \over {L_{y}}}\ldots,0,\ldots,{{1} \over {2d_{y}}}- {{1} \over {L_{y}}}\right],&\cr &&(13)}]$

where the inclusion of the negative Nyquist frequency is specific to even numbers for N_x and N_y.

The complex envelope of the incident beam should be given for the entrance surface of the grid [i.e. E₀(d_xi_x,d_yi_y,0)] and the components transformed using a discrete Fourier transform (DFT) to give the initial condition in reciprocal space where the amplitude of the scattered beam is zero at the entrance surface, i.e. $[\tilde{E}_{h}(q_{x},q_{y},0) = 0]$ . This initial condition is constructed by appending the components of these two arrays into a single vector: $[{\bf E}(0)]$ .

We can rewrite the equations (9) to the form used by exponential Runge–Kutta methods: $[({\partial} / {\partial z}){\bf E}(z)]$ = $[A{\bf E}(z)]$ + $[{\bf B}[z,{\bf E }(z)]]$ , where A is a diagonal matrix containing the coefficients in the square brackets of equation (9) evaluated at the reciprocal-space grid defined by equation (13) and $[{\bf B}]$ contains the convolution terms, i.e. the last terms of equations (9). The function is implemented by applying an inverse DFT on $[{\bf E}]$ , multiplying by the scattering coefficients, that must be given on the real-space grid at the slice z, and transforming back to reciprocal space:

$[\eqalignno{{\bf B}(z,{\bf E})& = \Bigg[-{{ik} \over {2\cos(\alpha_{0})}} {\rm DFT}\left\lbrace{\rm DFT}^{-1}\lbrace\tilde{E}_{h}\rbrace\chi^{\prime}_{\overline{ h}}(x,y,z)\right\rbrace &\cr & -{{ik} \over {2\cos(\alpha_{h})}}{\rm DFT}\left\lbrace{\rm DFT}^{-1} \lbrace\tilde{E}_{0}\rbrace\chi^{\prime}_{h}(x,y,z)\right\rbrace\Bigg].&(14)}]$

5. Exponential Runge methods and convergence behaviour

The Fourier-transformed TTEs can be solved by exponential Runge–Kutta methods of the type given by Hochbruck & Ostermann (2010). To test the convergence of this approach, we utilized two different exponential integrators. The first is an archetypal exponential integrator based on the explicit Euler scheme, given by

$[{\bf E}(z+h) = \exp(hA){\bf E}(z)+h(hA)^{-1}[\exp(hA)-1]{\bf B}\left[z, {\bf E}(z)\right].\eqno(15)]$

The second exponential integrator is an explicit second-order method based on Heun's method, given by the steps (Friedli, 1978 )

$[\eqalignno{ E^{*}_{1}& = {\bf E}(z),&\cr b_{1}& = {\bf B}(z_{n},E_{1}^{*}),&\cr E^{*}_{2}& = \phi_{0}E_{1}^{*}+h\phi_{1}b_{1},&\cr b_{2}& = {\bf B}(z_{n}+h,E_{2}^{*}),&\cr {\bf E}(z+h)& = \phi_{0}E_{1}^{*}+{{h} \over {2}}[(2 \phi_{1}-\phi_{2})b_{1}+\phi_{2}b_{2}],&(16)}]$

where the ϕ functions are given by $[\phi_{0} = \exp(hA)]$ and $[\phi_{n} = n(hA)^{-1}(\phi_{n-1}-1)]$ , and exp(hA) is in general the matrix-exponential which here is equal to the element-wise exponential because A is diagonal. We choose this scheme based on Heun's method because it only evaluates the B function on the same regular intervals where the field is calculated, and therefore only requires the value of the scattering function on the same grid where the fields are evaluated.

For comparison with existing methods, we also implemented a normal finite difference method based on a recent publication by Shabalin et al. (2017) using the half-step finite difference for the derivatives. A derivation of this method is given in Appendix A.

To evaluate the convergence behaviour of these exponential methods, we generated a virtual sample consisting of a perfect single crystal with a single edge dislocation with Burger's vector (100) close to the path of the direct beam. Plots of the displacement field as well as the amplitudes of the converged solution are shown in Fig. 2.

Figure 2
Plots of the sample and calculated fields used in the second convergence test. (a) Displacement field in units of the lattice constant, a. (b) Transmitted field on a logarithmic scale. (c) Scattered field on a logarithmic scale. The calculated fields are for the case $[\beta = 0]$ .

The fields are simulated under low absorption and highly dynamical conditions, and we simulated only a single slice in the y direction with the dimensions 50 × 115 µm at a point 1 µm from the dislocation core. The incident beam is a narrow Gaussian of width σ = 0.2 µm, with parameters chosen corresponding to the (111) reflection of a diamond crystal with a molybdenum source, given as: $[\chi_{0}]$ = $[7.6\times 10^{-6}-1.4i\times 10^{-9}]$ , $[\chi_{h} = 5.0\times 10^{-6}-\,0.7i\times 10^{-9}]$ and $[\chi_{\overline{h}} = \chi_{h}]$ , where i is the imaginary unit. We set k = 8.85 Å⁻¹ and $[2\theta]$ = 20°. The incident beam has a Gaussian envelope with width σ_x = 0.2 µm.

In order to accommodate the comparison with existing methods, we utilized a grid with step sizes $[\Delta z = h]$ and $[\Delta x = 2\tan(\theta)h]$ for the exponential methods and a grid with the same density of points for the normal finite difference method. Then, to check the convergence of the methods, we calculate the fields on progressively finer grids from a first grid consisting of 101 × 41 points. The error represents the deviation of the exit surface amplitudes to that computed using the normal finite difference approach on a very fine grid of 10 241 × 25 601 steps (evaluated at the points where the coarse and fine grids coincide).

Fig. 3 shows the convergence of the three integration methods. While all methods show the expected convergence on a perfect sample (a), the traditional half-step method does not show the expected second-order convergence with the edge dislocation sample (b). The first-order exponential Euler method suffers from an exponential instability and only gives a qualitatively correct result when impractically small step sizes are utilized.

Figure 3
Convergence of the new exponential integrators and a traditional finite difference scheme. The black lines mark first- and second-order convergence. All errors are calculated relative to the solution using the traditional half-step method with 10 241 steps. We tested integration schemes on two different samples. One (a) is a perfect crystal, the other (b) is the edge dislocation type sample shown in Fig. 2

. All calculations are in the case $[\beta = 0]$ .

To demonstrate that our method is not limited only to the exact Bragg condition, we also test the convergence of our method at the tails of the rocking curve. Fig. 4(a) shows the convergence of the screw-dislocation test case at a rocking angle of $[\phi]$ = 300 µrad, where the two methods reach the expected first- and second-order convergence. At this point, the diffraction is approximately kinematical. Figs. 4(c), 4(d) show that the transmitted beam is approximately undisturbed by the crystal and that the beam is only significantly scattered in small areas close to the surface and the dislocation. Fig. 4(b) shows the profile of the scattered beam at the exit surface as a function of rocking angle. We see that at low angles, the scattering is clearly dynamical and the profile shows Pendellösung fringes. At higher angles, the bulk of the crystal scatters less strongly and the scattering pattern is dominated by the defects. Because the divergence of the incident radiation is larger than the characteristic `Darwin width' of dynamical diffraction, the crystal approximately acts as an analyser, and maps out the angular spectrum of the incident beam, which is further illustrated in Fig. 4(e).

Figure 4
Testing of the integration algorithms for finite rocking angles. (a) Convergence plot for the screw-dislocation test case at $[\phi]$ = 300 µrad. (b) Spatial profile of the scattered beam as a function of rocking angle. (c) Transmitted beam in a cross section of the crystal at $[\phi]$ = 300 µrad. (d) Scattered beam in a cross section of the crystal at $[\phi ]$ = 300 µrad. (e) Integrated rocking curve of the dislocation test case overlaid on the same curve calculated for a perfect crystal and the angular spectrum of the incident beam.

6. Discussion and conclusion

We have described and demonstrated a finite difference scheme capable of integrating the TTEs on an orthogonal grid with few restrictions on the choice of grid. We achieve this by implicitly utilizing Fourier interpolation at the level of the individual finite difference step. The method results in approximately the same error as the traditional half-step finite difference scheme.

The method utilizes FFTs (fast Fourier transforms) at each step and has to perform in total four 2D Fourier transforms of the entire sample volume (the Heun method makes two evaluations of the $[{\bf B}]$ function at each step), which is expected to have a computational cost compared with the existing methods. In certain geometries, these may be replaced by 1D Fourier transforms. Our experience here, using unoptimized code, is that this increase amounts to about a factor of 4, which we believe should be unimportant in most cases of practical interest.

The ability to freely choose the computational grid makes implementation of this approach easier, especially when it needs to be combined with other numerical modelling methods, for example if the input for either the crystal microstructure or the incident field is given by a numerical simulation, or if the scattered fields should be propagated through image-forming optics.

APPENDIX A

Traditional finite difference scheme

For comparison with the exponential integrators presented in this paper, we also present calculations performed with a second-order implicit finite difference scheme based on a central difference estimate for the derivatives, which is often applied for dynamical scattering calculations. The derivation here follows the one given by Shabalin et al. (2017) with small changes to the notation.

We limit our attention to a symmetric geometry defined by

$[{\bf k}_{0} = k\left[\matrix{\sin\theta\cr 0\cr \cos\theta}\right]\quad{\rm and}\quad{\bf k}_{h} = k\left[\matrix{- \sin\theta\cr 0\cr \cos\theta}\right].\eqno(17)]$

Starting from equation (5) we introduce the coordinates

$[\left[\matrix{s_{0}\cr s_{h}}\right] = \left[\matrix{\sin\theta&\cos\theta\cr -\sin\theta&\cos\theta}\right]\left[\matrix{x\cr z}\right]\eqno(18)]$

and arrive at a well known form of the TTE (suppressing the y dependence):

$[\eqalignno{{{\partial E_{0}(s_{0},s_{h})} \over {\partial s_{0}}}& = {{k} \over {2i}}\left[\chi_{0}E_{0}(s_{0},s_{h})+\chi_{\overline{h}}^{\prime}(s_{0 },s_{h})E_{h}(s_{0},s_{h})\right]&\cr {{\partial E_{h}(s_{0},s_{h})} \over {\partial s_{h}}}& = {{k} \over {2i}} \left[\chi_{0}E_{h}(s_{0},s_{h})+\chi_{h}^{\prime}(s_{0},s_{h})E_{0}(s_{0},s_{ h})\right].&(19)}]$

To avoid truncation errors due to the complex rotation caused by the $[\chi_{0}]$ terms, we introduce the scaled fields $[E_{0}^{\prime} = \exp[\chi_{0}({ik} / {2})s_{0}]]$ and $[E_{h}^{\prime} = \exp[\chi_{0}({ik} / {2})s_{h}]]$ . Plugging in and simplifying some terms gives:

$[\eqalignno{{{\partial E_{0}^{\prime}(s_{0},s_{h})} \over { \partial s_{0}}} &= {{k} \over {2i}}\exp\left[\chi_{0}{{ik} \over {2}}(s_{0}-s_{h})\right] \chi_{\overline{h}}^{\prime}(s_{0},s_{h})E_{h}^{\prime}(s_{0},s_{h})&\cr {{\partial E_{h}^{\prime}(s_{0},s_{h})} \over {\partial s_{h}}}& = {{k} \over {2i}}\exp\left[\chi_{0}{{ik} \over {2}}(s_{h}-s_{0})\right]\chi_{h}^{\prime }(s_{0},s_{h})E_{0}^{\prime}(s_{0},s_{h}).&\cr &&(20)}]$

We now introduce a rectangular grid in the original (x,y) coordinates with step size h in the z direction and $[h\tan(\theta)]$ in the x direction. With this choice of grid a subset consisting of every second grid point constitutes a sheared grid aligned with the s₀ and s_h directions with both step sizes equal to $[p = h[1+\tan(\theta)^{2}]^{1/2}]$ . This allows us to calculate the fields using the finite difference methods and the exponential integrators on grids with the same density of grid points and that coincide on every second plane. We therefore have to choose a grid with an odd number of grid points in the x direction so that we can compare the result on the final slice (see Fig. 5).

Figure 5
Computational grids used for the traditional half-step approach marked with red dots and for the exponential integrators marked with blue circles. The recurrence relation for the half-step method is drawn with green arrows.

We denote the discretized envelope fields by $[E^{\prime}(x_{j},z_{i}) = E^{i,j}]$ . The recurrence relation is obtained by the centred first-order approximation for the derivatives and a similar centred approximation for the right-hand sides in equation (20) to arrive at the equations

$[\eqalignno{{{E_{0}^{i,j}-E_{0}^{i-1,j-1}} \over {p}}& = {{k} \over {2i}}B {{E_{h}^{i,j}+E_{h}^{i-1,j-1}} \over {2}}&\cr {{E_{h}^{i,j}-E_{h}^{i-1,j+1}} \over {p}} &= {{k} \over {2i}}D{{E_{0}^{i, j}+E_{0}^{i-1,j+1}} \over {2}},&(21)}]$

where

$[B = \exp\left[\chi_{0}{{ik} \over {2}}(s^{i,j}_{0}-s^{i,j}_{h}-p/2)\right]\chi_{ \overline{h}}\left(s^{i,j}_{0}-p/2,s^{i,j}_{h}\right)]$

and

$[D = \exp\left[\chi_{0}{{ik} \over {2}}(s^{i,j}_{h}-s^{i,j}_{0}-p/2)\right]\chi_{h} \left(s^{i,j}_{0},s^{i,j}_{h}-p/2\right).]$

Introducing the constant A = 4i / kp we can finally write

$[\eqalignno{ E_{0}^{i,j}& = E_{0}^{i-1,j-1}+B/AE_{h}^{i,j}+B/AE_{h} ^{i-1,j-1}&\cr E_{h}^{i,j} &= E_{h}^{i-1,j+1}+D/AE_{0}^{i,j}+D/AE_{0}^{i-1,j+1} &(22)}]$

which are the implicit recurrence relations used in the calculations. Furthermore, we need the boundary conditions, that the fields are both zero at the top and bottom surfaces: $[E_{0/h}^{i,0} = E_{0/h}^{i,N_{x}-1} = 0]$ .

A1. Availability

A Python3 implementation of the described algorithm is available at https://github.com/Multiscale-imaging/dynamical_diffraction and is given as supporting information.

Supporting information

Python3 implementation of the described algorithm. DOI: https://doi.org/10.1107/S2053273322004934/iv5022sup1.exe

Footnotes

¹The choice of $[{\bf k}_{0}]$ is arbitrary, but leads to different versions of the TTEs; the other common choice is for $[{\bf k}_{0}]$ to be the wavevector of the refracted wave inside the crystal.

Funding information

MC and HS acknowledge funding from ERC Starting Grant No. 804665.

References

Authier, A., Malgrange, C. & Tournarie, M. (1968). Acta Cryst. 24, 126–136. CrossRef IUCr Journals Web of Science Google Scholar
Bremer, J. (1984). Acta Cryst. 40, 283–291. CSD CrossRef CAS Web of Science IUCr Journals Google Scholar
Carnis, J., Gao, L., Labat, S., Kim, Y. Y., Hofmann, J. P., Leake, S. J., Schülli, T. U., Hensen, E. J. M., Thomas, O. & Richard, M.-I. (2019). Sci. Rep. 9, 17357. Google Scholar
Epelboin, Y. (1981). Acta Cryst. A37, 132–133. CrossRef CAS IUCr Journals Web of Science Google Scholar
Friedli, A. (1978). Numerical Treatment of Differential Equations (Proc. Conf. Math. Forschungsinst. Oberwolfach, 1976), pp. 35–50. Lecture Notes in Mathematics, Vol. 631. Heidelberg: Springer. Google Scholar
Hare, A. R. & Morrison, G. R. (1994). J. Mod. Opt. 41, 31–48. CrossRef Web of Science Google Scholar
Hochbruck, M. & Ostermann, A. (2010). Acta Numer. 19, 209–286. CrossRef Google Scholar
Holstad, T. S., Raeder, T. M., Carlsen, M., Bergbäck Knudsen, E., Dresselhaus-Marais, L., Haldrup, K., Simons, H., Nielsen, M. M. & Poulsen, H. F. (2022). J. Appl. Cryst. 55, 112–121. Web of Science CrossRef CAS IUCr Journals Google Scholar
Honkanen, A.-P., Ferrero, C., Guigay, J.-P. & Mocella, V. (2018). J. Appl. Cryst. 51, 514–525. Web of Science CrossRef CAS IUCr Journals Google Scholar
Kolosov, S. & Punegov, V. (2005). Crystallogr. Rep. 50, 357–362. Web of Science CrossRef CAS Google Scholar
Li, K., Wojcik, M. & Jacobsen, C. (2017). Opt. Express, 25, 1831. Google Scholar
Li, P., Maddali, S., Pateras, A., Calvo-Almazan, I., Hruszkewycz, S. O., Cha, W., Chamard, V. & Allain, M. (2020). J. Appl. Cryst. 53, 404–418. Web of Science CrossRef CAS IUCr Journals Google Scholar
Osterhoff, M. (2012). Wave Optical Simulations of X-ray Nano-focusing Optics, Vol. 009 of Göttingen Series in X-ray Physics. Goettingen: Universitaetsverlag Goettingen. Google Scholar
Pedersen, A. F., Chamard, V. & Poulsen, H. F. (2018). Opt. Express, 26, 23411–23425. Web of Science CrossRef CAS PubMed Google Scholar
Rodriguez-Fernandez, A., Diaz, A., Iyer, A. H. S., Verezhak, M., Wakonig, K., Colliander, M. H. & Carbone, D. (2021). Phys. Rev. Lett. 127, 157402. Web of Science PubMed Google Scholar
Shabalin, A. G., Yefanov, O. M., Nosik, V. L., Bushuev, V. A. & Vartanyants, I. A. (2017). Phys. Rev. B, 96, 064111. Web of Science CrossRef Google Scholar
Takagi, S. (1962). Acta Cryst. 15, 1311–1312. CrossRef CAS IUCr Journals Web of Science Google Scholar
Takagi, S. (1969). J. Phys. Soc. Jpn, 26, 1239–1253. CrossRef CAS Web of Science Google Scholar
Taupin, D. (1967). Acta Cryst. 23, 25–35. CrossRef CAS IUCr Journals Web of Science Google Scholar
Yan, H. & Li, L. (2014). Phys. Rev. B, 89, 014104. Web of Science CrossRef Google Scholar

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

FOUNDATIONS
ADVANCES

ISSN: 2053-2733

Volume 78| Part 5| September 2022| Pages 395-401

https://doi.org/10.1107/S2053273322004934

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

A finite difference scheme for integrating the Takagi–Taupin equations on an arbitrary orthogonal grid

1. Introduction

2. The Takagi–Taupin equations

3. Rewriting the TTEs in terms of transverse Fourier transforms

4. Discretization and numerical evaluation

5. Exponential Runge methods and convergence behaviour

6. Discussion and conclusion

APPENDIX A

Traditional finite difference scheme

A1. Availability

Supporting information

Footnotes

Funding information

References

research papers