research papers
Experimentally obtained and computersimulated Xray asymmetric eightbeam pinhole topographs for a silicon crystal
^{a}NanoEngineering Research Center, Institute of Engineering Innovation, Graduate School of Engineering, The University of Tokyo, 21116 Yayoi, Bunkyoku, Tokyo 1138656, Japan, ^{b} Japan Synchrotron Radiation Research Institute, SPring8, 111 Kouto, Mikazukicho, Sayogun, Hyogo 6795198, Japan, and ^{c}Rigaku Corporation, 3912 Matsubara, Akishimashi, Tokyo 1968666, Japan
^{*}Correspondence email: okitsu@soyak.t.utokyo.ac.jp
This article is dedicated to Professor K. Kohra, who passed away on 29 January 2019.
In this study, experimentally obtained eightbeam pinhole topographs for a silicon crystal using synchrotron Xrays were compared with computersimulated images, and were found to be in good agreement. The experiment was performed with an asymmetric allLaue geometry. However, the Xrays exited from both the bottom and side surfaces of the crystal. The simulations were performed using two different approaches: one was the integration of the nbeam Takagi–Taupin equation, and the second was the fast Fourier transformation of the Xray amplitudes obtained by solving the eigenvalue problem of the nbeam Ewald–Laue theory as reported by Kohn & Khikhlukha [Acta Cryst. (2016), A72, 349–356] and Kohn [Acta Cryst. (2017), A73, 30–38].
Keywords: Xray diffraction; dynamical theory; multiple reflection; computer simulation; nbeam reflection; phase problem; silicon; protein crystallography.
1. Introduction
We previously reported a derivation of the nbeam Takagi–Taupin (TT) equation and an algorithm to integrate it (Okitsu, 2003; Okitsu et al., 2006). We verified these by comparing computersimulated and experimentally obtained topographs using a sixbeam case (Okitsu et al., 2003, 2006, 2011) and three, four, five, six, eight and 12beam 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. 2006, O et al. 2011 and OIY 2012, respectively.
In OIY 2012, the nbeam TT equation was derived from the nbeam Ewald–Laue (EL) 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 Xray rocking curves that were obtained by fast Fourier transformation of the Xray amplitude in a threebeam topograph, and compared them with those computed by solving the eigenvalue problem of the threebeam EL theory. In contrast, Heyroth et al. (2001) reported Xray threebeam topographs experimentally obtained and computer simulated by coherently superposing the Xray amplitude calculated based on the EL theory.
Recently, Kohn & Khikhlukha (2016) and Kohn (2017) reported computersimulated nbeam topographs (n = 6) obtained by fast Fourier transformation of the rocking amplitude calculated using the nbeam EL theory. The present article reports the efficiency of Kohn's method by comparing computersimulated eightbeam pinhole topographs (EL&FFT simulation) with the experimentally obtained and computersimulated ones based on the nbeam TT equation (TT simulation) published in OIY 2012 {Figs. 5 [S_{x}(TT)], 5 [E_{x}], 6 [S_{x}(TT)] and 6 [E_{x}] ()}.
2. 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 sixbeam pinhole topographs. However, the goniometer axes were adjusted such that the 000 forwarddiffracted (FD) and 004, 026, 066, 084, 080, and transmittedreflected (TR) Xrays 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 Xrays with a photon energy of 18.245 keV at BL09XU of SPring8 was controlled by using a rotating fourquadrant 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 Xray beam whose dimensions were limited to 25 × 25 µm 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 µm was placed 47.5 mm behind the crystal, such that its surface was approximately perpendicular to the [100] direction of the crystal.
3. Computer simulation
3.1. Integration of the nbeam Takagi–Taupin equation
The TT 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 µm. The height of the octagonal Borrmann pyramid was assumed to be 19.015 mm as calculated by mm, where , and are unit vectors whose directions are as drawn in Figs. 1 and 4. The integration of the nbeam TT 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) in O et al. 2006 with the lattice displacement term omitted. Here, D_{i}^{(l)} and D_{j}^{(m)} are the amplitudes of the ithnumbered lpolarized and jthnumbered mpolarized Xrays, where and l, m ∈ {0, 1}. is the thorder Fourier coefficient of the electric susceptibility. C_{i, j}^{(l, m)} is the polarization factor as defined later in equation (7). The TT 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 E52680v3) 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, 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.

