Experimentally obtained and computer-simulated X-ray asymmetric eight-beam pinhole topographs for a silicon crystal

Experimentally obtained eight-beam pinhole topographs for a silicon crystal were compared with computer simulations based on the n-beam Takagi–Taupin equation and Ewald–Laue theory.


Introduction
We previously reported a derivation of the n-beam Takagi-Taupin (T-T) equation and an algorithm to integrate it (Okitsu, 2003;Okitsu et al., 2006). We verified these by comparing computer-simulated and experimentally obtained topographs using a six-beam case (Okitsu et al., , 2006(Okitsu et al., , 2011 and three-, four-, five-, six-, eight-and 12-beam cases (Okitsu et al., 2012). Hereafter, Okitsu et al. (2006), Okitsu et al. (2011) and Okitsu et al. (2012) are denoted by O et al. 2006O et al. , O et al. 2011O et al. and OIY 2012 In OIY 2012, the n-beam T-T equation was derived from the n-beam Ewald-Laue (E-L) theory, and vice versa by their Fourier transformation, which explicitly revealed a simple relationship between them described by a Fourier transform. Ishiwata et al. (2010) reported X-ray rocking curves that were obtained by fast Fourier transformation of the X-ray amplitude in a three-beam topograph, and compared them with those computed by solving the eigenvalue problem of the three-beam E-L theory. In contrast, Heyroth et al. (2001) reported X-ray three-beam topographs experimentally obtained and computer simulated by coherently superposing the X-ray amplitude calculated based on the E-L theory.
Recently, Kohn & Khikhlukha (2016) and Kohn (2017) reported computer-simulated n-beam topographs (n = 6) obtained by fast Fourier transformation of the rocking amplitude calculated using the n-beam E-L theory. The present article reports the efficiency of Kohn's method by comparing computer-simulated eight-beam pinhole topographs (E-L&FFT simulation) with the experimentally obtained and computer-simulated ones based on the n-beam T-T equation (T-T

Experimental
The optics used in the present work, shown in Fig. 1(a), were fundamentally the same as those in Fig. 2, which is a reproduction of Fig. 7 of O et al. 2006, showing the experimental setup used when taking the six-beam pinhole topographs. However, the goniometer axes were adjusted such that the 000 forward-diffracted (FD) and 004, 026, 066, 084, 080, 062 and 022 transmitted-reflected (TR) X-rays were simultaneously strong, and the vector product of the 000 FD and 066 TR beam directions was horizontal. The polarization state of the incident synchrotron X-rays with a photon energy of 18.245 keV at BL09XU of SPring-8 was controlled by using a rotating four-quadrant phase retarder system. Its schematic and photograph are shown in Figs. 3(a) and 3(b), respectively [these are reproductions of Figs. 3(a) and 3(b) of OIY 2012]. Its usage was described in O et al. 2006. An X-ray beam whose dimensions were limited to 25 Â 25 mm was incident on a position on the entrance surface of a floating zone (FZ) silicon crystal with a thickness of 9.6 mm. The incident point was 16.5 mm from the corner of the crystal block, as shown in Fig.  1(a). The orientation of the crystal is also shown in Fig. 1. An imaging plate (IP) with a pixel size of 50 Â 50 mm was placed 47.5 mm behind the crystal, such that its surface was approximately perpendicular to the [100] direction of the crystal.

Integration of the n-beam Takagi-Taupin equation
The T-T simulations were performed in a similar manner to the approach described in O et al. 2006, except that the crystal was divided into small octagonal pyramids, as shown in Fig.  4(b). The calculated value of l 1 in Fig. 4 was 29.493 mm. The height of the octagonal Borrmann pyramid was assumed to be 19.015 mm as calculated by ½z b ð16:5x c þ 9:6z c Þ mm, where z b , x c and z c are unit vectors whose directions are as drawn in Figs. 1 and 4. The integration of the n-beam T-T equation was performed by solving equation (1) [see Fig. 4(a)] layer by layer whose thickness and normal direction were (19.015/4000) mm and [100], respectively: where n = 8. Equation (1) is a reproduction of equation (8)  Reproduction of Fig. 7 in Okitsu et al. (2006), showing the six-beam pinhole topographs. An identical floating zone silicon crystal with a thickness of 9.6 mm was also used in the present experiment. However, the angles of the goniometers were adjusted such that the 000 forwarddiffracted (FD) and 004, 026, 066, 084, 080, 062 and 022 transmittedreflected (TR) X-rays were simultaneously strong, as shown in Fig. 1(a). l-polarized and jth-numbered m-polarized X-rays, where i; j 2 f0; 1; 2; . . . ; 7g and l, m 2 {0, 1}. h i Àh j is the ðh i À h j Þth-order Fourier coefficient of the electric susceptibility. C i, j (l, m) is the polarization factor as defined later in equation (7). The T-T simulation took 15 h using one node of the supercomputer system 'Sekirei' of the Institute for Solid State Physics of the University of Tokyo. Each node (Intel Xeon E5-2680v3) had 24 cores. The program was then parallelized using Fortran 90 with MPI (Message Passing Interface).
As shown in Fig. 1(a), two vacuum parts are included in the Borrmann pyramid. Then, h i Àh j ði; j 2 f0; 1; 2; . . . ; 7gÞ are assumed to be zero in the vacuum regions, as described in O et al. 2011. The reflection parameters as calculated using XOP 2.3 (Sanchez del Rio & Dejus, 1998) are used for the simulation, as summarized in Table 1.

Fast Fourier transformation of the rocking amplitude calculated using the Ewald-Laue theory
The n-beam X-ray amplitude when the crystal is rotated two-dimensionally in the vicinity of the exact n-beam condition, with an incidence of planewave X-rays, can be obtained by solving the eigenvalue problem of the E-L theory.
In the T-T simulations reported by Okitsu and co-authors (O et al. 2006, O et al. 2011and OIY 2012, a boundary condition that the incident X-rays have a nonzero amplitude only at the incident point on the crystal was assumed, i.e. the incident X-ray amplitude of the delta function was assumed as the boundary condition. The coincidence between the experimentally obtained and T-T simulated results implied that the physical properties of the pinhole topography could be approximated by the boundary condition of the delta function. Because a function of unity with an identical phase in reciprocal space is obtained by the Fourier transformation of the delta function in real space, the X-ray amplitudes of the n-beam pinhole topographs can also be obtained by Fourier transforming the X-ray amplitudes computed by solving the eigenvalue problem of the E-L theory. E-L&FFT simulations of the n-beam section topography for n = 6 were reported by Kohn & Khikhlukha (2016) and Kohn (2017) for a symmetric six-beam Laue case. In this case, s i Á n z (i 2 f0; 1; . . . ; n À 1g; n ¼ 6) have identical values, where s i is a unit vector in the direction of the ith-beam propagation and n z is a unit vector in the direction of the downward surface normal to the crystal. However, the eight-beam pinhole topography reported in the present work was performed for asymmetric Laue geometry, as shown in Fig. 1(a). Furthermore, the exit surface of the crystal was not a single plane. Accordingly, the procedure of the E-L&FFT simulation in the present work can be described as follows.
When plane-wave X-rays are incident on a parallel plate crystal, a Bloch waveD D is excited as follows: Here, n is the number of waves, i 2 f0; 1; 2; . . . ; n À 1g, 0 is the common initial point of the wavevector of the Bloch wave and r is the location vector. Laue's fundamental equation of the dynamical theory (Laue, 1931;Authier, 2005) restricts the amplitude and wavevector of the Bloch wave as follows: Here, K = 1/, where is the wavelength of the X-rays in vacuum and ½D j ?k i is the component vector of D j perpendi-cular to k i . By applying the approximation k i + K ' 2K, equation (3) becomes When e ðlÞ i is the direction of polarization of the ith-numbered l-polarized X-ray beam, e ðlÞ i was defined in the present work such that e ð0Þ the n-beam E-L theory that is applicable to asymmetric Laue geometry is found in Section 7.1.2 of Chang's book (Chang, 2004). An algorithm to solve the eigenvalue problem was first developed by Colella (1974) and is more complex as it considers the second-order term. A simpler eigenvalue problem with a linear approximation found in Chang's book can be described as follows: Here, the polarization factors C and S are defined as Equation (6) is the standard form of an eigenvalue problem, where cos Â i ¼ s i Á n z . The values of Â i for Fig. 1(b) and 1(c) are summarized in Table 1. n z is the downward surface normal of the exit surface (but not the entrance surface as discussed later) of the crystal. The directions of n z for Figs. 1(b) and 1(c) are ½111 and ½211, respectively, as shown in the upper-right corner of Fig. 1. (0) and (1) are angular deviations of P 1 from the Laue point, La, i.e. P 1 La ! ¼ Kð ð0Þ e ð0Þ 0 þ ð1Þ e ð1Þ 0 Þ. P 1 is a point that exists on a spherical surface whose distance from H 0 is K, where H 0 is the origin of the reciprocal space. This spherical surface is assumed to be an approximate plane whose distance from H 0 is K in the present work. The geometry around the Laue point is the same as that of Fig. 4 in Okitsu et al. (2019) (hereafter denoted OIY 2019). Equation (6) can be described using the matrix A and vector D as follows: Here, D is a 2n-order column vector whose qth (q = 2j + m + 1) element is D q , and A is a 2n Â 2n matrix whose element of the pth (p = 2i + l + 1) row and qth column, A p;q , is as follows:   Here, p,q is the Kronecker delta. When the X-ray amplitudes in 1 and 2 , and those in 1 and 2 , in Figs. 1(b) and 1(c), respectively, are calculated, Â i () and Â i () (summarized in Table 1) were substituted for Â i in equations (6) and (9). In general, equation (8) has 2n couples of eigenvalues k and eigenvectors D k (k 2 f1; 2; Á Á Á ; 2ng). After obtaining these, the following equation is solved to satisfy the boundary condition that depends on the polarization state of the incident X-rays: Here, D is a 2n Â 2n matrix whose element of the pth row and kth column is D p;k ð¼ D ðlÞ i;k Þ, C is a column vector whose kth element is C k , and B is a column vector of the boundary condition that depends on the polarization state of the incident X-rays. For 0-and 1-polarized incident X-rays, B is ð1; 0; 0; . . . ; 0Þ T and ð0; 1; 0; . . . ; 0Þ T , respectively.
Incidentally, when P 0 1,k is the common initial point of the kth Bloch wave, P 0 1;k P 1 ! ¼ k n z . Here, let point P 1 00 be defined such that P 1 P 00 Ák y e y , whereas e x and e y are unit vectors defined such that e x ¼ e ð0Þ 0 and e x , e y and n z form a right-handed orthogonal system in this order (see Fig. 4 in OIY 2019). The total wavefield D D total ðrÞ is excited by the incident plane-wave X-rays, and is given byD where r is a location vector in the crystal. The amplitude of the ith-numbered X-ray with polarization state l on the exit surface, r exit , of the crystal, D ðlÞ i ðr exit Þ, is given by Here, r exit ¼ xe x þ ye y þ T z n z where T z is the thickness of the crystal. Let the amplitude D ðlÞ i ðx; yÞ depend on two scalar values, x and y, then Since P 0 1;k P 1 ! ¼ k n z , P 1 P 00 1 ! ¼ 0 n z and P 00 1 La ! ¼ Ák x e x þ Ák y e y (see Fig. 4 in OIY 2019), then Here, D ðlÞ0 i ðÁk x ; Ák y Þ has been defined as follows: The wavefield D i (l) (x, y) excited by the incident X-rays whose wavefront is the delta function is a coherent superposition of the wavefield D ðlÞ i ðx; yÞ which is excited by the incident plane-wave X-rays with an amplitude of unity at the entrance surface of the crystal. Therefore, The term 1=ðs 0 Á n z Þ in the above equation is necessary to separately calculate the X-ray amplitudes of 1 and 1 in Fig.  1(a), as shown in Figs. 8( 1 ) and 8( 1 ) {from which Figs. 5 [S v (E-L)] and 6 [S v (E-L)] have been obtained}. The following assumptions are made: n z ¼ z c and T z = 9.6 mm for Fig. 1(b); and n z ¼ x c and T z = 16.5 mm for Fig. 1(c). The term 1=ðs 0 Á n z Þ exists in equation (17) because the projection of the twodimensional integration element over a plane normal to s 0 should be the same, even for different directions of n z . The calculated values of 1=ðs 0 Á n z Þ are 1.60518 and 1.74046 for Figs. 1(b) and 1(c), respectively. Since r exit for D ðlÞ i ðr exit Þ has a nonzero value only if it is inside the Borrmann pyramid, an infinitesimal spatial resolution is not necessary, and equation (17) can be replaced with a discrete Fourier transform as follows: Here, n x , n y , k x and k y are integers, N is an even integer with which the lower and upper limits of the summations are determined, and L Â L is the field size of a square that includes the topograph image. Therefore, Ák x = k x /L, Ák y = k y /L, x = n x L/N and y = n y L/N. If D ðlÞ0 i ðÁk x ; Ák y Þ is a twodimensional periodic function with a period 1/L, the contents in the summation of the right-hand side of equation (19) are evidently also a two-dimensional periodic function with a period N. Therefore, D i (l) (n x L/N, n y L/N) in the left-hand side of equation (18) is also a two-dimensional periodic function with a period of L. Then,equation (19) can be replaced by the following equation: Equation (20) Figs. 5, 6, 7, 8( 1 ) and 8( 1 ) were obtained using the above procedure with L = 60 mm and N = 4096.
The eigenvalue/eigenvector problem described by equations (8) and (9) was solved using ZGEEV of LAPACK. The Fourier transform described by equation (20)   , which shows an evident discrepancy owing to the polarization state of the incident X-rays. An HpSP is also found in an elliptical region (Ellip) in Fig. 7 [E h ]. HpSPs and Ellips also exist in all of the experimentally obtained and computer-simulated images of Fig. 7. However, in Fig. 7, the HpSP for the incidence of vertically polarized X-rays is evidently fainter than that of the horizontally polarized X-rays.

Results
In Fig. 7 [E x ] and [S x (E-L)] (x 2 {h, v}), the central region with a relatively strong X-ray intensity seems to be surrounded by a 'veil' with small X-ray intensities (Veil). However, such a faint region like a 'veil' is absent in the T-T simulated topographs.
As described in the last paragraph of Section 3, the calculation speed of the E-L&FFT simulation was over 100 times faster than the T-T simulation. The calculation speed of the E-L&FFT simulation is constant and independent of the crystal thickness, whereas that of the T-T simulation is proportional to the reciprocal of the thickness cubed. The superiority of this method reported by Kohn & Khikhlukha (2016) and Kohn (2017) was verified for a perfect crystal.

Discussion
In Fig. 6 [S h (T-T)], a sharp line similar to a knife edge (KEL) is observed. However, a KEL is not seen in Fig. 6 The width of the KEL is extremely narrow. In the case of the T-T simulation, a boundary condition of the incident X-rays whose amplitude is the delta function was assumed. Then the incident X-rays have a plane-wave component whose initial point of the wavevector was far from the Laue point. However, the incident X-rays used in the experiment have a finite angular width. In addition, in the E-L&FFT simulation, the integration range is finite. This is probably the reason for the absence of the KEL in Fig. 6 With regard to the Veil, this pattern may be explained by the following hypothesis. When the crystal is thick in the case of two-beam section topography in general, dark areas are observed in the forward-diffracted and transmitted-reflected topographs on both sides of the central bright region due to Enlargements of the 000 forward-diffracted images in Fig. 5. the Borrmann effect. The Veil may correspond to these dark areas excited by the incident X-ray plane-wave component whose initial point of the wavevector is distant from the Laue point. The Veil can be observed in all h i k i l i -diffracted images (i 2 f0; 1; 2; . . . ; 7g) in Fig. 5 [E x ] and [S x (E-L)] (x 2 {h, v}). However, it cannot be found in Fig. 5 [S x (T-T)] (x 2 {h, v}). This feature of the Veil may be explained by the weak or zero intensity of the X-ray plane-wave component whose incident angle is far from the exact eight-beam condition.
When the rocking curves of the X-rays are discussed in the E-L theory in general, the downward surface normal is assumed to be perpendicular to the entrance surface of the crystal. Its directions for Figs. 1(b) and 1(c) are perpendicular to each other. However, Figs. 8( 1 ) and 8( 1 ) are smoothly linked as shown in Fig. 6 [S v (E-L)]. When performing the E-L&FFT simulation, the direction of the entrance surface of the crystal does not have to be considered, whereas the direction of the exit surface was important. Thus, the X-ray intensity on the exit surface of the crystal does not depend on the shape of the entrance side of the crystal.
The X-ray n-beam rocking amplitude from a planar perfect crystal can also be obtained by solving the n-beam T-T equation, which will be reported in the near future. Since the result of the E-L&FFT simulation is verified to be independent of the shape of the entrance surface, the same result as reported in the present work should be obtained by the fast Fourier transform of the rocking amplitude calculated by the T-T equation.

Summary
Experimentally obtained and computer-simulated (both T-T simulated and E-L&FFT simulated) asymmetric eight-beam pinhole topographs, which were in good agreement, were reported. This is for a case where the exit surface was not a single plane. It was verified that the X-ray wavefield could be computed not only based on the n-beam T-T equation but also on the n-beam E-L theory.
The present work has provided the first demonstration of the E-L&FFT simulation overcoming difficulties when calculating the X-ray intensities diffracted from such a complexshaped crystal as shown in Fig. 9 to verify the first hypothesis concerning an excessively large R factor in a protein crystal structure analysis.    ( 1 ) are computed separately under the assumption of vertically polarized incident X-rays. These figures have been computed by projecting intensities of the 000 forward-diffracted X-rays on the exit planes 1 and 1 in Fig. 1(a) on the IP whose surface was normal to the [100] direction. X-ray intensities of 2 and 2 in Figs. 1(b) and 1(c) have been erased.

Figure 9
Outline drawing of a tetragonal lysozyme crystal with 12 facets.