Dynamical effects in Bragg coherent x-ray diffraction imaging on finite crystals

TL;DR

This paper develops a numerical simulation approach based on dynamical diffraction theory to analyze artifacts in Bragg CXDI reconstructions of finite crystals, highlighting limitations of the kinematical approximation.

Contribution

It introduces a dynamical theory-based simulation method for Bragg CXDI of finite crystals, accounting for strain, shape, and artifacts, advancing understanding of diffraction effects.

Findings

Dynamical effects cause artifacts in Bragg CXDI reconstructions.

The developed model can simulate diffraction for arbitrary crystal shapes.

A criterion to distinguish between kinematical and dynamical regimes is proposed.

Abstract

We present simulations of Bragg Coherent X-ray Diffractive Imaging (CXDI) data from finite crystals in the frame of the dynamical theory of x-ray diffraction. The developed approach is based on numerical solution of modified Takagi-Taupin equations and can be applied for modeling of a broad range of x-ray diffraction experiments with finite three-dimensional crystals of arbitrary shape also in the presence of strain. We performed simulations for nanocrystals of a cubic and hemispherical shape of different sizes and provided a detailed analysis of artifacts in the Bragg CXDI reconstructions introduced by the dynamical diffraction. A convenient way to treat effects of refraction and absorption supported by analytical derivations is described. Our results elucidate limitations for the kinematical approach in the Bragg CXDI and suggest a natural criterion to distinguish between kinematical and dynamical cases in coherent x-ray diffraction on a finite crystal.

Equations

$$A(\mathbf{q})\propto \frac{F_h}{V_{u.c.}}\int S_h(\mathbf{r})e^{-i\mathbf{q}\cdot \mathbf{r}}d\mathbf{r}.$$

$$S_h(\mathbf{r})=s_h(\mathbf{r})e^{i{\phi }_h(\mathbf{r})},{\phi }_h(\mathbf{r})=-\mathbf{h}\cdot \mathbf{u}(\mathbf{r}),$$

$${\mathbf{k}}_h={\mathbf{k}}_0+\mathbf{h}{|}_{\Delta \theta =0}.$$

$$\mathbf{E}(\mathbf{r})=\sum _s\left[{\mathbf{e}}_{0s}{E}_{0s}(\mathbf{r})e^{i{\mathbf{k}}_0\cdot \mathbf{r}}+{\mathbf{e}}_{hs}{E}_{hs}(\mathbf{r})e^{i{\mathbf{k}}_h\cdot \mathbf{r}}\right],$$

