Many synchrotron beamlines that produce hard X-rays (photon energy ≥ 3 keV) select a suitable wavelength by using diffraction from perfect silicon crystal monochromators. To select photon energies up to about 20 keV, diffraction in the Bragg case is used. In this case the diffracted beam emerges from the surface upon which the incident beam falls. By using perfect crystals one can obtain a resolving power δE/E ≈ 10-4 in the photon energy E, which is suitable for many types of experiments. In the following, practical examples of the application of the Bragg case will be discussed.
X-ray diffraction from large perfect crystals must be treated with the theory of dynamical diffraction, as the simpler kinematic theory does not account for extinction of the incident wave or for the coherent coupling of the incident and diffracted waves inside the crystal. Dynamical diffraction of plane waves from perfect crystals was treated at length in the early twentieth century by Darwin,1 Prins,2 Ewald3 and von Laue4; reviews of their theory have been provided by Zachariasen5 and by Batterman and Cole6 among others. Dynamical diffraction of spherical waves from absorbing perfect crystals was subsequently treated by Kato7,8 in the Laue case and by Saka, Katagawa and Kato9 in the Bragg case. Although spherical-wave dynamical diffraction has not been widely applied to synchrotron beamline modeling, the formulas for plane-wave reflectivity are routinely used for calculations of rocking curves by XOP10 and in ray traces by SHADOW.11 However, despite this long history, there is no widely available software to calculate the complete complex amplitude distribution of a synchrotron X-ray beam diffracted by a crystal. Not only the widespread use of perfect crystal monochromators, but also the increasing need for diffraction-limited focusing, make it worthwhile to fill this gap.
The software package “Synchrotron Radiation Workshop,” or SRW, was written by Chubar and Elleaume. SRW propagates X-ray wavefields through any combination of drift spaces, apertures, lenses or short mirrors. However, SRW could not include perfect crystals until now. Here the integration of a new module for Bragg-case dynamical diffraction from a perfect crystal is presented. This module has been tested in the extreme cases of approximately plane-wave and spherical-wave incident beam diffracted either vertically or horizontally by Bragg reflections of arbitrary asymmetry.
The perfect crystal propagator is implemented in the SRW C/C++ library. SRW-Python, the latest version of SRW, operates under Python 2.7 and higher. Python 3.3.3 was used to run the simulations of this article. The perfect crystal is defined in SRW-Python by the class SRWLOptCryst. The user supplies the following parameters as input:
• _d_sp, the interplanar spacing d of the selected Bragg reflection
• _psi0r, 4 π × the real part of the Fourier component Ψ0 of the electric susceptibility
• _psi0i, 4 π×the imaginary part of the Fourier component Ψ0 of the electric susceptibility
• _psi_hr, 4 π×the real part of the Fourier component ΨH of the electric susceptibility
• _psi_hi, 4π×the imaginary part of the Fourier component ΨH of the electric susceptibility
• _psi_hbr, 4π×the real part of the Fourier component Ψ-H of the electric susceptibility
• _psi_hbi, 4π×the imaginary part of the Fourier component Ψ -H of the electric susceptibility
• _tc, the crystal thickness t
• _ang_as, the Bragg reflection’s asymmetry angle α. Note that α > 0 reduces the angle of the incident beam with the crystal surface, while α < 0 reduces the angle of the diffracted beam with the crystal surface. If α = 0, the incident and diffracted beams make the same angle with the crystal surface.
• _dTh, the deviation δθ of the incidence angle from the center of the rocking curve.
For the tests of this article, the electric susceptibilities were calculated by the Web application code X0h.13
Only cases in which the incident and diffracted beams lie in a common plane perpendicular to the crystal surface are calculated here. Also, only horizontal and vertical beam deflection are considered as these cases permit the use of standard Fast Fourier Transform algorithms without interpolation (see §2.2). Nonetheless, the great majority of crystal monochromators installed at synchrotron beamlines can now be simulated.
A Fourier transform is performed on the wavefront. As each Fourier component is a plane wave, plane-wave dynamical diffraction theory is used to generate the corresponding component of the diffracted wavefront. The formulas used here for evaluating the complex amplitude DH and the wave vector kH of a diffracted plane wave from an incident plane wave of complex amplitude D0 and wave vector k0 are given in Zachariasen.5 An inverse Fourier transform is then applied over all the diffracted plane waves to produce the final wavefront.
The application of the discrete Fourier transform to the sampled points of the incident and diffracted wave-fronts requires some modifications for asymmetric reflections in order to convert correctly between the reciprocal and the real representation. The following set of orthogonal unit vectors is defined as fixed to the crystal:
Another set of orthogonal unit vectors and is fixed to the incident beam’s central wave vector ko such that , where λo is the wavelength. Note that the incident wavefront is assumed to be monochromatic. These three vectors define the laboratory coordinate system. If the crystal is oriented so that = , then = ŝ and ŷ = . From here the crystal would be rotated about to set it to the Bragg reflection. Roll rotations (about ) and yaw rotations (about ŷ) could also subsequently be applied. Then the components of ŝ, and in the laboratory coordinate system are given by ŝ = (sx, sy, sz), = (nx,ny,nz) and = (tx,ty,tz), and the components of a general vector A can be transformed from the crystal frame to the laboratory frame and vice versa by an orthogonal transformation matrix:
Because diffraction from crystals deflects X-ray beams through a large angle, it is convenient to choose a new laboratory coordinate system for the diffracted beam. The orthogonal unit vectors x’ and y’ and z’ are defined starting from a form of Snell’s Law, kHT = koT + hT, where
are the components of ko, h and kH parallel to the crystal surface. Because |kH| = 1/λ0,
Let ’ = λokH so that it is parallel to the diffracted beam. The remaining vectors are given by ’ = h x ’/|h x ’| and ŷ’ = ’ x ’. Equation 1 can be used to calculate the components of these vectors in the incident beam laboratory system, or an analogous equation
can be used to calculate their components in the diffracted beam coordinate system.
Now, without loss of generality, each of the two independent orthogonally polarized components of the frequency-domain transverse electric field of the incident wavefront may be treated separately as a scalar. Let Uo(r) be one of these components given in the (x, y, z) coordinate system. The wavefront will then be represented in reciprocal space by the Fourier transform:
Rcryst(K0) is the crystal’s complex reflectivity function as given by Zachariasen5 for a plane wave of wave vector K0 having the polarization of U0(r). KH is the diffracted beam wave vector determined from K0 by Snell’s Law. r’ = x’’ + y’ ŷ’. The transformation from K0 to KH is simplified here by assuming that , , so that a matrix RH0 can be defined as follows:
then to a good approximation the elements of RH0 are
Thus the diffracted wavefront in real space is calculated by an inverse Fourier transform of the product of the plane-wave components of the incident wavefront with the respective complex reflectivities of the crystal. However, although the usual rectangular mesh of sampled points r of the incident wavefront will be transformed by Eq. 8 into a rectangular mesh of points (K0x, K0y) in reciprocal space, Eq. 13 will in general create an oblique mesh of points (x’,y’) that is not suitable for standard 2-D Fast Fourier Transform algorithms. The real-space mesh of the diffracted wavefront will be rectangular if a12 = a21 = 0 while a11 ≠ 0 and a22 ≠ 0; that is, if the deflection of the beam is either vertical (in the yz-plane) or horizontal (in the xz-plane). For vertical deflection, a11 = 1 and a22 = −|b|, where b is the asymmetry factor. For horizontal deflection, a11 = −|b| and a22 = 1. We note that the diffracted wavefront will also be rectangular if a11 = a22 = 0 while a12 ≠ 0 and a21 ≠ 0, corresponding to a 90° coordinate rotation that exchanges x’ and y’. Otherwise, the oblique mesh must be interpolated onto a rectangular mesh.
All crystals are silicon, the most common material for X-ray monochromators. The thickness t of the crystal was set to 10 mm. Because this is many times larger than the reflections’ extinction lengths, there is no significant energy flow from the bottom surface of the crystal.
Collimated Gaussian beam: Rocking curve comparison with XOP
Incident wave properties
First, a linearly horizontally polarized Gaussian source of 600 µm FWHM width in both horizontal (x) and vertical (y) directions is modeled. The source produces monochromatic X-rays of 9131 eV. The Rayleigh length zR of a Gaussian beam is given by
where λ is the wavelength of the radiation and w0 is the radius of the Gaussian wave at its waist, here assumed to be located at the source. The radius w(z) at any distance z along the beam axis from the waist is defined here as the distance from the axis at which the intensity drops to 1/e2 its central value. w0 is thus the FWHM at the waist divided by , or double the sigma value. Therefore, w0 = 509.6 µm and zR = 6008 m. The radius of curvature R of the wavefront is given by
and the beam radius at arbitrary distances by
If the crystal is placed 19 m downstream from the source, then R = 1.900 x 106 m and w is essentially the same as w0. Therefore, this source produces a highly collimated beam that can be used to test the crystal’s calculations for the plane-wave case.
SRW rocking curves: calculation method
The incident Gaussian wavefront at the crystal 19 m downstream from the source was calculated on a grid of points over a region twice as large as the FWHM of the Gaussian beam, which is large enough to ensure that the incident beam intensity is negligible at the edges and therefore does not cause artefacts under Fourier transformation. The number Nx x Ny of points on the grid was determined automatically by SRW. The intensities at all points on the grid were summed to obtain the total incident flux Tinc. Then the crystal propagator was applied to this wavefront, and the resulting diffracted wavefront was propagated through a 5 m drift space without any semi-analytical treatment of the leading quadratic terms of the phase. (Attempts to treat these phase terms semi-analytically created artefacts in the phase of the diffracted wavefront.) The number of points on the diffracted wave mesh remained the same as for the incident wave mesh, but the range of x or y was rescaled appropriately by the factor 1/|b|, where b is the asymmetry factor of the Bragg reflection. The intensities at all points on the grid of the final wavefront were summed to obtain the total diffracted flux Tdiff, from which SRWs value for the reflectivity, RSRW, was determined:
By stepping through a set of deviations δθ from the center of the rocking curve, it was possible to have SRW generate a rocking curve that could be compared to the results of XOP’s module for flat perfect crystals, XINPRO.10 However, while XOP assumes that the origin of δθ is the angle θB that strictly fulfills Bragg’s Law,
SRW sets the origin of _dTh at the true center of the rocking curve, which is shifted from θB because the diffracting crystal’s index of refraction is slightly different from the vacuum’s value of 1. The angular shift is given by Zachariasen:5
Therefore, in the figures of this section, Δθrefr was added to all values of _dTh in the SRW trials in order to shift the angular values to those of XOP. This serves as a check on SRWs ability to calculate the refraction shift correctly.
Comparison of rocking curves of SRW and XOP
Information on the diffraction cases for which rocking curves were calculated by SRW is displayed in Table 1. For X-rays of energy 9131 eV, silicon has Ψ0 = -1.173103 × 10-5 + 2.138459 × 10-7i and, because its crystal lattice is centrosymmetric, Ψ-H = ΨH.
Diffraction cases of silicon crystals for which rocking curves calculated by SRW (see § 3.2.2) and XOP10 were compared. The SRW rocking curves were calculated using a large (600 µm FWHM) collimated Gaussian beam that approximates a plane wave (see § 3.2.1). The XOP rocking curves were calculated using standard plane-wave dynamical diffraction theory. The source is linearly horizontally polarized and produces monochromatic radiation at 9131 eV.
|Reflection||(1 1 1)||(1 1 1)||(6 4 2)||(6 4 2)|
|Re ΨH||-6.201480 × 10-6||-6 201480 × 10-6||-3.112563 × 10-6||-3.112563 × 10-6|
|Im ΨH||+1.493016 × 10-7||+1 493016 × 10-7||+1.686663 × 10-7||+1.686663 × 10-7|
|Nx × Ny||2688 × 2688||2688 × 2688||2688 × 2688||2688 × 2688|
|Inc. range (mm)||-1.2 ≤ x ≤ +1.21.2 ≤y≤ +1.2||-1.2 ≤ x ≤ +1.2-1.2 ≤y≤ +1.2||-1.2 ≤ x ≤ +1.2-1.2 ≤y≤ +1.2||-1.2 ≤ x ≤ +1.2-1.2 ≤y≤ +1.2|
|Diff. range (mm)||-0.428 ≤ x’ ≤ +0.428- 1.2 ≤y’≤ +1.2||- 1.2 ≤ x’ ≤ +1.2-3.36 ≤ y’ ≤ +3.36||- 1.2 ≤ x’ ≤ +1.2- 14.6 ≤y’≤ +14.6||- 14.6 ≤ x’ ≤ +14.6- 1.2 ≤ y’ ≤ +1.2|
Figures 1 and 2 show the comparisons made between SRW and XOP for Si (1 1 1), which is used in most double crystal monochromators at X-ray synchrotron beamlines, and for Si (6 4 2), which is to be applied to a high-energy-resolution X-ray monochromator on the upcoming NSLS-II inelastic X-ray scattering beamline 10-ID. Note the good agreement between SRW and XOP in all four cases. That for the (6 4 2) reflection the angular width and peak reflectivity are lower in horizontal deflection than in vertical deflection is a result of the beam polarization relative to the “diffraction plane” that contains k0 and kH. Because the source is linearly horizontally polarized, vertical deflection (Fig. 2 (a)) results in p-polarization, in which the electric field points perpendicular to the diffraction plane, while horizontal deflection (Fig. 2(b)) results in p-polarization, in which the electric field lies within the diffraction plane. In the latter case, the strength of the Bragg reflection, which is given by ΨH, is effectively reduced by cos 2θB. Therefore the incident wave penetrates the crystal more deeply and suffers increased absorption.
Divergent Gaussian beam: real-space distributions of intensity and phase at various distances from crystal
Incident wave properties
Eqs. 14-16 in §3.2.1 are now applied to a Gaussian source of 2 µm FWHM width in both horizontal and vertical directions. The diffracting silicon crystal is still placed 19 m downstream from the source. The quantities that define the Gaussian beam are as follows:
This beam is strongly divergent and hence approximates a spherical wave. The linear horizontal polarization and the monochromatic spectrum at 9131 eV are as before.
To check SRW’s calculations, both the intensity and the phase of the Gaussian wavefront at the crystal position are analyzed. A contour plot of the wavefront’s intensity is shown in Fig. 3, and the profiles through the center of the wavefront along x and y are shown in Fig. 4. The FWHM of the best-fit Gaussian is 568.1 µm in both directions, matching the theoretical value very well. Figure 5 shows that the phase distribution through the center of the wavefront is quadratic in both x and y, as expected for a Gaussian beam. If the phase ϕ(x) is fit to a quadratic polynomial Aphx2 + Bphx + Cph, then the radius of the wavefront is π/(λAph) = 19.0005 m, again in good agreement with theory. The same quadratic fit applied in y gives the same curvature.
Bragg reflections examined with divergent Gaussian beam
The silicon crystal’s Bragg reflection is set to (6 4 2), whose spacing d and susceptibility coefficient ΨH have already been listed in Table 1. The examples are shown in Table 2. δθ = 0 to set the crystal orientation at the true center of the rocking curve.
Diffraction cases of silicon crystals for which diffraction of a divergent Gaussian beam was calculated. The source has a FWHM width of 2 µm in both horizontal and vertical directions, is linearly horizontally polarized, and produces monochromatic radiation at 9131 eV.
|Reflection||(6 4 2)||(6 4 2)|
|Darwin width (µrad)||9.89||34.5|
|Nx x Ny||6048 × 6048||6048 × 6048|
|Inc. range (mm)||- 1.8 ≤ x ≤ +1.8-1.8 ≤ y ≤ +1.8||- 1.8 ≤ x ≤ +1.8- 1.8 ≤ y ≤ +1.8|
|Diff. range (mm)||- 1.8 ≤ x ≤ +1.8- 1.8 ≤ y ≤ +1.8||- 1.8 ≤ x ≤ +1.8-21.91 ≤y≤ +21.91|
This case was chosen because its small Darwin width acts as an angular filter on the incident wavefront. The height of the beam after diffraction should be roughly (9.89 µrad)(19 m) = 0.188 mm, which is noticeably smaller than the incident wavefront’s height. The filtering is clear in Fig. 6, which shows the diffracted wavefront after propagation through drift spaces of various lengths.
Cross sections taken through the center of the diffracted wavefront provide further detail:
1. Because the deflection is vertical, the horizontal intensity profile is expected to be that of the unmodified incident Gaussian wavefront. Figure 7, which displays these profiles and their widths as a function of the length of the drift space after the crystal, bears this out. The width of the profile calculated by SRW agrees with that calculated using Eq. 16 except for a slight overall shift that is probably due to roundoff error.
2. The horizontal phase profile should also match that of the unmodified incident Gaussian wavefront. Figure 8 shows that this is the case. SRW calculates that the radius of the diffracted wavefront in x’ should be equal to the total distance along the beam axis from the source, as would be the case for a Gaussian wave at distances far beyond its Rayleigh length.
3. The vertical intensity profile will be cut off by the limited Darwin width, and because it will be affected by the angle-dependent reflectivity of the Bragg reflection, it will not be a Gaussian. This is shown in Fig. 9. However, note that the vertical FWHM of the diffracted wavefront is 0.177 mm, very close to the rough value of 0.188 mm calculated at the beginning of this section. The vertical size of the diffracted wavefront continues to grow linearly with drift space length.
4. For the vertical phase profile, SRW’s results are shown in Fig. 10. Notice that the radius of the wavefront is that of the unmodified Gaussian wave. This is a result of the absence of asymmetry for this case (α = 0). When the best-fit quadratic polynomial is subtracted from SRW’s calculated phase distribution after a 1 mm drift space, a phase jump of 2.82 rad emerges within the region of significant diffraction (Fig. 11). This is consistent with the results of plane-wave dynamical diffraction, which predict that the diffracted beam will shift in phase from 0 to π as the incidence angle is stepped through the Darwin curve.
This example will show how a wavefront is collimated by a highly asymmetric Bragg reflection. Because of the large asymmetry of this reflection, the vertical width of the beam that will be accepted by this reflection will be roughly (34.5 /xrad)(19 m) = 0.656 mm, which is greater than the vertical size of the incident wavefront. Figure 12 shows SRW’s calculation of the diffracted wavefront after propagation through drift spaces of various length. Note that SRW does indeed predict, in accordance with experiment, that a much larger portion of the beam will be diffracted. Also note that SRW’s diffracted wavefront spreads horizontally but not vertically with increasing distance from the crystal, indicating already that SRW predicts the expected collimation.
Cross sections taken through the center of the diffracted wavefront provide the following information:
1. The horizontal intensity profile is that of the unmodified incident Gaussian wavefront as shown in Fig. 13.
2. The horizontal phase profile is that of the unmodified incident Gaussian wavefront as shown in Fig. 14.
3. The vertical intensity profile is considerably broader than it is when α = 0. This is shown in Fig. 15. Notice that the collimation is manifested in the very slow increase in beam width with distance from the crystal.
4. The results for the vertical phase profile are shown in Fig. 16. After a drift space of 1 mm, the radius of the wavefront calculated by SRW is 2815.2 m. This is very close to the value 2823.6 m obtained by multiplying the 19 m radius of the incident Gaussian wave by 5-2. It results from the simultaneous stretching of the beam’s vertical cross section by a factor and reduction of the beam’s angular divergence by a factor , both of which are predicted by dynamical diffraction theory and observed in experiments. The phase jump found after a drift space of 1 mm when the best-fit quadratic polynomial is subtracted from the phase profile is much less here than at α = 0 because less of the Darwin curve is traversed across the beam’s width.
A new wavefront propagator for crystal diffraction has been integrated into SRW and tested under several different Bragg reflections, crystal orientations and asymmetry angles in conjunction with drift spaces. The incident wavefront was simulated as either a highly collimated Gaussian for comparisons with plane-wave dynamical diffraction or as a highly divergent Gaussian for comparisons with spherical-wave dynamical diffraction. The tests of this report have yielded results in line with theory.
Beyond the results reported here, the short-term priority is to examine further X-ray diffraction cases with fully and partially coherent radiation beams of XFEL and storage ring sources. In the medium term, more economical crystal propagators requiring lower sampling densities and less CPU time will be investigated. Future additions to the SRW package could include, first of all, a new crystal propagator module for the Laue case, which could be written with only slight modifications to the present propagator for the Bragg case. This would extend the functionality of SRW to X-ray phase plates and to crystal monochromators designed for high photon energies. In the longer term, support for diffraction from bent crystals could be added using already developed models such as Penning-Polder14 or Takagi-Taupin.15, 16
The authors thank Kawal Sawhney of Diamond Light Source for his constant support and David Laundy of Diamond Light Source for stimulating discussions. The present work was supported by the US Department of Energy, Contract No. DE-AC02-98CH10886.
Zachariasen, W. H., [Theory of X-Ray Diffraction in Crystals], John Wiley & Sons, New York (1945).Google Scholar
Sánchez del Río, M. and Dejus, R. J., “Status of XOP: an x-ray optics software tookit,” Proc. SPIE 5536, 171–174 (2004).Google Scholar
Chubar, O. and Elleaume, P., “Accurate and Efficient Computation of Synchrotron Radiation in the Near Field Region,” in [Proceedings of EPAC ‘98, Sixth European Particle Accelerator Conference], Liljeby, L., Myers, S., Petit-Jean-Genaz, C., Poole, J., and Rensfelt, K.-G., eds., 1177–1179 (1998).Google Scholar
Penning, P. and Polder, D., “Anomalous transmission of x-rays in elastically deformed crystals,” Philips Res. Rep. 16, 419–440 (1961).Google Scholar
Taupin, D., “Théorie dynamique de la diffraction des rayons X par les cristaux déformés,” Bull. Soc. Franç. Minér. Cryst. 87, 469–511 (1964).Google Scholar