3.2. Fast Fourier transformation of the rocking amplitude calculated using the Ewald–Laue theory
The nbeam Xray amplitude when the crystal is rotated twodimensionally in the vicinity of the exact nbeam condition, with an incidence of planewave Xrays, can be obtained by solving the eigenvalue problem of the EL theory.
In the TT simulations reported by Okitsu and coauthors (O et al. 2006, O et al. 2011 and OIY 2012), a boundary condition that the incident Xrays have a nonzero amplitude only at the incident point on the crystal was assumed, i.e. the incident Xray amplitude of the delta function was assumed as the boundary condition. The coincidence between the experimentally obtained and TT 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 is obtained by the Fourier transformation of the delta function in real space, the Xray amplitudes of the nbeam pinhole topographs can also be obtained by Fourier transforming the Xray amplitudes computed by solving the eigenvalue problem of the EL theory.
EL&FFT simulations of the nbeam section topography for n = 6 were reported by Kohn & Khikhlukha (2016) and Kohn (2017) for a symmetric sixbeam Laue case. In this case, () have identical values, where is a unit vector in the direction of the ithbeam propagation and is a unit vector in the direction of the downward surface normal to the crystal. However, the eightbeam 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 EL&FFT simulation in the present work can be described as follows.
When planewave Xrays are incident on a parallel plate crystal, a Bloch wave is excited as follows:
Here, n is the number of waves, , , is the ithnumbered reflection vector, P_{1}′ is the common initial point of the wavevector of the Bloch wave and is the location vector. Laue's fundamental equation of the (Laue, 1931; Authier, 2005) restricts the amplitude and wavevector of the Bloch wave as follows:
Here, K = 1/λ, where λ is the wavelength of the Xrays in vacuum and is the component vector of perpendicular to . By applying the approximation k_{i} + K ≃ 2K, equation (3) becomes
When is the direction of polarization of the ithnumbered lpolarized Xray beam,
was defined in the present work such that = and . A form of the nbeam EL 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 secondorder 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 . The values of Θ_{i} for Fig. 1(b) α and 1(c) β are summarized in Table 1. is the downward surface normal of the exit surface (but not the entrance surface as discussed later) of the crystal. The directions of for Figs. 1(b) and 1(c) are and , respectively, as shown in the upperright corner of Fig. 1.
β^{(0)} and β^{(1)} are angular deviations of P_{1} from the Laue point, La, i.e. . 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 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 and vector as follows:
Here, is a 2norder column vector whose qth (q = 2j + m + 1) element is , and is a 2n × 2n matrix whose element of the pth (p = 2i + l + 1) row and qth column, , is as follows:
Here, δ_{p,q} is the Kronecker delta. When the Xray 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 (). After obtaining these, the following equation is solved to satisfy the boundary condition that depends on the polarization state of the incident Xrays:
Here, is a 2n × 2n matrix whose element of the pth row and kth column is , is a column vector whose kth element is , and is a column vector of the boundary condition that depends on the polarization state of the incident Xrays. For 0 and 1polarized incident Xrays, is and , respectively.
Incidentally, when P′_{1,k} is the common initial point of the kth Bloch wave, . Here, let point P_{1}′′ be defined such that and , whereas and are unit vectors defined such that and , and form a righthanded orthogonal system in this order (see Fig. 4 in OIY 2019). The total wavefield is excited by the incident planewave Xrays, and is given by
where is a location vector in the crystal. The amplitude of the ithnumbered Xray with polarization state l on the exit surface, , of the crystal, , is given by
Here, where T_{z} is the thickness of the crystal. Let the amplitude depend on two scalar values, x and y, then
Since , and (see Fig. 4 in OIY 2019), then
Here, has been defined as follows:
The wavefield D_{i}^{(l)}(x, y) excited by the incident Xrays whose wavefront is the delta function is a coherent superposition of the wavefield which is excited by the incident planewave Xrays with an amplitude of unity at the entrance surface of the crystal. Therefore,
The term in the above equation is necessary to separately calculate the Xray amplitudes of α_{1} and β_{1} in Fig. 1(a), as shown in Figs. 8(α_{1}) and 8(β_{1}) {from which Figs. 5 [S_{v}(EL)] and 6 [S_{v}(EL)] have been obtained}. The following assumptions are made: and T_{z} = 9.6 mm for Fig. 1(b); and and T_{z} = 16.5 mm for Fig. 1(c). The term exists in equation (17) because the projection of the twodimensional integration element over a plane normal to should be the same, even for different directions of . The calculated values of are 1.60518 and 1.74046 for Figs. 1(b) and 1(c), respectively. Since for 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 is a twodimensional periodic function with a period 1/L, the contents in the summation of the righthand side of equation (19) are evidently also a twodimensional periodic function with a period N. Therefore, D_{i}^{(l)}(n_{x}L/N, n_{y}L/N) in the lefthand side of equation (18) is also a twodimensional periodic function with a period of L. Then, equation (19) can be replaced by the following equation:
Equation (20) has a general form of a twodimensional fast Fourier transform (FFT) (Cooley & Tukey, 1965). Therefore, D_{i}^{(l)}(n_{x}L/N, n_{y}L/N) can be obtained using an FFT as defined by equation (20). However, before performing the FFT, [− N/2 ≤ {k_{x}, k_{y}} ≤ N/2 − 1] in equation (19) should be replaced with [0 ≤ {k_{x}, k_{y}} ≤ N − 1] in equation (20). Similarly after performing the FFT, D_{i}^{(l)}′(n_{x}L/N, n_{y}L/N) [0 ≤ {n_{x}, n_{y}} ≤ N − 1] in equation (20) should be replaced with D_{i}^{(l)}(n_{x}L/N, n_{y}L/N) [− N/2 ≤ {n_{x}, n_{y}} ≤ N/2 − 1] in equation (19). [S_{h}(EL)] and [S_{v}(EL)] shown in 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) was calculated with the FFT routine in the Intel Math Kernel Library (MKL). It took 470 s (280 s to solve the eigenvalue problem, 20 s for the FFT and 170 s to write the topographs to the hard disk) to obtain eight topograph images as shown in Fig. 5 [S_{h}(EL)] or [S_{v}(EL)] using one node (Intel Xeon E52680v3) of the supercomputer system `Sekirei' of the Institute for Solid State Physics at the University of Tokyo.
4. Results
Fig. 5 [E_{x}] (x ∈ {h, v}) shows the experimentally obtained pinhole topograph images recorded on the IP for the incidence of the horizontally polarized (x = h) and vertically polarized (x = v) Xrays. Fig. 5 [S_{x}(TT)] and [S_{x}(EL)] show the TT and EL&FFT simulated images corresponding to Fig. 5 [E_{x}]. In Figs. 6 and 7, enlargements of the 000 FD and 066 TR images, respectively, are shown. Fig. 6 [E_{h}] and [E_{v}] are reproductions of Fig. 11 [S(a)] and [S(b)] in OIY 2012. Fig. 6 [S_{h}(TT)] and [S_{v}(TT)] are obtained by solving the nbeam TT equation layer by layer with a thickness of (19.015/4000) mm, whereas the number of layers was 3600 for Fig. 11 [S(a)] and [S(b)] in OIY 2012. In Figs. 6 and 7, there is good agreement between the experimentally obtained and computersimulated topographs (both the TT and EL&FFT simulated ones).
A `harpshaped' pattern (HpSP), a pattern whose shape is similar to the alphabetical character `Y' (YSP) and `nailshaped' patterns (NSP) in Fig. 6 [E_{h}] are also shown in both Fig. 6 [S_{h}(TT)] and [S_{h}(EL)], revealing the equivalence between the TT and EL&FFT simulations. NSPs were also found in Fig. 6 [S_{v}(TT)], [E_{v}] and [S_{v}(EL)]. The HpSPs in Fig. 6 [S_{v}(TT)], [E_{v}] and [S_{v}(EL)] are practically the same but fainter than those in Fig. 6 [E_{h}], [S_{h}(TT)] and [S_{h}(EL)], which shows an evident discrepancy owing to the polarization state of the incident Xrays. 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 computersimulated images of Fig. 7. However, in Fig. 7, the HpSP for the incidence of vertically polarized Xrays is evidently fainter than that of the horizontally polarized Xrays.
In Fig. 7 [E_{x}] and [S_{x}(EL)] (x ∈ {h, v}), the central region with a relatively strong seems to be surrounded by a `veil' with small Xray intensities (Veil). However, such a faint region like a `veil' is absent in the TT simulated topographs.
As described in the last paragraph of Section 3, the calculation speed of the EL&FFT simulation was over 100 times faster than the TT simulation. The calculation speed of the EL&FFT simulation is constant and independent of the crystal thickness, whereas that of the TT 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.
5. Discussion
In Fig. 6 [S_{h}(TT)], a sharp line similar to a knife edge (KEL) is observed. However, a KEL is not seen in Fig. 6 [E_{h}] or [S_{h}(EL)]. The width of the KEL is extremely narrow. In the case of the TT simulation, a boundary condition of the incident Xrays whose amplitude is the delta function was assumed. Then the incident Xrays have a planewave component whose initial point of the wavevector was far from the Laue point. However, the incident Xrays used in the experiment have a finite angular width. In addition, in the EL&FFT simulation, the integration range is finite. This is probably the reason for the absence of the KEL in Fig. 6 [E_{h}] and [S_{h}(EL)].
With regard to the Veil, this pattern may be explained by the following hypothesis. When the crystal is thick in the case of twobeam section topography in general, dark areas are observed in the forwarddiffracted and transmittedreflected topographs on both sides of the central bright region due to the The Veil may correspond to these dark areas excited by the incident Xray planewave 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 () in Fig. 5 [E_{x}] and [S_{x}(EL)] (x ∈ {h, v}). However, it cannot be found in Fig. 5 [S_{x}(TT)] (x ∈ {h, v}). This feature of the Veil may be explained by the weak or zero intensity of the Xray planewave component whose incident angle is far from the exact eightbeam condition.
When the rocking curves of the Xrays are discussed in the EL 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}(EL)]. When performing the EL&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 on the exit surface of the crystal does not depend on the shape of the entrance side of the crystal.
The Xray nbeam rocking amplitude from a planar perfect crystal can also be obtained by solving the nbeam TT equation, which will be reported in the near future. Since the result of the EL&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 TT equation.
6. Summary
Experimentally obtained and computersimulated (both TT simulated and EL&FFT simulated) asymmetric eightbeam 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 Xray wavefield could be computed not only based on the nbeam TT equation but also on the nbeam EL theory.
The present work has provided the first demonstration of the EL&FFT simulation overcoming difficulties when calculating the Xray 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 analysis.
Acknowledgements
The SGI ICE XA supercomputer system `Sekirei', consisting of Intel Xeon E52680v3 processors, of the Institute for Solid State Physics at the University of Tokyo was used for the computer simulations. The authors are indebted to Dr X.W. Zhang and Dr T. Oguchi for their technical support in the present experiments, and also to Professor Emeritus S. Kikuta for his encouragement and fruitful discussions with regard to the present work.
Funding information
The theoretical component and computer simulation of the present work were supported by the Nanotechnology Platform Project (No. 12024046) of the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. The preliminary experiments were performed at ARBL03A of the Photon Factory AR under the approval of the Photon Factory Program Advisory Committee (Proposal Nos. 2003G202 and 2003G203). The main experiments were performed at BL09XU of SPring8 under the approval of the Japan Synchrotron Radiation Research Institute (JASRI) (Proposal No. 2005B0714).
References
Authier, A. (2005). Dynamical Theory of Xray Diffraction, reprinted with revisions 2004, 2005. Oxford University Press. Google Scholar
Chang, S.L. (2004). Xray MultipleWave Diffraction, Theory and Application. Berlin, Heidelberg, New York: Springer. Google Scholar
Colella, R. (1974). Acta Cryst. A30, 413–423. CrossRef CAS IUCr Journals Web of Science Google Scholar
Cooley, J. W. & Tukey, J. W. (1965). Math. Comput. 19, 297–301. CrossRef Web of Science Google Scholar
Heyroth, F., Zellner, J., Höche, H., Eisenschmidt, C., Weckert, E. & Drakopoulous, M. (2001). J. Phys. D Appl. Phys. 34, A151–A157. CrossRef CAS Google Scholar
Ishiwata, G., Okitsu, K. & Ishiguro, M. (2010). Acta Cryst. A66, 484–488. Web of Science CrossRef IUCr Journals Google Scholar
Kohn, V. G. (2017). Acta Cryst. A73, 30–38. CrossRef IUCr Journals Google Scholar
Kohn, V. G. & Khikhlukha, D. R. (2016). Acta Cryst. A72, 349–356. Web of Science CrossRef IUCr Journals Google Scholar
Laue, M. von (1931). Ergeb. Exakten Naturwiss. 10, 133–158. CrossRef Google Scholar
Okitsu, K. (2003). Acta Cryst. A59, 235–244. Web of Science CrossRef CAS IUCr Journals Google Scholar
Okitsu, K., Imai, Y., Ueji, Y. & Yoda, Y. (2003). Acta Cryst. A59, 311–316. Web of Science CrossRef CAS IUCr Journals Google Scholar
Okitsu, K., Imai, Y. & Yoda, Y. (2012). Recent Advances in Crystallography, pp. 67–86. IntechOpen. https://dx.doi.org/10.5772/47846. Google Scholar
Okitsu, K., Imai, Y. & Yoda, Y. (2019). Acta Cryst. A75, 483–488. CrossRef IUCr Journals Google Scholar
Okitsu, K., Yoda, Y., Imai, Y. & Ueji, Y. (2011). Acta Cryst. A67, 550–556. Web of Science CrossRef IUCr Journals Google Scholar
Okitsu, K., Yoda, Y., Imai, Y., Ueji, Y., Urano, Y. & Zhang, X. (2006). Acta Cryst. A62, 237–247. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sanchez del Rio, M. & Dejus, R. J. (1998). Proc. SPIE, 3448, 340–345. CrossRef Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.