$$\begin{array}{r}\frac{\partial E_{0s}(\mathbf{r})}{\partial s_0}=\frac{i\pi }{\lambda }[{\chi }_0{E}_{0s}(\mathbf{r})+C{\chi }_{\bar{h}}{e}^{-i\Delta \mathbf{q}\cdot \mathbf{r}+i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_{hs}(\mathbf{r})],\\ \frac{\partial E_{hs}(\mathbf{r})}{\partial s_h}=\frac{i\pi }{\lambda }[{\chi }_0{E}_{hs}(\mathbf{r})+C{\chi }_h{e}^{i\Delta \mathbf{q}\cdot \mathbf{r}-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_{0s}(\mathbf{r})],\end{array}$$

$$L_{ex}=\frac{\lambda \sqrt{{\gamma }_0|{\gamma }_h|}}{\pi Re\left[\sqrt{{\chi }_h{\chi }_{\bar{h}}}\right]},$$

$$E_0(\mathbf{r})=E_{in}(\mathbf{r}),\text{at}{s}_0=s_0^{in}\text{and}{E}_h(\mathbf{r})=0,\text{at}{s}_h=s_h^{in}.$$

$$\begin{array}{r}A(q_{\perp },q_y,\Delta \mathbf{q})=\iint E_h^{ESW}(s_{\perp },s_y,\Delta \mathbf{q})e^{-i{q}_{\perp }{s}_{\perp }-i{q}_y{s}_y}d{s}_{\perp }d{s}_y,\end{array}$$

$$S_h(\mathbf{r})=s_h(\mathbf{r})e^{i{\phi }_h(\mathbf{r})}=\frac{1}{{(2\pi )}^3}\int A(\mathbf{q})e^{i\mathbf{q}\cdot \mathbf{r}}d\mathbf{q},$$

$${\theta }_{ref}=\mp \frac{{\chi }_{0r}(1\pm \beta )}{2\beta \sin 2{\theta }_B},$$

$$\begin{array}{r}\frac{\partial E_0(\mathbf{r})}{\partial s_0}=\left(\frac{i\pi }{\lambda }\right){\chi }_0{E}_0(\mathbf{r}),\\ \frac{\partial E_h(\mathbf{r})}{\partial s_h}=\left(\frac{i\pi }{\lambda }\right)[{\chi }_0{E}_h(\mathbf{r})+{\chi }_h{e}^{i\Delta \mathbf{q}\cdot \mathbf{r}-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_0(\mathbf{r})].\end{array}$$

$$E_0(\mathbf{r})=E_0^{in}\exp \left[i\frac{{\chi }_0}{2}{\mathbf{k}}_0(\mathbf{r}-{\mathbf{R}}_{in})\right],$$

$$\frac{\partial E_h(\mathbf{r})}{\partial s_h}=\left(\frac{i\pi }{\lambda }\right)[{\chi }_0{E}_h(\mathbf{r})+{\chi }_h{e}^{i\Delta \mathbf{q}\cdot \mathbf{r}-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{e}^{i\frac{{\chi }_0}{2}{\mathbf{k}}_0(\mathbf{r}-{\mathbf{R}}_{in})}].$$

$$E_h(\mathbf{r})=E_h^{\prime }(\mathbf{r})e^{i\frac{{\chi }_0}{2}{\mathbf{k}}_0\mathbf{r}},$$

$$E_h^{ESW}({\mathbf{R}}_{out})=\left(\frac{i\pi }{\lambda }\right){\chi }_h\int S_h(\mathbf{r})f_c(\mathbf{r})e^{i\Delta \mathbf{q}\cdot \mathbf{r}}d{s}_h,$$

$$f_{c}(\mathbf{r})=\left|f_{c}(\mathbf{r})\right|e^{i\varphi_{c}(\mathbf{r})}=\exp\left[i\frac{\chi_{0}}{2}\mathbf{k}_{0}\cdot(\mathbf{r}-\mathbf{R}_{in})+i\frac{\chi_{0}}{2}\mathbf{k}_{h}\cdot(\mathbf{R}_{out}-\mathbf{r})\right]\$$

$$\left|f_c(\mathbf{r})\right|=\exp \left[-\frac{{\chi }_{0i}}{2}{\mathbf{k}}_0\cdot (\mathbf{r}-{\mathbf{R}}_{in})-\frac{{\chi }_{0i}}{2}{\mathbf{k}}_h\cdot ({\mathbf{R}}_{out}-\mathbf{r})\right]$$

$${\phi }_c(\mathbf{r})=\frac{{\chi }_{0r}}{2}{\mathbf{k}}_0\cdot (\mathbf{r}-{\mathbf{R}}_{in})+\frac{{\chi }_{0r}}{2}{\mathbf{k}}_h\cdot ({\mathbf{R}}_{out}-\mathbf{r}).$$

$$\Delta\mathbf{E}(\mathbf{r})-$\mathrm{g}\mathrm{r}\mathrm{a}\mathrm{d}$($\mathrm{d}\mathrm{i}\mathrm{v}$\mathbf{E}(\mathbf{r}))+\frac{\omega^{2}}{c^{2}}[1+\chi(\mathbf{r})]\mathbf{E}(\mathbf{r})=0\ ,$$

$${\chi }^{(id)}(\mathbf{r})=\sum _h{\chi }_h^{(id)}{e}^{i\mathbf{h}\cdot \mathbf{r}},$$

$${\chi }_h(\mathbf{r})={\chi }_h^{(id)}{e}^{-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}.$$

$$\begin{array}{r}\frac{\partial E_0(\mathbf{r})}{\partial s_0}=\frac{i\pi }{\lambda }[{\chi }_0{E}_0(\mathbf{r})+{\chi }_{\bar{h}}{e}^{i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_h(\mathbf{r})],\\ \frac{\partial E_h(\mathbf{r})}{\partial s_h}=\frac{i\pi }{\lambda }[{\chi }_0{E}_h(\mathbf{r})+{\chi }_h{e}^{-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_0(\mathbf{r})].\end{array}$$

$${\chi }^{\prime }(\mathbf{r})=\sum _h[{\chi }_h{e}^{i({\mathbf{h}}^{\prime }-\mathbf{h})\cdot \mathbf{r}}]e^{i\mathbf{h}\cdot \mathbf{r}}.$$

$$\begin{array}{r}{\chi }_h\to {\chi }_h{e}^{i({\mathbf{h}}^{\prime }-\mathbf{h})\cdot \mathbf{r}}\to {\chi }_h{e}^{i\Delta \mathbf{q}\cdot \mathbf{r}}\\ {\chi }_{\bar{h}}\to {\chi }_{\bar{h}}{e}^{i(\mathbf{h}-{\mathbf{h}}^{\prime })\cdot \mathbf{r}}\to {\chi }_{\bar{h}}{e}^{-i\Delta \mathbf{q}\cdot \mathbf{r}}.\end{array}$$

$$\begin{array}{r}{E}_{0,h}(\mathbf{r})=E_{0,h}^{\prime }(\mathbf{r})e^{i\frac{{\chi }_0}{2}{\mathbf{k}}_{0,h}\cdot \mathbf{r}},\end{array}$$

$$\begin{array}{r}\frac{\partial E_0^{\prime }(\mathbf{r})}{\partial s_0}=\frac{i\pi }{\lambda }{\chi }_{\bar{h}}{e}^{-i\Delta {\mathbf{q}}^{\prime }\cdot \mathbf{r}+i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_h^{\prime }(\mathbf{r}),\\ \frac{\partial E_h^{\prime }(\mathbf{r})}{\partial s_h}=\frac{i\pi }{\lambda }{\chi }_h{e}^{i\Delta {\mathbf{q}}^{\prime }\cdot \mathbf{r}-i\mathbf{h}\cdot \mathbf{u}(\mathbf{r})}{E}_0^{\prime }(\mathbf{r}),\end{array}$$

$$\Delta {\mathbf{q}}^{\prime }=\Delta \mathbf{q}+\frac{{\chi }_0}{2}({\text{k}}_0-{\text{k}}_h).$$

$$E_{0,h}^{{}^{\prime }vac}({\mathbf{R}}_b)=E_{0,h}^{{}^{\prime }cryst}({\mathbf{R}}_b)e^{i\frac{{\chi }_0}{2}{\mathbf{k}}_{0,h}\cdot {\mathbf{R}}_b},$$

$$\frac{\partial }{\partial x}f(x-\frac{p}{2},y)=\frac{f(x,y)-f(x-p,y)}{p}.$$

$$f(x-\frac{p}{2},y)=\frac{f(x,y)}{2}+\frac{f(x-p,y)}{2}.$$

Full text

I Introduction

Since its first demonstration Robinson et al. (2001); Vartanyants and Robinson (2001); Williams et al. (2003) Bragg Coherent x-ray Diffractive Imaging (CXDI) has become a powerful technique for analysis of microstructure and strain distribution in submicron crystalline samples Pfeifer et al. (2006); Robinson and Harder (2009); Harder and Robinson (2013); Ulvestad et al. (2015); Robinson et al. (2016); Ulvestad et al. (2016); Dzhigaev et al. (2016). Recently this approach was extended to imaging ultrafast dynamics in nanocrystals using free-electron lasers Clark et al. (2013). Nowadays there are several actively exploited experimental approaches based on the Bragg CXDI concept, among those are Bragg ptychography Godard et al. (2011); Zhu et al. (2015); Dzhigaev et al. (2017) and Fourier transform holography Chamard et al. (2010) (see for review of Bragg CXDI methods Ref. Vartanyants and Yefanov (2015)). In Bragg CXDI technique a finite crystalline sample is illuminated by intense coherent x-ray beam and an interference pattern in the vicinity of a single or several Bragg reflections is recorded Shabalin et al. (2016). An inversion of such data from reciprocal to real space by means of three-dimensional (3D) Fourier transformation provides a high-resolution image of a continuous scattering density distribution in the crystal. The phase of this complex function represents the projection of a local deformation field on the reciprocal lattice vector Vartanyants and Yefanov (2015).

For most of such experiments the dimensions of considered specimens are rather small, therefore the approximation of a single scattering event for hard x-rays is typically used. In the theory of x-ray diffraction by crystals this approach is commonly referred to as the kinematical approximation which is valid while the intensity of the diffracted radiation is small in comparison to the intensity of the incident wave Als-Nielsen and McMorrow (2011). Kinematical description provides a simple expression which allows to calculate scattered amplitude from a finite size crystal as a Fourier transform of its electron density. Such simplification is not applicable for larger crystals, with the sizes bigger than the so-called extinction length Batterman and Cole (1964); Pinsker (1978); Vartanyants and Kovalchuk (2001a); Authier (2003), where effects of cross-coupling between the diffracted and transmitted waves, together with refraction and absorption might become significant and affect Bragg CXDI reconstruction. These effects can be fully described in the frame of the dynamical theory which, however, does not provide a simple analytical expression for the scattering amplitude from a strained crystal of arbitrary shape. This theory has been extensively developed already for decades Batterman and Cole (1964); Pinsker (1978); Vartanyants and Kovalchuk (2001a); Authier (2003), but the influence of the dynamical effects on the results of Bragg CXDI has not yet been fully studied up to now.

The dynamical theory of x-ray diffraction considers the interaction of the wave field with the periodic potential of the crystal lattice taking into account all multiple scattering effects. In this theory one of the most convenient ways to propagate the transmitted and diffracted components of the wave field through the weakly strained crystalline media is based on a set of differential equations with corresponding boundary conditions. This approach developed by Takagi and Taupin Takagi (1962); Taupin (1967); Takagi (1969) describes a general case of the two-beam dynamical diffraction on a perfect or weakly distorted crystal. An analytical solution of these equations is nontrivial and can be performed for a few specific cases only, such as a crystal plate finite in one dimension, but infinite in two other dimensions Afanas’ev and Kohn (1971); Authier (2003); Vartanyants and Kovalchuk (2001a). In a recent work Gorobtsov and Vartanyants (2016) an analytical solution of the Takagi-Taupin equations was found for the phase of the transmitted beam in a quasi-kinematical approximation. Methods of numerical integration of the Takagi-Taupin equations for simulations of the wave field distribution in the crystal were developed in Refs. Epelboin and Riglet (1979); Gronkowski (1991); Mocella et al. (2003). Rapidly increasing number of publications on coherent x-ray scattering experiments on finite size crystals in recent years resulted in a growing interest to understand the role of the dynamical scattering effects in these experiments. For example, Darwin recurrence formalism was applied to study dynamical scattering effects in reciprocal space mapping while scattering on a crystal of rectangular cross-section Kolosov and Punegov (2005); Punegov et al. (2014, 2016). A different approach to solve the Takagi-Taupin equations iteratively via a converging series for a finite size crystal was proposed in Ref. Yan and Li (2014). However, there was still no full analysis of the consequences of the dynamical effects on the reconstruction of the crystal shape and strain field.

In this work we present a general model based on a modification of the Takagi-Taupin equations optimized for geometry of Bragg CXDI measurement aiming to facilitate a numerical solution in a finite 3D crystal of an arbitrary shape in the presence of deformations. Using this model we performed a series of calculations of 3D maps of the complex scattered amplitude distribution in the surrounding of a Bragg reflection for cube- and hemisphere-shaped crystals of different sizes. After kinematical (Fourier) inversion of the simulated 3D reciprocal space data sets into real space the results were compared with their original ones thus revealing a character of the artifacts, introduced by the dynamical diffraction. Next, we analysed effects of refraction and absorption on the reconstructed shape and phase in real space. By neglecting coupling term for the transmitted and diffracted beams we found an analytical solution of the Takagi-Taupin equations that allows to separate the contributions of refraction and absorption. We also determined a correction function that allows to eliminate these effects in the Bragg CXDI. A similar approach was previously presented in Ref. Harder et al. (2007) where the refraction phase shift was calculated accordingly to the optical path for each position in the crystal and subtracted from the results of reconstruction. In this work we provide detailed analytical consideration and demonstrate our method of correction on simulations. We also discuss limitations of the kinematical approach in the Bragg CXDI method.

II Bragg CXDI technique

Typical geometry of a Bragg CXDI experiment assumes an isolated crystal fully illuminated by a coherent x-ray beam. The size of a crystal is generally significantly smaller than the beam and incoming wave field is considered to be a plane wave Vartanyants and Yefanov (2015). The diffracted intensities are recorded by a two-dimensional (2D) pixelated detector located in the far-field and series of diffraction patterns are measured by rotating the sample in the region around the Bragg peak. In kinematical approximation the complex scattered amplitude $A(\textbf{q})$ in the vicinity of a reflection with the corresponding reciprocal lattice vector $\mathbf{h}$ is given by a Fourier integral Vartanyants and Yefanov (2015)

[TABLE]

Here $F_{h}$ is the structure factor, $V_{u.c.}$ is the volume of the unit cell, the momentum transfer vector $\mathbf{q}$ is defined as $\mathbf{q}=\mathbf{Q}-\mathbf{h}$, where $\mathbf{Q}=\mathbf{k}_{h}-\mathbf{k}_{0}$. In kinematical approximation both incoming $\mathbf{k}_{0}$ and diffracted $\mathbf{k}_{h}$ vectors are defined in vacuum and have the magnitude $|\mathbf{k}_{0}|=|\mathbf{k}_{h}|=2\pi/\lambda$, where $\lambda$ is the wavelength of radiation. In Eq. (1) we introduced a complex crystalline function

[TABLE]

where its amplitude $s_{h}(\mathbf{r})$ is so-called shape function, that is defined as unity within the crystal and zero everywhere outside it and its phase $\varphi_{h}(\mathbf{r})$ is proportional to the local deformation field $\mathbf{u}(\mathbf{r})$ that describes displacement of atoms from the ideal lattice positions. In the case of a perfect crystal the intensity distribution function given by a square modulus of the expression (1) is centrosymmetric with respect to the specific reciprocal lattice nodes. However, in the presence of a deformation field this symmetry breaks down thus encoding information about the lattice deformations Vartanyants and Yefanov (2015).

Equation (1) is a basic concept of the Bragg CXDI method. In particular, it directly shows that the reconstructed complex crystalline function has its amplitude $s_{h}(\mathbf{r})$ that is determined by the shape function of the crystalline part of the sample and the phase $\varphi_{h}(\mathbf{r})$. Here we want to point out that the shape function introduced in Eq. (1) does not give any information about electron density of the sample. Such information can be deduced only from the CXDI forward scattering experiments (see for review Ref. Nugent (2010)). Variations of the values of the shape function inside the crystal describe rather modulations of atomic planes associated with the chosen reflection and not electron density modulations. To distinguish between these two cases in the following we will call the shape function also crystalline function. The phase introduced in Eq. (2) by its definition $\varphi_{h}(\mathbf{r})=-\mathbf{h}\cdot\mathbf{u}(\mathbf{r})$ can be attributed to the projection of the local displacement field on the reciprocal lattice vector h around which the measurements are performed. The negative sign of the phase reflects the fact that the positive displacement (expanded lattice) leads to the positional shift of the Bragg peak towards lower momentum transfer values of $\mathbf{Q}$. Taking into account that there is an ambiguity in the constant shift of the phase in the phase retrieval, typically it is the difference between the strained and relaxed parts of the crystal that is determined in CXDI experiment and not its absolute value. It is important to note here that in kinematical approximation described by Eq.  (1) effects of refraction and absorption are not included.

In a typical Bragg CXDI experiment at synchrotron sources 3D measurements of the scattered intensity in the vicinity of the Bragg peak are obtained by an angular scan of the sample with the fixed directions of the incident beam and detector. The concept of such measurement in reciprocal space is depicted in Fig. 1. If Bragg conditions are exactly satisfied the Ewald sphere crosses the selected reciprocal lattice node. At this specific angular position of the crystal the momentum transfer vector $\mathbf{Q}$ coincides with the reciprocal lattice vector $\mathbf{h}$. When the crystal is rotated by an angle $\Delta\theta$, end of the reciprocal lattice vector moves by $\Delta\mathbf{q}=\mathbf{h}^{\prime}-\mathbf{h}$, where $\mathbf{h}^{\prime}$ is the reciprocal lattice vector at the new crystal orientation. Typical values of the angular deviation in Bragg CXDI experiments do not exceed one degree, therefore the length of the vector $\Delta\mathbf{q}$ can be well approximated as $|\mathbf{h}|\Delta\theta$. In our formalism we assume that directions and magnitudes of the incident and diffracted wave vectors $\mathbf{k}_{0}$, $\mathbf{k}_{h}$ are constant during the rocking scan and they always form a constant angle $2\theta_{B}$ at all values of the angular deviation $\Delta\theta$. As such, the wave vector of diffracted field $\mathbf{k}_{h}$ is defined as a constant vector of magnitude $2\pi/\lambda$ pointing at that position on the Ewald sphere which crosses the reciprocal lattice node when the Bragg condition is exactly fulfilled

[TABLE]

III Model description

Conventional Bragg CXDI is based on Eq. (1) which is valid only in the frame of kinematical approximation. For large crystals the kinematical description breaks down and equation (1) cannot be used any more. In this section we will discuss how this simple approach can be modified when dynamical scattering effects are taken into account.

Here dynamical simulations of the Bragg CXDI will be performed in the geometry described in the previous section (see Fig. 1). A detailed sketch of the implemented numerical model is presented in Fig. 2. As shown in this figure, the crystal is embedded in a 3D shape rhombic prism formed by the directions of the incident and diffracted vectors $\mathbf{k}_{0}$ and $\mathbf{k}_{h}$ and centered around the crystal rotation axis. Simulations of the scattered amplitudes for each value of the angular deviation $\Delta\theta$ are performed in four steps. On the first step, the 2D distribution of the incoming wave field $E_{in}(\mathbf{r})$ is projected on the left facet of the rhombic prism (shown as $s_{0}^{in}$ in Fig. 2). Next, the wavefield is propagated through the scattering volume by numerical solution of the Takagi-Taupin equations Takagi (1962); Taupin (1967); Takagi (1969). As a result corresponding 2D distribution of the transmitted $E_{0}(\mathbf{r})$ and diffracted $E_{h}(\mathbf{r})$ amplitudes is obtained at the exit facets of the rhomb (at $s_{0}^{out}$ and $s_{h}^{out}$ facets in Fig. 2). At the third step, thus determined distribution of the diffracted amplitude $E_{h}(\mathbf{r})$ is projected on the plane, perpendicular to the diffracted wave vector $\mathbf{k}_{h}$, yielding the 2D exit surface wave (ESW) $E_{h}^{ESW}(\mathbf{r})$. Finally, on the fourth step, the ESW is propagated to the detector plane, which in the far-field (Fraunhofer) limit is obtained by applying 2D Fourier transform Not (a). When all series of 2D diffraction patterns as a function of the rocking angle $\Delta\theta$ are simulated, we merge them all together into a 3D scattered intensity map in reciprocal space. By this simulation of the intensity distribution in the far-field including dynamical effects in scattering are finalised.

In order to understand what kind of artifacts are introduced by the dynamical scattering the simulated intensity distribution have to be inverted to real space by applying the phase retrieval techniques Fienup (1982); Marchesini (2007). In our case, as soon as the complex amplitudes in the far-field are known, they can be directly inverted to real space without the need of the phase retrieval procedure. Characterization of the dynamical artifacts in real space is performed by a comparison of the output of these simulations with the original crystalline function $S_{h}(\mathbf{r})$ (see Eq. (2)). In the next sections we describe all these steps in more details.

III.1 Propagation of the wave field through the crystal

In our simulations of Bragg CXDI we used the laboratory coordinate system in which the direction of the incident beam and detector position are fixed during the angular scan and the sample is rotating. The origin of the coordinate system is chosen on the crystal rotation axis that is parallel to the y-axis (see Fig. 2). In the following we assume the two-beam diffraction conditions

[TABLE]

where $\mathbf{e}_{0s}$ and $\mathbf{e}_{hs}$ are the polarization unit vectors and $s$ is the polarization index. As a sample we consider a perfect or weakly deformed finite size crystal. To propagate the complex electric field through the volume of the three-dimensional crystal we introduce the symmetric form of the Takagi-Taupin equations (see Appendix A for details)

[TABLE]

which are supplemented by the boundary conditions. Here the partial derivatives $\partial/\partial s_{0},\partial/\partial s_{h}$ are taken along the directions of the wave vectors $\mathbf{k}_{0}$ and $\mathbf{k}_{h}$. Response of the crystal is described by the Fourier components of the susceptibility $\chi_{0}=\chi_{0r}+i\chi_{0i}$ and $\chi_{h,\bar{h}}=\chi_{hr,\bar{h}r}+i\chi_{hi,\bar{h}i}$. Real and imaginary parts of the zeroth component of the susceptibility $\chi_{0r}$,$\chi_{0i}$ describe effects of refraction and absorption, respectively. The term with $\chi_{h}$ describes diffraction of the transmitted component $E_{0s}(\mathbf{r})$ by a set of crystallographic planes with the reciprocal vector h. In its turn, the diffracted component $E_{hs}(\mathbf{r})$ undergoes diffraction by the same set of planes but from the opposite side, which is described by the term with $\chi_{\bar{h}}$. Vector $\Delta\mathbf{q}$ determines the angular deviation from the exact Bragg condition as shown in Fig. 1. In equations (5) $C$ stands for the polarization factor, which is equal to unity in the case of a $\sigma-$polarization and $\cos 2\theta_{B}$ in the case of a $\pi-$polarization. Without restricting the generality in the following we will assume only $\sigma-$polarization with $C=1$ and omit polarization index $s$ in the wavefield amplitudes.

The boundary conditions assume that the total electric field, represented by equation (4) is continuous everywhere on the crystal-vacuum boundary. In our formalism, similar to Refs Takagi (1962, 1969); Taupin (1967), we assume that the wave vectors $\mathbf{k}_{0,h}$ are the same inside and outside the crystal, therefore the amplitudes $E_{0,h}(\mathbf{r})$ are continuous functions on the boundary of the crystal. Consequently, equations (5) do not require any specific transformation of the amplitudes $E_{0,h}(\mathbf{r})$ on the crystal-vacuum boundary which is particularly convenient in the case of a three-dimensional crystal with an arbitrary shape. We note that the Fourier components $\chi_{0,h,\bar{h}}$ of susceptibility drop down to zero outside the crystal and thereby undergo discontinuity on the crystal boundary.

For a finite size crystal the evolution of the wave field depends on the crystals size, shape, and diffraction geometry. When exact Bragg conditions are satisfied, the transfer of energy from the transmitted beam into the diffracted beam is strongly enhanced due to constructive interference of the wavefield inside the crystal. At the same time at these conditions the wave field is not penetrating deep into crystal. This effect, known as extinction Pinsker (1978); Authier (2003), is described by the characteristic decay length of the wave field commonly referred to as the extinction length $L_{ex}$ Vartanyants and Kovalchuk (2001b); Not (b)

[TABLE]

where $\gamma_{0,h}=\cos(\mathbf{n}\cdot\mathbf{k}_{0,h})$ are the direction cosines and $\mathbf{n}$ is the inward normal to the entrance surface of the crystal. In the denominator of the expression (6) the real part of the complex valued square root $\sqrt{\chi_{h}\chi_{\bar{h}}}$ is used. The extinction length (6) is commonly referred to as a characteristic value to distinguish between the cases of the kinematical and dynamical diffraction. When the crystal size is much smaller than the extinction length the effects of coupling between the transmitted and diffracted components of the wave field are small and kinematical approximation can be used safely. When the crystal size is about or bigger than the extinction length these effects are becoming important and the dynamical theory has to be used.

III.2 Numerical solution of the Takagi-Taupin equations

In this work we perform numerical integration of the Takagi-Taupin equations (5) applying an approach similar to that described in Ref. Gronkowski (1991). To propagate the wavefield along the directions of partial derivatives $\partial/\partial s_{0},\partial/\partial s_{h}$, we introduce the laboratory coordinate system with the origin on the crystal rotation axis. The set of basis vectors $\{\mathbf{s}_{0},\mathbf{s}_{h},\mathbf{s}_{y}\}$ is represented by the unit vectors in the direction of the incident beam ($\mathbf{s}_{0}$), diffracted beam ($\mathbf{s}_{h}$) and normal to the scattering plane ($\mathbf{s}_{y}$) (see Fig. 2). Thus, the partial derivatives are taken along $\mathbf{s}_{0}$ and $\mathbf{s}_{h}$ vectors, and rotation is performed around the $\mathbf{s}_{y}$-axis. The angle between the vectors $\mathbf{s}_{0}$ and $\mathbf{s}_{h}$ is equal to $$2$\theta_{B}$, therefore the coordinate system generally is not orthogonal. Any position within the considered volume can be described by the radius vector $\mathbf{r}=s_{0}\mathbf{s}_{0}+s_{h}\mathbf{s}_{h}+s_{y}\mathbf{s}_{y}$, where $s_{0,h,y}$ are corresponding coordinates.

We perform the numerical integration over a rhombic prism in which the whole crystal is embedded, as shown in Fig. 2. More specifically, the prism is sliced to a set of layers, defined for different values of $s_{y}$ coordinate parallel to the scattering plane and Takagi-Taupin equations (5) are solved in the two-dimensional grid independently for each of these layers. Since directions $\mathbf{s}_{0,h}$ do not depend on the angular deviation $\Delta\theta$ in chosen coordinate system the whole grid remains invariable during the angular scan, while rotation transformations are applied to the susceptibility and shape function of the crystal (see Appendix A for details). The nodes, which belong to the crystal, are characterized by the values of Fourier components of the susceptibility, which are replaced by zeros for the nodes outside the crystal.

In the numerical integration method the complex amplitudes $E_{0,h}(\mathbf{r})$ are represented by a discrete set of values over all integration grid and the Takagi-Taupin equations (5) are transformed into a pair of recurrence relations (see Appendix C for details). The inset in Fig. 2 shows the recurrence property of the obtained equations for the neighboring nodes of the integration grid. For the node $(i,j)$ the values of amplitudes $E_{0,h}^{(i,j)}$ are calculated from the values $E_{0,h}^{(i-1,j)}$ and $E_{0,h}^{(i,j-1)}$ at the previous nodes $(i-1,j)$ and $(i,j-1)$. In such a way calculations proceed from node to node in the direction of the transmitted and diffracted beams. The values of the amplitudes $E_{0,h}(\mathbf{r})$ on the left ($s_{0}^{in}$) and bottom ($s_{h}^{in}$) sides of the prism are defined as

[TABLE]

Such form of boundary conditions is universal and particularly convenient to implement for numerical integration in the case of a three-dimensional crystal with an arbitrary shape. Once established these boundary conditions can be applied to any shape and orientation of the crystal embedded into the integration prism.

III.3 Propagation to the detector plane

Numerical integration of the Takagi-Taupin equations (5) over the rhombic prism results in the complex amplitude of the transmitted beam $E_{0}(\mathbf{r})$ at the right facet of the prism ($s_{0}^{out}$) and diffracted wave $E_{h}(\mathbf{r})$ at the upper facet of the prism ($s_{h}^{out}$). For the free space propagation to the detector we exploit an orthogonal coordinate system with the basis $\{\mathbf{s}_{\bot},\mathbf{s}_{h},\mathbf{s}_{y}\}$, where the vector $\mathbf{s}_{\bot}=\mathbf{s}_{h}\times\mathbf{s}_{y}$ is introduced (see Fig. 2). This vector is perpendicular to the direction of propagation for diffracted component and lies in the scattering plane, therefore the transition for the coordinate $s_{0}$ is performed by means of a simple projection $s_{\bot}=s_{0}\sin($2$\theta_{B})$. The result of such projection, applied to the calculated 2D distribution of the diffracted wave field will be further referred to as the exit surface wave $E_{h}^{ESW}(s_{\bot},s_{y},\Delta\mathbf{q})$. To determine the scattered amplitude at the detector plane in the far-field we apply 2D Fourier transform to the exit surface wave

[TABLE]

where $q_{\bot},q_{y}$ are corresponding coordinates in reciprocal space.

In generic Bragg CXDI experiment the measured diffraction pattern corresponds to the cut of reciprocal space by the Ewald sphere (see Fig. 1). Our model does not account for divergence of the wave field $E_{0,h}(\mathbf{r})$ while its propagation in a crystal (see Appendix A). In fact, this is similar to the projection approximation, when the simulated 2D diffraction pattern is attributed to the flat surface in reciprocal space (see Fig. 1). As such, the 2D distribution of the scattered amplitude $A(q_{\bot},q_{y},\Delta\mathbf{q})$ defined by Eq. (8) is determined in a plane in reciprocal space perpendicular to the direction of the diffracted wave $\mathbf{s_{h}}$ and corresponding to a fixed angular deviation $\Delta\mathbf{q}$. Changing the value of the rocking angle $\Delta\theta$ the full set of complex amplitudes $A(q_{\bot},q_{y},\Delta\mathbf{q})$ is determined in reciprocal space. By taking the square modulus of the amplitudes the final 3D distribution of the intensity $I(q_{\bot},q_{y},\Delta\mathbf{q})=\left|A(q_{\bot},q_{y},\Delta\mathbf{q})\right|^{2}$ is obtained. As a next step, this set of 2D images is interpolated on a 3D uniform grid with the orthogonal coordinates $q_{x},q_{y},q_{z}$ (see Fig. 3(a)).

As soon as phases of the scattered amplitudes are known in our simulations we perform 3D inverse Fourier transform of simulated 3D amplitudes

[TABLE]

where $S_{h}(\mathbf{r})$ is a complex crystalline function defined in (2). In kinematical approximation it should reproduce the crystal shape by its amplitude and be proportional to the projected strain field by its phase (compare with Eq. (2)). According to this approach, if a crystal is unstrained, the inversion of the scattered amplitude by Eq. (9) should give a real shape function with the constant amplitude values. We will see in the following how dynamical scattering may affect these results.

IV Results

In order to illustrate general features of the dynamical scattering effects in the Bragg CXDI we considered first a simple object in the form of a cubic-shaped gold crystal without strain. A schematics of the diffraction geometry in real and reciprocal space and the orthogonal coordinate system with the $x,y,z$-axes oriented along the cube edges is shown in Fig. 3. We assume that a cubic unit cell (with a lattice parameter a=$$4.078\text{\,}\mathrm{\SIUnitSymbolAngstrom}) is also aligned along the same coordinate axes. In our simulations we considered the incident plane wave with $8\text{\,}\mathrm{keV}$ photon energy (wavelength \lambda=$$1.55\text{\,}\mathrm{\SIUnitSymbolAngstrom}) and 004 reflection conditions. In this scattering geometry reciprocal space vector $\mathbf{h}_{004}$ is parallel to $q_{z}-$axis in reciprocal space (see Fig. 3(a)) and scattering plane is parallel to $xz-$plane in real space (see Fig. 3(b)). The Bragg angle in these conditions is $\theta_{B}=49.47^{\circ}$ and values of the extinction length $L_{ex}$ (6) are $711$ nm and $607$ nm in the Bragg and Laue geometry, respectively.

We performed simulations for two crystal sizes of $100\text{\,}\mathrm{nm}$ and $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$. For $100\text{\,}\mathrm{nm}$ crystal the angular scan was performed covering the angular range from $-3.3^{\circ}$ to $+3.3^{\circ}$, with the angular increment of $6\cdot 10^{-3}$ degree. The exit surface wave field distribution $E_{h}^{ESW}(s_{\bot},\Delta\mathbf{q}=0)$ obtained by a solution of Takagi-Taupin equations at exact Bragg conditions is presented in Fig. 4. The amplitude of the exit wave calculated in the frame of the dynamical theory (red curve) is compared to the results of the kinematical theory (black curve) obtained by setting $\chi_{0,\bar{h}}=0$ (see AppendixD). Nearly complete coincidence in simulations for a $100\text{\,}\mathrm{nm}$ crystal (see Fig. 4(a)) suggests that the cross coupling between the diffracted and transmitted waves is not strong enough to have any significant effect on the scattering and, therefore, kinematical approximation provides a rather accurate result. The calculated phase profile (green curve) shows a small phase shift, which can be attributed to refraction. The phase distribution is shown relative to the phase of the incoming wave that was set to zero, so that the phase of the diffracted wave at the top left corner of the crystal (see Fig. 3) appears to be zero as well. As a result of our simulations we can see that the phase due to refraction accumulates more for the waves propagating the longest distance in the crystal from its depth and finally reaches its minimum value of $-0.25\text{\,}\mathrm{rad}$.

Similar simulations performed for a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au crystal are presented in Fig. 4(b). For a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ crystal the angular scan was performed covering the angular range from $-0.83^{\circ}$ to $+0.83^{\circ}$, with the angular increment of $3.3\cdot 10^{-3}$ degree. The dynamical calculations revealed a considerably lower amplitude profile in comparison to the kinematical prediction, which can be attributed to the attenuation of the transmitted wave due to extinction. This affects mostly the lower part of the crystal. We would like to note that contribution to attenuation due to normal absorption is much lower than extinction effect. Indeed, taking into account that normal absorption length for gold at the considered photon energy is $2.9\text{\,}\mathrm{\SIUnitSymbolMicro m}$, we obtain attenuation of the x-ray amplitude only by 16% on the length of a Au particle of $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in size. The phase distribution, in fact, reproduces major features of the phase for $100\text{\,}\mathrm{nm}$ crystal, which supports our observation that characteristic phase gradient originate mostly from refraction. At the same time, the phase profile at the lower part of the cube reveals slight but noticeable bending, which cannot be attributed to refraction, since the refraction phase is linear. We also notice that due to a bigger crystal size the observed phase shift is about one order of magnitude larger than in the case of $100\text{\,}\mathrm{nm}$ crystal and reaches the value of $-2.17$ rad. Below we will analyse results of inversion obtained for the two different crystal sizes separately.

IV.1 100 nm Au crystal of a cubic shape

Results of inversion of the whole 3D reciprocal space dataset for a $100\text{\,}\mathrm{nm}$ Au crystal obtained by the dynamical simulations are presented in Fig. 5. A 2D distribution of the crystalline amplitude function $s_{h}(\mathbf{r})$ in $xz$-slice taken through the center ($y$=0) of the crystal is shown in Fig. 5(a); the line profiles along $x,y,z$-axes are given in Fig. 5(b). The phase distribution $\varphi_{h}(\mathbf{r})$ is presented in a similar way by the corresponding slice in Fig. 5(c) and line profiles in Fig. 5(d), respectively. Outside of the cube the amplitude of the reconstructed complex density function is rapidly going down (see Fig. 5(a,b)). In this region the phases are not defined, therefore, the phase distribution presented in Fig. 5(c,d) was cropped by the cube edges.

Distribution of the crystalline function reveals, as expected, well defined cubic structure of our model sample. We should note here that due to plane facets of a cubic sample crystal truncation rods Als-Nielsen and McMorrow (2011) are extending quite far in reciprocal space and induce observed oscillations in the crystalline amplitude function obtained by Fourier inversion. Therefore, slight periodic variations of its values are due to truncation of reciprocal space intensities imposed by the limited range where simulations were performed.

The most intriguing and not expected result was obtained for the phase $\varphi_{h}(\mathbf{r})$ of the crystalline function (see Fig. 5(b)). Instead of being a uniform function inside of an unstrained crystal it shows slight variation of the phase going down to the values of about $-0.3$ rad. When attributed to strain, these values of the phase would give rise to the displacement of about $0.049$ Å and associated strain for a Au crystalline sample of $~{}1.2\cdot 10^{-2}$. As we will show in the following these variations of the phase can be attributed to refraction effects, that are not considered in the conventional kinematical theory. Indeed, on the top facet neither incident nor scattered wave experience refraction, therefore the phase shift is zero. When radiation penetrates deep in the crystal the phase shift due to refraction is accumulated on its way in and out of the crystal. Since refraction index for x-rays is less then one the accumulated phase is negative. For the lower facet of the cube the phase shift reaches its minimum value of about $-0.3$ rad. This value corresponds to an optical path length of x-ray beam going in and out of the crystal, which for the Bragg angles smaller than $63.43^{\circ}$ gives for the phase shift due to refraction $\varphi_{refr}=-(2\pi/\lambda)\delta(d/\cos\theta_{B})\simeq-0.29$ rad, where $\delta$ is the real part of the refraction index. We will introduce later a correction function that will compensate these effects completely and will allow to determine correct values of the phase that can be attributed to strain.

From these simulations we can see that even in the case of very small crystalline samples when dynamical effects should not play any role refraction effects introduce certain phase variations in the reconstructed crystalline function that could lead to a wrong statements about the strain field in the sample.

IV.2 1 $\mu$m Au crystal of a cubic shape

As a next step, we performed simulations for a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au crystal of a cubic shape. Simulated 2D distribution of the modulus of the scattered amplitude $\left|A(q_{x},q_{z})\right|$ taken through the center of reciprocal space and obtained by using the kinematical and dynamical approaches is presented in Fig. 6 (a,b). Two sets of crystal truncation rods perpendicular to the direction of the facets of the crystal as well as a regular structure of the square speckles due to coherent scattering on a cubic shape crystal are well seen in this Figure. At the same time we see a significant difference between simulations performed with the kinematical and dynamical approaches. The later ones show lower contrast and noticeable aberrations in the position and magnitude of the fringes. We also observed an additional intensity in the form of a diagonal cross in the case of the dynamical theory simulations (see Fig. 6(b)) that was also noticed in simulations performed in Ref. Punegov et al. (2016).

The difference in the position and intensity of the speckles is clearly seen in a linear scan of the amplitude $\left|A(q_{z})\right|$ taken along the central rod (see Fig. 6(c)). A comparison of the kinematical (black line) and dynamical (red line) results show a displacement of the whole profile and particularly Bragg peak position in the positive direction of $q_{z}$-axis for the case of the dynamical theory simulations. This result is well known in the dynamical theory Authier (2003) and is due to refraction effect. According to the dynamical theory the angular position of the maximum of the reflectivity curve is shifted from the exact Bragg position to positive values by

[TABLE]

where parameter $\beta=\gamma_{0}/|\gamma_{h}|$ for Bragg and $\beta=\gamma_{0}/\gamma_{h}$ for Laue geometries and the upper sign corresponds to Bragg diffraction and the lower one to Laue diffraction. For Au(004), $8\text{\,}\mathrm{keV}$ and symmetric Bragg geometry ($\gamma_{0}=|\gamma_{h}|$) the equation (10) provides $19.6\text{\,}\mathrm{\SIUnitSymbolArcsecond}$ angular shift which is equivalent to $$5.0\text{,}\mathrm{\SIUnitSymbolMicro m}$^{-1}$ of the positional displacement of the Bragg peak in reciprocal space along $q_{z}$ axis (compare with the similar results obtained in Ref. Punegov et al. (2016)). At the same time, in symmetric Laue geometry ($\gamma_{0}=\gamma_{h}$) no positional shift of the reflectivity curve is observed. In the considered case of a cubic crystal the diffraction geometry is represented by a mixture of symmetric Bragg and Laue cases, therefore, refraction effects characteristic for these two geometries are superimposed. That is revealed in a smearing of the central speckle in the direction of $q_{z}$-axis together with a positional shift of the maximum by $3.5\text{\,}{\mathrm{\SIUnitSymbolMicro}\mathrm{m}}^{-1}$ (see Fig. 6(c)). Although small angular displacement of the whole diffraction pattern due to refraction can be precisely determined in simulations it is rather challenging to consider it experimentally. In most of experiments these effects are neglected and the maximum of the Bragg peak is assumed to be at an exact position of the reciprocal lattice node and is used as a reference position.

It is also well seen in Fig.  6(c) that due to the dynamical scattering contrast of the diffraction pattern is significantly reduced. In experiment this might be erroneously attributed to lack of the transverse coherence and consequently partial coherence illumination Vartanyants and Robinson (2001); Williams et al. (2007), or vibrations of the sample stage. While these effects may be compensated in reconstruction by the multimode decomposition Clark et al. (2012); Thibault and Menzel (2013) and attributed to the incoming field, however their physical origin is quite different and is due to the dynamical scattering effects.

In Fig. 6(d) the corresponding $q_{z}$-profiles of the phase distributions for the kinematical (black line) and dynamical (red line) calculations are presented. Similar to the amplitude profiles shown in Fig. 6(c) a comparison between the kinematical and dynamical results shows a positional displacement of the phase profile in the positive direction of $q_{z}$-axis in the case of the dynamical theory simulations. In addition, the symmetry with respect to the positive and negative directions is broken and more complex structure of the profile is observed.

The most intriguing were results of inversion performed for a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ size crystal. In contrast to the previous case of a small crystal, results of inversion for a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ crystal (see Fig. 7) clearly show visible artifacts in the crystalline amplitude $s_{h}(\mathbf{r})$ (a,b) and phase $\varphi_{h}(\mathbf{r})$ (c,d) distribution in real space. One strong effect, well visible in Fig. 7(a,b), is depletion of the crystalline amplitude towards the bottom of the crystal. This is an expected effect of the dynamical theory. Due to coupling of the incoming and diffracted waves at Bragg conditions the wave field is not propagating inside the crystal. We observed that in our case the values of the amplitude dropped by more than 50% (see Fig. 7(b)) instead of being uniform and constant on the level of one as in the case of a small crystal (see Fig. 5(b)). Another unexpected effect was appearance of an additional intensity, which extends below the bottom of the crystal (see Fig. 7(a,b)). In our simulations for larger crystals (not shown) we observed that this artifact becomes stronger with the increase of the ratio of the crystal size to extinction length. If the crystal shape is unknown before CXDI experiment such dynamical effects can result in a wrong reconstruction of the crystal shape as well as the values of the crystalline amplitude function.

Importantly, our simulations have revealed that the phase profile inside a crystal has a complicated distribution (see Fig. 7(c,d)). We want to remind that initially we were considering Au crystals without any deformation. At the same time, we obtained strong variations in the phase of the inverted crystalline function that should not be interpreted as originating from the crystal lattice deformation. We will see in the following that some features of this phase distribution can be compensated by taking into account refraction effects. Without such corrections the values of the strain field obtained from the CXDI reconstruction could be significantly different from the ones in the sample under investigation and in this way could bring to a wrong interpretation of the results in the Bragg CXDI experiment.

V Treatment of refraction and absorption

Here we will analyse how effects of refraction and absorption could be taken into account. We will perform analysis in the semi-kinematical approximation, when coupling between the incident and diffracted waves could be neglected but refraction and absorption effects will be taken into account specifically (see also Ref. Gorobtsov and Vartanyants (2016)). We want to point out here that in conventional kinematical theory the incident and diffracted waves have no attenuation due to absorption and refraction effects are also neglected.

To take all this into account, we will consider Takagi-Taupin equations (5) in which the coupling term in the first equation, proportional to $\chi_{\bar{h}}$, is eliminated that leads to the following system of equations

[TABLE]

The first equation can be easily solved as

[TABLE]

where $E_{0}^{in}$ is the incoming wavefield that will be put to unity in the following and $\mathbf{R}_{in}=\mathbf{R}_{in}(\mathbf{r})$ is a radius vector of the point, where the incoming beam enters the crystal for the considered element of volume (see sketch in Fig. 8(a)). Substituting this result in the second equation of the system of equations (11) we obtain

[TABLE]

This equation for the diffracted wave can be solved by the following substitution

[TABLE]

which leads finally to the following analytical expression for the exit surface wave

[TABLE]

where as before $S_{h}(\mathbf{r})$ is a complex crystalline function and $\mathbf{R}_{out}=\mathbf{R}_{out}(\mathbf{r})$ is a radius vector of the position on a crystal surface where the diffracted beam exits the crystal for the considered element of volume (see sketch in Fig. 8(a)). We also introduced here a correction function

[TABLE]

with its modulus due to absorption

[TABLE]

and the phase due to refraction

[TABLE]

We want to point out again that in this treatment dynamical scattering effects are completely neglected. Purely kinematical scattering can be directly obtained from Eqs. (15-16) by putting $\chi_{0}=0$ and consequently assuming that $f_{c}(\mathbf{r})\equiv 1$.

As it follows from equations (15-16) the inversion of the reciprocal space dataset to real space should result in a complex function, which is represented by $S_{h}(\mathbf{r})f_{c}(\mathbf{r})$. Therefore, in principle, to determine correctly the shape and strain field in a crystalline particle the correction function (16) should be applied after inversion from reciprocal space. In this correction function an optical path along the incident $\mathbf{k}_{0}\cdot(\mathbf{r}-\mathbf{R}_{in})$ and diffracted $\mathbf{k}_{h}\cdot(\mathbf{R}_{out}-\mathbf{r})$ beams should be calculated for each position $\mathbf{r}$ in a crystal (see sketch in Fig. 8(a)). An estimate is performed for a fixed angular position neglecting small variations of the optical path while the rocking scan. Since the crystal shape and directions of the vectors $\mathbf{k}_{0}$ and $\mathbf{k}_{h}$ are known, the correction function (16) can be evaluated numerically in most of the cases.

Results of such correction applied to the complex electron density distribution obtained for a $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Au crystalline particle are shown in Fig. 8. The correction was performed only for positions inside a cubic volume, leaving the exterior part below the cube unchanged. Comparison of the amplitudes reveals no significant changes in the upper part of the crystal. In the lower part, where effects of absorption due to extinction are stronger, the values of the amplitude are increased from $0.38$ to $0.52$ (compare Fig. 7(a,b) and Fig. 8(a,b)). The correction revealed also a noticeable bump on the $z$-profile of the amplitude (see Fig. 8(b)), which was barely pronounced in Fig. 7(b). Still, major artifacts in the amplitude distribution, such as the depletion of the amplitude of the crystalline function in the bottom part of the crystal remained almost unchanged. By that we conclude that remaining artifacts in the amplitude are related to purely dynamical effects in scattering.

At the same time, by applying correction function in the phase distribution we observed that a strong gradient of phase present in Fig. 7(c,d) is effectively removed (see Fig. 8(c,d)). Small residual aberrations in the range from $0\text{\,}\mathrm{rad}0.3\text{\,}\mathrm{rad}$ are apparently connected to the dynamical effects Gorobtsov and Vartanyants (2016). More specifically, we determined that they can be attributed to the imaginary part of the Fourier components of the susceptibility $\chi_{h}$ and $\chi_{\bar{h}}$, which introduce a small phase shift when the wave is reflected by a crystalline plane. To illustrate this we performed simulations for the same Au crystalline particle in which imaginary parts of the susceptibilities were eliminated from the Takagi-Taupin equations (5) by setting $\chi_{hi}=\chi_{\bar{h}i}=$0$$. Results of these simulations after inversion to real space and applying correction by the function $f_{c}(\mathbf{r})$ (16) are shown in Fig. 9. As we can see from this figure all residual artifacts in the phase distribution were completely removed which approves our conclusion about the origin of these features.

To be sure that our correction function takes into account entire contribution of refraction and absorption we performed complementary simulations (not shown) where the corresponding terms were completely eliminated by setting $\chi_{0}=$0$$ in Takagi-Taupin equations (5). A comparison with the simulations performed by the fully dynamical case and corrections applied by the function $f_{c}(\mathbf{r})$ (16), presented in Fig. 8, showed that both results entirely coincide with each other. This can be explained by a suggestion that contributions due to the dynamical scattering effects are completely decoupled from contributions originating from refraction and absorption.

VI Simulations for a $Pb$ particle of a hemispherical shape

Results of simulations for a perfect cubic Au crystalline particle have shown that the dynamical diffraction can lead to an appearance of artifacts in the real space reconstruction. In order to estimate the contribution of the dynamical effects for a practical case, we considered experimental parameters described in Refs Pfeifer et al. (2006); Harder et al. (2007). In that experiment 3D reconstruction of the Bragg CXDI data was used to characterize the strain distribution in a lead nanocrystal of a hemispherical shape of $0.75\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in diameter. The crystal was coherently illuminated by a monochromatic x-ray beam of $1.38\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ wavelength and Pb (111) reflection was selected. In these experimental conditions the Bragg angle was $13.97\text{\,}\mathrm{\SIUnitSymbolDegree}$ and values of the extinction length $L_{ex}$ were $0.28\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and $1.14\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for the Bragg and Laue geometries, respectively.

For our simulations we considered a shape function represented by a sphere truncated from one side by $1\text{/}3$ of its diameter, as the closest model. Following the description of the experiment, we oriented the truncation plane to form an angle of $27\text{\,}\mathrm{\SIUnitSymbolDegree}$ with respect to the (111) crystallographic plane. The diffraction geometry from two perspective views is schematically shown in Fig. 10. Note orientation of the coordinate axes: similar to the case of simulations for a cubic Au particle the $x$- and $z$-axes lie in the scattering plane and the $y$-axis is orthogonal to them. According to the chosen geometry (see Fig. 10) the cut along $z$-axis is not symmetric with respect to the center but covers the range from -d/2+d/3=$$-125\text{\,}\mathrm{nm} to d/2=$$375\text{\,}\mathrm{nm}.

Results of simulations performed by the dynamical theory are presented in Fig. 11. A series of diffraction patterns were calculated in the angular range from $-0.83\text{\,}\mathrm{\SIUnitSymbolDegree}$ to $0.83\text{\,}\mathrm{\SIUnitSymbolDegree}$ with the $3.3\cdot 10^{-3}$ degrees angular increment. They were merged into a 3D reciprocal space dataset, and then inverted to real space. The amplitude distribution shown in Fig. 11(a, b) reveals slight depletion of the crystalline amplitude function $s_{h}(\mathbf{r})$ in the central part, which corresponds to attenuation of the incident and diffracted waves in the bulk of the crystal. However, this artifact appears to be relatively small (about $10\text{\,}\mathrm{\char 37\relax}$ of the average value) which supports applicability of the kinematical approach in this case. At the same time, in the phase distribution (see Fig. 11(c,d)) a considerable phase gradient with the maximum deviation of the phase about $0.7\text{\,}\mathrm{rad}$ is observed.

To reveal the origin of this phase gradient we applied the correction function $f_{c}(\mathbf{r})$ (16) to the complex crystalline function $S_{h}(\mathbf{r})$ shown in Fig.  11. Results of this correction are presented in Fig. 12. We do not observe any significant changes in the amplitude distribution (see Fig. 12(a, b)), as soon as contribution due to absorption is comparably small for this particle. At the same time, correction due to refraction removed a major part of the gradient in the phase distribution (see Fig. 12(c, d)). Leftover residual variations were on the level of $0.03\text{\,}\mathrm{rad}$ and can be neglected. These values are much less than the values of the maximum phase deviation, which were estimated in Ref. Harder et al. (2007) to be about $1.15\text{\,}\mathrm{rad}$ after correction for refraction effects.

Our results demonstrate, first, that corrections of the phase due to refraction are important even in the case of kinematical scattering and can not be neglected. Second, our theoretical results demonstrate that an approach proposed in Ref.  Harder et al. (2007) can be safely applied in the case when dynamical scattering effects are negligible and kinematical scattering approximation can be used.

VII Conclusions

We present a general model based on a specific form of the Takagi-Taupin equations optimized for geometry of the Bragg CXDI measurement and with the aim to facilitate a numerical solution in a finite 3D crystal of an arbitrary shape in the presence of deformations. As a result, the complex amplitude distributions of the transmitted and diffracted waves on the exit surface are calculated. Propagation to the far-field provides the amplitude and phase distributions of the diffraction pattern that corresponds to a specific cross-section in reciprocal space. By performing a series of such calculations for different values of rotation angle a full 3D reciprocal space dataset in the vicinity of the corresponding reciprocal lattice node can be constructed. The complex crystalline function of the object in real space is obtained by the inverse Fourier transform.

Using this model we performed simulations of the dynamical diffraction on a perfect crystal of gold of a cubic shape of $100\text{\,}\mathrm{nm}$ and $1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ in size. For a small crystal results of our calculations were in full agreement with the kinematical theory. In the simulations for the large crystal artifacts introduced by the dynamical scattering effects were observed in real as well as in reciprocal spaces. We analyzed the contributions of different phenomena, such as refraction, absorption and cross-coupling between the diffracted and transmitted waves. Based on the analytical derivations we developed an approach which corrects the results of reconstructions for the effects of refraction and absorption. Such correction, applied to the results of the simulations, demonstrates a complete removal of corresponding contributions in the real space reconstruction. The residual artifacts in the amplitude and phase distributions are attributed to the dynamical effects of scattering in the crystal. Additional simulation for a practical case of a Bragg CXDI experiment with a hemispherical Pb particle of $750\text{\,}\mathrm{nm}$ in size was also performed. By applying the correction for refraction and absorption we demonstrate that remaining dynamical artifacts were small and did not affect results of the reconstruction.

We conclude that limitations of kinematical approach in the Bragg CXDI experiments depend on the relative values of the crystal size $d$ and extinction length $L_{ex}$. We suggest the following critera. If scattering conditions (crystal shape and orientation) are predominantly Bragg (Laue) than the size of the crystal $d$ should be compared with the corresponding Bragg (Laue) extinction length. When the crystal size is smaller than the corresponding extinction length, dynamical effects should not affect reconstruction significantly. However, even in this case effects of absorption and especially refraction should be specially analyzed. If necessary, correction function should be applied to determine correct values of strain in the sample. In other cases the dynamical theory should be applied.

Finally, we think that our findings will be of high importance for all groups working in the fast developing field of coherent scattering and imaging in Bragg scattering conditions.

Acknowledgements.

Authors are grateful to E. Weckert for fruitful discussions and support of the project. We acknowledge careful reading of the manuscript and helpful discussions with O. Gorobtsov and I. Sergeev. This research was partially supported by the Virtual Institute VH-VI-403 of the Helmholtz Association.

Appendix A Derivation of modified Takagi-Taupin equations

Here we describe derivation of the Takagi-Taupin equations in symmetric form. The properties of the electric field vector $\mathbf{E}(\mathbf{r})$ inside a crystal are described by the following wave propagation equation

[TABLE]

where $\omega$ is the frequency of the wave field and $c$ is the speed of light. In equation (19) $\chi(\mathbf{r})$ is the susceptibility of the crystal. We assume in the following the case of $\sigma-$polarization, so the electric field will be further considered as a scalar field.

In the case of a perfect crystal the susceptibility $\chi(\mathbf{r})$ is a periodic function with the period of the crystal lattice that can be expanded as a Fourier series

[TABLE]

where $\mathbf{h}$ is the reciprocal lattice vector. In equation (20) the summation is carried out over all reciprocal lattice vectors. In the case of weak deformations, when relative displacements are small the susceptibility of the crystal $\chi(\mathbf{r})$ is defined from that of a perfect one according to the relation Takagi (1962), $\chi(\mathbf{r})=\chi^{(id)}(\mathbf{r}-\mathbf{u}(\mathbf{r}))$. The Fourier components of the susceptibility in the weakly deformed crystal now depend on the coordinate $\mathbf{r}$ and can be defined as

[TABLE]

The solution of equation (19) may be found in the form of an expansion analogous to Bloch waves (see equation (4)) and leads to a well known form of the Takagi-Taupin equations (see for example Authier (2003); Vartanyants and Kovalchuk (2001a)).

We will consider now that the orientation of the crystal satisfies the exact Bragg conditions. In this case we have for the Takagi-Taupin equations

[TABLE]

We assume now that the crystal is rotated by an angle $\Delta\theta$, and denote reciprocal lattice vector at this new orientation as $\mathbf{h}^{\prime}=\mathbf{h}+\Delta\mathbf{q}$. Then, Fourier decomposition of the susceptibility $\chi^{\prime}(\mathbf{r})$ for this new angular position of a crystal may be written as

[TABLE]

Comparing this expression with the decomposition (20) we conclude that in the equations (22) the following substitutions should be made

[TABLE]

Following this approach, the Takagi-Taupin equations for this new crystal orientation can be finally written in the form (5). The angular dependence in the set of equation (5) is represented by the phase exponent $\exp(i\Delta\mathbf{q}\cdot\mathbf{r})$ which leads to convenient and symmetric form of the Takagi-Taupin equations used in the simulations.

Appendix B Takagi-Taupin equations for the modified amplitudes

In Takagi-Taupin equations (5) two coupling terms with $\chi_{h,\bar{h}}$ are responsible for dynamical diffraction effects. If we eliminated these terms from equations (5), they turn into a linear independent differential equations, which describe an independent transmission of the waves $E_{0,h}(\mathbf{r})$ through a crystal without diffraction. In this case, the analytical solution for each equation is represented by an exponential function $\exp(i\pi\chi_{0}s_{0,h}/\lambda)$.

Following this approach we substitute the amplitudes $E_{0,h}(\mathbf{r})$ in the Takagi-Taupin equations (5) by the new ones defined as (compare with Ref.Gorobtsov and Vartanyants (2016))

[TABLE]

which lead us to a new system of equations

[TABLE]

where the complex vector $\Delta\mathbf{q}^{\prime}$ is defined as

[TABLE]

This approach allows one to consider products in exponential factors as additives to the wave vectors and treat those as complex values with the directional properties given only by their real parts Authier (2003). As such the refraction and absorption of both diffracted and transmitted waves are included in their definition by applying boundary conditions for these amplitudes on a crystal surface. As it naturally follows from expressions (25), the modified amplitudes $E^{\prime}_{0,h}(\mathbf{r})$ differ from the amplitudes $E_{0,h}(\mathbf{r})$ inside the material, but are the same in a vacuum. Therefore, the boundary conditions for equations (26) should be expressed as

[TABLE]

where $\mathbf{R}_{b}$ is the radius vector of a considered point at the crystal-vacuum boundary, $E_{0,h}^{{}^{\prime}vac}(\mathbf{r})$ and $E_{0,h}^{{}^{\prime}cryst}(\mathbf{r})$ are the values of amplitudes in a vacuum and inside the crystal. In fact, that is equivalent to the condition of continuity of the tangential component of the electric field at the interface. Such approach allows to treat effects of refraction and absorption while propagating the wavefields in a crystal in a simple way.

Appendix C Numerical matrix form of the Takagi-Taupin equations

In the numerical integration method the complex amplitudes $E_{0,h}(\mathbf{r})$ are represented by a discrete set of values over all integration grid and the Takagi-Taupin equations are transformed to a recurrence matrix form, similar to Ref. Gronkowski (1991). Relying upon the symmetry of equations between the transmitted and the diffracted amplitudes, we take the same elementary integration step $p$ for both directions. For any smooth, slowly varying function $f(x,y)$ one can use the finite difference approximation to estimate the partial derivative from values at two neighboring points

[TABLE]

The value of the function in this middle point is given by a half sum

[TABLE]

When these formulas are applied to the differential equations (5) those are transformed to the following set

[TABLE]

with substitutions

[TABLE]

All the considered points belong to the same scattering plane, therefore in further derivations we simply omit the $s_{y}$ coordinate in the aid of shortness. The set of equations (31) can be further reorganized to

[TABLE]

Here two additional substitutions were made

[TABLE]

and it was assumed that $C\neq 0$, which is evidently true for any real positive $p$ as far as $\chi_{0r}\neq 0$. Solving this system with respect to $E_{0}(s_{0},s_{h})$ and $E_{h}(s_{0},s_{h})$ we obtain

[TABLE]

These relations can be also written in the matrix form

[TABLE]

where coefficients of matrix $M$ are expressed as

[TABLE]

Appendix D Kinematical limit of the Takagi-Taupin equations

Here we derive an analytical solution of the Takagi-Taupin equations for the purely kinematical case. In equations (5) we neglect coupling between the transmitted and diffracted components of the wave field, which is described by the term $\chi_{\bar{h}}$, we also neglect effects of refraction and absorption which are described by the term $\chi_{0}$. This leads to the following form of the Takagi-Taupin equations (5)

[TABLE]

The first equation can be easily solved as $E_{0}(\mathbf{r})=1$, where we assumed that the amplitude of the incoming beam is equal to unity. The amplitude of the diffracted wave field at the exit surface of the crystal can be obtained from the second equation in (38)

[TABLE]

where the following representation of the position vector $\mathbf{r}=s_{\bot}\mathbf{s}_{\bot}+s_{y}\mathbf{s}_{y}+s_{h}\mathbf{s}_{h}$ is used (see Fig. 2) and $S_{h}(\mathbf{r})$ is a complex crystalline function. As it was discussed before the far-field diffraction pattern can be obtained by the 2D Fourier transform of the exit surface wave

[TABLE]

where $q_{\bot},q_{y}$ are the reciprocal space coordinates. By substituting an expression (39) in equation (40) we obtain

[TABLE]

Since $ds_{\bot}ds_{y}ds_{h}=d\mathbf{r}$, the integral on the right side of equation (41) is the 3D Fourier transform of a complex crystalline function $S_{h}(\mathbf{r})=s_{h}(\mathbf{r})\exp(-i\mathbf{h}\cdot\mathbf{u}(\mathbf{r}))$. Comparison of equation (41) with equation (1) shows that they completely coincide, which gives a confidence that our approximations to the Takagi-Taupin equations indeed correspond to the kinematical diffraction case.