Partially coherent X-ray wavefront propagation simulations including grazing-incidence focusing optics

X-ray beamlines in modern synchrotron radiation sources make extensive use of grazing-incidence reflective optics, in particular Kirkpatrick–Baez elliptical mirror systems. These systems can focus the incoming X-rays down to nanometer-scale spot sizes while maintaining relatively large acceptance apertures and high flux in the focused radiation spots. In low-emittance storage rings and in free-electron lasers such systems are used with partially or even nearly fully coherent X-ray beams and often target diffraction-limited resolution. Therefore, their accurate simulation and modeling has to be performed within the framework of wave optics. Here the implementation and benchmarking of a wave-optics method for the simulation of grazing-incidence mirrors based on the local stationary-phase approximation or, in other words, the local propagation of the radiation electric field along geometrical rays, is described. The proposed method is CPU-efficient and fully compatible with the numerical methods of Fourier optics. It has been implemented in the Synchrotron Radiation Workshop (SRW) computer code and extensively tested against the geometrical ray-tracing code SHADOW. The test simulations have been performed for cases without and with diffraction at mirror apertures, including cases where the grazing-incidence mirrors can be hardly approximated by ideal lenses. Good agreement between the SRW and SHADOW simulation results is observed in the cases without diffraction. The differences between the simulation results obtained by the two codes in diffraction-dominated cases for illumination with fully or partially coherent radiation are analyzed and interpreted. The application of the new method for the simulation of wavefront propagation through a high-resolution X-ray microspectroscopy beamline at the National Synchrotron Light Source II (Brookhaven National Laboratory, USA) is demonstrated.

1. Introduction

In modern third- and fourth-generation synchrotron radiation sources, i.e. storage rings and free-electron lasers (FELs), X-ray beam micro-focusing is routinely used in a large number of experimental techniques. The appearance of new low-emittance storage ring sources (Ozaki et al., 2007; Eriksson et al., 2013; Reich, 2013), which generate considerable amounts of coherent X-ray flux, the start-up of operation of the X-ray FEL, producing highly coherent X-ray pulses (Emma et al., 2010; Ishikawa et al., 2012), and the recent substantial improvement of mirror surface quality allowing nearly diffraction-limited resolution to be achieved (Mimura et al., 2010; Yamauchi et al., 2011), make it particularly important to use accurate wave-optics-based methods (Bahrdt, 1997; Chubar & Elleaume, 1998; Chubar et al., 2002) for the simulation, optimization and design of such mirror systems and entire beamlines containing them.

At low demagnification, when the longitudinal extent of a focusing mirror is much smaller than the distance from the mirror to the focused radiation spot, the focusing mirror can be reasonably approximated by an ideal thin lens (Bowler et al., 2008). This may not be the case for mirror systems performing strong focusing, when the distance to the focused spot becomes comparable with the mirror length. Therefore, a special numerical method had to be developed to enable the simulation of such cases within the framework of wave optics. In the pioneering simulations of wave propagation through a beamline containing grazing-incidence optics, Bahrdt (1997) used an approach combining the stationary phase approximation (see, for example, Born & Wolf, 1999) with a transport matrix formalism. This approach enabled the simulation of wave propagation through an entire optical system, when the system did not contain apertures limiting the beam and causing diffraction. Taking into account the diffraction with this approach was not straightforward, because the description of the diffraction requires a different type of asymptotic expansion, other than the stationary phase method. Chubar & Elleaume (1998) and Chubar et al. (2002) implemented a different approach to propagate the wave through an entire beamline, which is based on Fourier optics methods and employs the numerical apparatus of fast Fourier transforms (FFT). This method `natively' takes into account both wave aberrations and diffraction at apertures in a very CPU-efficient manner, and it easily allows high spatial frequencies to be taken into consideration (e.g. surface height profile or slope errors; De Andrade et al., 2011). However, the Fourier optics method, in its traditional implementation, operates only with `thin' optical elements (Goodman, 2005), i.e. it is not directly applicable for an accurate simulation of grazing-incidence optics with large longitudinal extent along the optical axis. In some cases this limitation can be overcome using the `split-operator' approach, allowing one to simulate an extended `thick' optical element by a sequence of `thin' elements and drift spaces (Hoekstra, 1997). However, this may slow down the simulation, compromising one of the main advantages of the Fourier optics method, its CPU efficiency.

In this paper we describe the implementation and extensive benchmarking of a `hybrid' method, which preserves all the advantages and overcomes the limitations of the Fourier optics method and the stationary phase approximation, profiting from their complementarity. Our method is still essentially Fourier optics based; it natively treats diffraction, and yet allows for simulation of extended grazing-incidence mirrors almost as efficiently as `thin' optical elements. The new method has been implemented in the Synchrotron Radiation Workshop (SRW) code (Chubar & Elleaume, 1998; Chubar et al., 2002).

The benchmarking simulations are performed for the Kirkpatrick–Baez (KB) type mirror system (Kirkpatrick & Baez, 1948) composed of two orthogonal elliptical cylinders. This popular type of X-ray optical mirror system is commonly used both in `high throughput' and in `high demagnification' modes. In the former, flux at a sample is favored over spatial resolution; in the latter, the flux in the focused spot is sacrificed in favor of the spatial resolution. Another attractive feature of this type of optics, that makes it particularly efficient for microspectroscopy-type experiments, is the absence of dispersion (i.e. the absence of strong dependence of the focusing properties on radiation photon energy, in contrast with diffractive and refractive optical elements).

The parameters of the optical schemes being considered are relatively close to those used for several beamlines of the National Synchrotron Light Source II (NSLS-II), the new low-emittance storage ring source constructed at Brookhaven National Laboratory, Upton, NY, USA. The first of the two beamlines uses the KB mirror system in the `high throughput' mode, and the second in the `high resolution' mode. The results of the wave-optics simulations performed for the fully and partially coherent radiation, using the new method implemented in the SRW code, are compared with the results obtained simulating the elliptical cylinders by `ideal lenses' (still using SRW) and with geometrical ray-tracing results obtained using the SHADOW code (Sanchez del Rio et al., 2011) for the same parameters of the source and optical schemes. The general agreement and some differences in the results, associated with the wave-optical phenomena, are discussed and interpreted.

2. Combining the Fourier optics and the stationary phase methods

The propagation of the transverse components of the frequency-domain electric field of radiation in free space or through an optical element can be described by an integral operator, i.e. `propagator', which, for the majority of practically important cases, can be formally represented by

where is a kernel specific to the jth optical element (j = 1, 2,…). In the general case, is a matrix (tensor) with its elements being dependent on the frequency ω and on the transverse coordinates in a plane before (x_j-1,y_j-1) and a plane after (x_j,y_j) the optical element. and are the transverse electric field components in the planes before and after the optical element, which, besides the transverse coordinates and frequency, depend also on a set of phase-space coordinates of a source . For synchrotron radiation emitted by one relativistic electron moving in an external (in the general case non-uniform) magnetic field, this is a set of the electron coordinates in six-dimensional phase space of the electron beam (see, for example, Chubar et al., 2011): = , where are the electron initial transverse coordinates and angles, is the electron energy, z_e is the initial longitudinal position. We note that usually are small and that the z_e dependence can often be omitted at simulations of spontaneous emission in storage rings.

Following the Huygens–Fresnel principle, the kernel of equation (1) for the propagation in free space between two parallel planes in the paraxial approximation is (see, for example, Born & Wolf, 1999)

where L is the distance between the planes, c is the speed of light and I is the unit matrix, i.e. in this case equation (1) is reduced to a scalar convolution-type relation for each of the two transverse components of the electric field. This allows the convolution theorem and the digital FFT to be applied for its fast numerical calculation. Often, the square-root in the argument of the exponent in equation (2) can be approximated by the first two terms of its Taylor expansion, assuming L , L ,

This approximation is convenient for various numerical manipulations aiming to increase the efficiency and to reduce the memory requirement for the calculation of the integrals in equation (1). One of such manipulations is the analytical treatment of the quadratic phase term, allowing for the free-space propagation simulations with a reduced sampling rate (Chubar et al., 2008), which dramatically improves technical feasibility and robustness of the Fourier optics method. Furthermore, if the dependence on x_j-1²,y_j-1² in the phases of components of the product can be ignored, as usually happens at propagation over a large distance L from or to a waist, then

and equation (1) reduces to a simple Fourier transform.

For a `thin' optical element, the kernel can be represented by

where is a complex transmission (matrix) function and is the δ-function. Equation (3) reduces equation (1) to a simple product. For an isotropic `thin' optical element, = , where is a scalar complex transmission function.

We should point out that equations (1)–(3) not only represent a good approximation for the simulation of coherent radiation propagation through an optical system but they also allow for accurate description of processes taking place in a large number of experiments, such as coherent diffraction imaging, coherent scattering, interferometry and others. Hence a robust numerical implementation of this Fourier-optics based `apparatus' can find a very large number of applications, well beyond the accurate simulation of radiation propagation through X-ray beamline optics.

Equation (1) when combined with (2) [or (2′)] and (3) can also be used for simulation of radiation propagation through an optical element with a large extent along the optical axis using the `split-operator' or `beam propagation' technique (Hoekstra, 1997). However, the numerical efficiency of such a method applied to two-dimensional wavefronts with a large number of steps versus longitudinal position may be not very high.

Alternatively, for such longitudinally extended optical elements, and in particular for grazing-incidence mirrors, one can use a numerical procedure which can be formally described by equation (1) with the following kernel,

where is a matrix function defining local transformations of amplitudes of electric field components (e.g. at reflection from a mirror surface and propagation between input and output planes within the mirror extent), are scalar functions defining the transformation of coordinates for points in transverse planes before and after the optical element, and is a scalar function defining the optical path between the points in the input and output planes. These functions can be found using the stationary phase approximation applied to the Kirchhoff integral over the mirror surface, or simply by applying locally the laws of geometrical optics and satisfying boundary conditions for the electric field components.

It is usually straightforward to find numerically for a given point in the input plane (x_j-1,y_j-1) the corresponding one in the output plane (x_j,y_j). This can be done by `sending a ray' from the point (x_j-1,y_j-1) (before the mirror) in the direction provided by the local gradient of the input radiation field phase, finding the intersection with the optical element (e.g. mirror surface), generating the reflected ray, and finding its intersection with the output plane (after the mirror).

The local `ray-tracing' does not guarantee that the points in the output plane will form a regular rectangular mesh, which is a requirement for using FFT to propagate the electric field from that plane further on. Hence, as a final step of the local propagation, the electric field components are re-calculated on the required rectangular mesh, using a two-dimensional interpolation.

The approach just described allows an electric field to be propagated from a transverse plane before an optical element to a transverse plane after it. The process can be reiterated by applying a sequence of propagators, according to equation (1), with the kernels corresponding to drift spaces and individual optical elements [see equations (2)–(4)]. In many cases these propagators can have CPU-efficient implementations based on FFT and simple multiplications. Therefore the entire simulation, even for a wide wavefront sampled on a two-dimensional mesh with thousands of points in both transverse directions, may take no longer than a few seconds on a moderately equipped desktop computer. In the cases of high or full transverse coherence, such as for the emission of UV and X-ray FELs, or spontaneous infrared emission in third-generation sources, the simulation procedure does not require any other processing. Additional information such as the spectral flux per unit surface, polarization and other radiation characteristics can be extracted from the propagated electric field components. However, in storage rings and in energy-recovery linacs, the synchrotron emission by the entire electron beam is known to be only partially coherent for the most part of the spectrum. For its accurate simulation, one has to take into account that the radiation generated by electrons having different phase-space coordinates is not temporally coherent for the most part of the spectrum, and that the wavefronts produced by these electrons may propagate differently through the optical elements of a beamline (Chubar et al., 2011), making the simulation more complex and CPU-intensive.

The simulation of both the fully and partially coherent radiation propagation through beamline optical systems using the approximations of the `thin' and extended optical elements is possible with the recent Open Source versions of the SRW code (https://github.com/ochubar/SRW ).

3. Simulation results

In order to verify the validity and reliability of the new wavefront propagation method implemented for grazing-incidence mirrors, we benchmarked it against two other well known methods: the Fourier optics, that simulates accurately the wave-optical effects (such as the diffraction at the apertures) however uses the `ideal thin lens' approximation for simulating the grazing-incidence mirrors, and the geometrical ray-tracing, that models accurately the reflection from an extended mirror surface, yet does not take into account the wave-optical effects. The three methods were compared by simulating the propagation of equivalently defined X-ray beams through beamlines having the optical schemes shown in Figs. 1 and 6. In §3.1 and §3.2 we describe the simulation results obtained for fully and partially coherent (at a low degree of coherence) Gaussian beams, and in §3.3 for a partially coherent undulator radiation source in a low-β straight section of NSLS-II. The simulations were performed using the recently implemented Python interfaces of the SRW and SHADOW codes. This avoided redefining the parameters of the sources and the optical schemes used in the simulations.

3.1. Fully coherent Gaussian beam

Before discussing the simulations performed for a coherent Gaussian beam, let us first note that the use of the geometrical ray-tracing method for such simulations requires special considerations about the parameters of a source and/or an input radiation beam. The geometrical ray-tracing procedure does not `natively' take into account the laws of wave optics, and formally allows sources or beams to have any size and divergence. Therefore the existing physical constraints on possible source parameters have to be imposed `artificially' before starting the ray-tracing simulations. In particular, the well known fundamental relation between the RMS source size and divergence of a coherent Gaussian radiation beam,

where λ is the radiation wavelength, should be taken into account when defining the Gaussian source for the ray-tracing simulations. Once this condition is met, the ray-tracing method, which is known to preserve the phase-space volume of the propagating beam, can in principle be used for simulating the propagation even of a coherent beam, at least for optical schemes without apertures at which diffraction may occur.

Note that, for the wavefront propagation simulation method, it is sufficient to define, for example, only the size of a coherent Gaussian source at its waist, and the correct divergence, in complete agreement with equation (5), will `develop' automatically in the process of simulation, since it is imposed by equations (1), (2′) and (2′′). Of course, the numerical wavefront propagation method allows for using any distribution of input electric field in equation (1), not just the Gaussian one. On the other hand, since for a non-Gaussian field distribution equation (5) converts to the inequality , the use of the geometrical ray-tracing with the RMS source sizes and beam divergences related by equation (5) may lead to incorrect results.

The simulations presented in this section were performed for a coherent monochromatic Gaussian beam source at 12 keV photon energy, with the horizontal and vertical RMS waist sizes equal to ∼1.84 µm and the RMS angular divergences equal to ∼4.47 µrad, in accord with equation (5). The optical scheme, for which the first set of simulations has been performed, is illustrated in Fig. 1. The scheme includes two mirrors focusing in the vertical and in the horizontal planes (VKB and HKB, respectively). The mirror surfaces are elliptical cylinders; their length in the tangential direction is 0.5 m; the incidence angle at their centers is 3.5 mrad; the ellipses' semi-axes lengths are equal to a = 30.5 m and b = 27.1 mm for VKB and a = 30.5 m and b = 19.25 mm for HKB. The VKB is located at 60.0 m from the source, the HKB at 60.5 m from the source; and the system focuses the beam at 61.0 m from the source. A rectangular aperture can be optionally introduced in front of the VKB to test the impact of diffraction.

Figure 1 Optical scheme used for simulations with fully coherent Gaussian beam.

The calculations by the three different methods were initially performed without the aperture. The distributions of intensity (i.e. flux per unit surface) obtained at three different locations, i.e. the focal point of the elliptical mirrors, 5 cm upstream and 5 cm downstream from the focal point, are presented in Fig. 2. As one can see from Fig. 2(a), the distributions obtained by the different methods are almost indistinguishable between each other at the location of the focus. The source, optical scheme and the focal spot parameters used and/or obtained in this and other simulations discussed in this paper are summarized in Table 1. As seen in the table, the ratio between the full width at half-maximum (FWHM) and the root mean square (RMS) focal spot sizes is close to 2.35 in the three cases as expected for a Gaussian beam. The horizontal and vertical RMS (, ) and skewness (, ) values specified in Table 1 were estimated numerically from the radiation intensity distributions in the focal spots F(x,y) normalized to 1,

Table 1
Basic parameters of the sources, the optical schemes and the focal spots obtained using different simulation methods for the cases described in §3.1–§3.3

Source				Focal spot (calculated by three methods†)
RMS size, H/V (µm)	RMS divergence, H/V (µrad)	Photon energy (keV)	De-magnification, H/V	Method #	FWHM size‡, H/V (nm)	FWHM to RMS§ ratio, H/V	Skewness§, H/V
§3.1, Fig. 2(a): fully coherent Gaussian beam at lower demagnification without apertures
1.84/1.84	4.47/4.47	12	121/60	1	37.7/75.4	2.27/2.43	0.0/0.0
				2	38.7/78.7	2.37/2.48	0.0/0.0
				3	35.0/73.0	2.35/2.43	0.0/0.0
§3.1, Figs. 3(a) and 3(b): fully coherent Gaussian beam at lower demagnification with 0.5 mm × 0.5 mm aperture
1.84/1.84	4.47/4.47	12	121/60	1	92.6/179	1.53/1.48	0.0/0.0
				2	92.6/179	1.53/1.48	0.0/0.0
				3	36.2/73.3	2.45/2.45	0.0/0.0
§3.1, Fig. 4: fully coherent Gaussian beam at higher demagnification without apertures
10/10	0.822/0.822	12	546/233	1	41.5/101	2.29/2.34	0.0/0.0
				2	42.9/101	2.34/2.35	0.0/0.0
				3	42.4/101	2.25/2.31	0.0/0.0
60/60	0.137/0.137	12	546/233	1	244/600	2.30/2.42	−0.3/−0.1
				2	257/602	2.37/2.33	0.0/0.0
				3	246/603	2.36/2.44	−0.3/−0.1
§3.2, Fig. 5(a): low-coherence Gaussian beam at higher demagnification
137/5.53	7.95/4.73	12	546/233	1	445/67.8	1.80/1.43	−0.1/0.0
				2	566/68.9	2.35/1.34	0.0/0.0
				3	457/66.9	1.76/2.32	−0.1/0.0
§3.3, Fig. 7: high-resolution microspectroscopy beamline with horizontal secondary source aperture of 50 µm
33.5/∼4.7	18.3/∼8.5	7.18	378/346	1	187/121	1.27/1.44	0.0/0.0
				2	196/126	1.32/1.46	−0.1/0.0
				3	171/109	2.50/2.27	0.0/−0.1
§3.3, Fig. 7: high-resolution microspectroscopy beamline with horizontal secondary source aperture of 20 µm
33.5/∼4.7	18.3/∼8.5	7.18	378/346	1	91.6/120	0.380/1.45	0.0/0.0
				2	95.3/126	0.418/1.46	−0.1/0.0
				3	88.0/108	2.84/2.27	0.0/−0.1
§3.3, Fig. 7: high-resolution microspectroscopy beamline with horizontal secondary source aperture of 10 µm
33.5/∼4.7	18.3/∼8.5	7.18	378/346	1	68.3/118	0.170/1.38	−0.1/0.0
				2	70.8/123	0.269/1.40	−0.1/0.0
				3	45.4/114	2.94/2.27	0.0/−0.1
§3.3, Fig. 7: high-resolution microspectroscopy beamline with horizontal secondary source aperture of 5 µm
33.5/∼4.7	18.3/∼8.5	7.18	378/346	1	64.3/118	0.134/1.33	−0.1/0.0
				2	67.4/122	0.129/1.49	−0.1/0.0
				3	23.2/105	2.92/2.27	0.0/0.0

†The three calculation methods used are: wavefront propagation with extended grazing-incident mirror models (method #1); wavefront propagation with mirrors simulated by ideal lenses (method #2); geometrical ray-tracing with extended grazing-incident mirror models (method #3).
‡The estimated error of the FWHM values for different methods varies between ∼0.1 and several percent.
§The RMS and skewness values were estimated numerically as described by equation (6); the relative error of the RMS values is ∼5–10%; the relative error of the skewness values is larger, especially in the cases when the absolute skewness values are small.

Figure 2 Horizontal (left) and vertical (right) cuts (along the horizontal and vertical planes containing the optical axis) of the intensity distributions obtained by three different calculation methods: wavefront propagation with grazing-incidence elliptical mirrors extended along optical axis (1), wavefront propagation with the focusing mirrors simulated by ideal thin lenses (2), and geometrical ray-tracing with the grazing-incidence elliptical mirrors (3): (a) at the longitudinal position of the focal point; (b) 5 cm upstream from the focal point; (c) 5 cm downstream from the focal point.

where

and are chosen so that

These constraints have been used to ensure the convergence of the integrals in equation (6) even for the intensity distributions strongly dominated by diffraction (see below).

A slight difference between the distributions develops as the observation plane moves off the focal point [see Figs. 2(b) and 2(c)]: the distributions obtained for the extended elliptical mirror, using both the wavefront propagation (curve 1) and the geometrical ray-tracing (curve 3) methods become slightly asymmetric, while staying in good agreement with each other, whereas the distribution obtained using the ideal thin lens approximation (curve 2) remains symmetric with respect to the optical axis. It is noteworthy that the asymmetry is stronger in the horizontal plane, which is explained by the stronger focusing taking place in that plane: the demagnification in the horizontal plane is 121, whereas in the vertical plane it is 60. As Figs. 2(b) and 2(c) show, the asymmetry of the distributions is different before and after the focal point.

The agreement between the wavefront propagation and the geometrical ray-tracing results ceases when the fully coherent Gaussian beam is truncated by an aperture, such as the one located before the VKB, as shown in Fig. 1. Fig. 3 shows the horizontal and vertical cuts of the intensity distributions at and after the focus, obtained using the three calculation methods with the aperture closed down to 0.5 mm in the horizontal and vertical directions. We note that this aperture size is smaller than the size of both KB mirrors projected on the transverse plane (1.75 mm) and also smaller than the FWHM size of the radiation beam at the location of the mirrors (∼0.63 mm).

Figure 3 Horizontal (left) and vertical (right) cuts of intensity distributions by the horizontal and vertical planes containing the optical axis, obtained by the three different calculation methods as described in the caption to Fig. 2: (a) at the longitudinal position of the focal point; (b) close-up of the same distribution at the focal point; (c) 5 cm downstream from the focal point.

As one can see from Fig. 3, the wavefront propagation method, both with the extended mirror (curve 1) and with the thin lens approximation (curve 2), gives considerably larger distribution sizes (by a factor of ∼2.5 for FWHM values) and lower intensity peaks (by a factor of ∼8) than the geometrical ray-tracing method (curve 3). The FWHM and RMS sizes of the focal spot obtained by the different methods are given in Table 1. We note that the curves 1 and 2 show secondary maxima (see Fig. 3b), revealing the wave-optics (diffraction) origin of the effects under consideration. Obviously, these effects cannot be simulated by pure geometrical ray-tracing, which ignores the diffraction occurring at the aperture located before the mirrors. The diffraction fringes are clearly present in the intensity distributions obtained by the wavefront propagation at the observation plane located 5 cm downstream from the focal point [see curves 1 and 2 in Fig. 3(c)], even though in this case the sizes of all distributions are comparable. One can also notice in Fig. 3(c) a small asymmetry of the distributions obtained for the extended mirrors, as was the case in the simulations without the aperture.

The next series of simulations aimed to explore how the intensity distributions in the focal spot vary when increasing the demagnification of the KB system and the source size of a coherent Gaussian beam. The parameters of the optical system were changed as follows: the KB mirror lengths were reduced to 15 cm, the VKB center was relocated at ∼61.238 m, the HKB at 61.388 m and the focus at 61.5 m from the source. The photon energy was kept at 12 keV but the RMS waist sizes of the round Gaussian beam source were set to 10 µm and 60 µm (for two different cases of simulation). The corresponding RMS angular divergences, following equation (5), were respectively ∼0.822 µrad and ∼0.137 µrad. Let us note that such cases of low-divergence coherent hard X-ray beams are not purely academic: these parameters are not far from those of radiation beams in existing X-ray FELs (Emma et al., 2010; Ishikawa et al., 2012).

The results of the simulated intensity at the focus are presented in Fig. 4. As one can see, the three different methods give very close results for the 10 µm source size and ∼0.822 µrad divergence Gaussian beam (Fig. 4a). However, the situation changes dramatically for the larger, 60 µm, source size and the lower, ∼0.137 µrad, beam divergence [see Figs. 4(b) and 4(c)]. The intensity distributions at the geometrical focal point, obtained using the extended grazing-incidence mirror methods, i.e. both the wavefront propagation and the geometrical ray-tracing, become considerably asymmetric, with the asymmetry being more significant in the horizontal plane, where the demagnification is stronger (∼546 against ∼233 for the vertical plane). The remarkable agreement between the wavefront propagation and the ray-tracing simulations using the extended mirror models in that case is also illustrated by the two-dimensional image plots in Fig. 4(c). Such good agreement is not surprising since there is no diffraction at the mirror apertures (the beam is not truncated by any mirror) and the source size and divergence of the beam for the ray-tracing simulations were chosen to obey the wave-optics relation expressed in equation (5). The asymmetry of the obtained intensity distributions is characterized by large absolute values of the skewness parameter in Table 1.

Figure 4 (a) Horizontal (left) and vertical (right) cuts of intensity distributions at the geometrical focal point of the KB system with ∼546 (∼233) horizontal (vertical) demagnification, obtained by the three different calculation methods, as described in the caption to Fig. 2, for a coherent Gaussian beam with 10 µm RMS waist size (at the source) and ∼0.822 µrad divergence. (b) Same for the Gaussian beam with 60 µm RMS waist size and ∼0.137 µrad divergence. (c) Image plots of the two-dimensional intensity distributions obtained in (b) using the extended grazing-incidence mirror methods: the wavefront propagation (left) and the geometrical ray-tracing (right).

We should point out that the significant deviation of the focal spot intensity distributions obtained for the grazing-incidence elliptical mirrors from the distribution obtained using the ideal thin lens approximation suggests that the elliptical mirror shape may not be the optimal one to focus a beam with such a large source size and low divergence. In fact, it can be shown that when the contribution to the size of the beam at the mirror from the divergence is much smaller than the source size, the appropriate focusing optical element to use is a parabolic mirror since the beam can be better described as collimated.

A special comment can be made about the CPU efficiency of calculations using the wavefront propagation and the geometrical ray-tracing methods in the case of a coherent radiation beam. The geometrical ray-tracing method, which uses the stochastic Monte Carlo approach for seeding the initial rays, requires huge statistics to achieve high-quality intensity profiles by slicing/binning a distribution along a line. The number of rays used in each of the ray-tracing simulations described above was 25 million. The corresponding ray-tracing computation time, which scales linearly with the number of rays, was about 10 min on a single core, whereas with the wavefront propagation method exploiting the two-dimensional FFT, with 1000–2000 equidistant points used for each transverse direction, the calculation on one core usually lasted about 10 s. Besides, no significant increase of the calculation time was observed in the case of the wavefront propagation simulations with the extended grazing-incidence mirrors compared to those with ideal lenses. This clearly shows that, for optical simulations with fully coherent radiation beams, the FFT-based wavefront propagation method (as implemented in the SRW code) is not only more accurate but also numerically more efficient than the geometrical ray-tracing (as implemented in the SHADOW code).

More generally, one can note that the numerical complexity of FFT scales as O[Nlog(N)], with N being the total number of data points (equal to the product of numbers of points in each in the two directions in the case of two-dimensional FFT). The complexity of other operations used in the Fourier optics method besides the FFT [e.g. those described by equation (1) with kernels given by equations (3) and (4)], scales approximately as O(N). The complexity of the geometrical ray-tracing method is formally also , where is the total number of rays used. However, whereas in the Fourier optics method each of the N data points is located on a regular rectangular mesh and directly carries information about the radiation electric field and so about the flux per unit surface area, in the Monte Carlo driven geometrical ray-tracing the knowledge of the density of ray intersection points with a surface is required for estimating the radiation flux per unit area of that surface. Therefore, for a comparable numerical accuracy between the two methods one has to use a number of rays for obtaining a sufficiently accurate value of the flux per unit surface at one point, so that the required total number of rays has to considerably exceed the number of final data points: N. This explains the smaller CPU time that was required for SRW simulation compared with the SHADOW simulation in the case of the fully coherent radiation beam considered in this section.

3.2. Low-coherence Gaussian beam

In this section, we present the wavefront propagation and geometrical ray-tracing simulation results obtained for a partially coherent X-ray beam focused by the KB mirror system used in the simulations illustrated in Fig. 4.

The partially coherent X-ray beam has been simulated using the Gaussian Schell model (Gori & Palma, 1978), allowing for approximate description of partially coherent undulator radiation in a synchrotron source (Coïsson, 1995). In this model the fully coherent single-electron undulator radiation is approximated by a coherent Gaussian beam, which, strictly speaking, is not fully accurate and results in a limited applicability of this model (Geloni et al., 2008); and the electron beam, in which different electrons are producing non-correlated temporally incoherent emission, is represented by an incoherent source with a Gaussian density distribution, the latter being a good approximation for the equilibrium state of an electron beam circulating in a storage ring.

The wave-optics-based simulation of the propagation of such a beam through an optical system is simulated in SRW code by propagating multiple uncorrelated coherent beams `originating' from different phase-space `points' of an extended incoherent source and summing their intensity distributions at the final plane. For the ray-tracing simulations, such beam is simulated simply by a Gaussian beam with the horizontal and vertical RMS sizes and divergences given by

i.e. resulting from the convolution of the coherent Gaussian beam with the RMS size and divergence [related by equation (5)], and an incoherent Gaussian beam with the horizontal and vertical RMS sizes and divergences . In the simulations, and values were taken to correspond to the central cone of single-electron radiation from a ∼2.6 m-long undulator at 12 keV photon energy ( ≃ 1.84 µm, ≃ 4.47 µrad); and , to correspond to the horizontal and vertical RMS sizes and divergences of electron beam in the middle of a high-β straight section of NSLS-II ( = 137.0 µm, = 5.22 µm, = 6.57 µrad, = 1.53 µrad). The horizontal and vertical FWHM dimensions of the radiation beam at the location of the KB mirrors are ∼1.44 mm and ∼0.68 mm, respectively, whereas the projected mirror size is ∼0.53 mm. The corresponding transverse coherence lengths are ∼8.5 µm and ∼0.22 mm, i.e. the X-ray beam is more coherent in the vertical direction. The transverse coherence lengths were estimated as separation distances between two slits in Young's interference scheme resulting in the visibility of fringes equal to 0.5.

The simulation results for this partially coherent X-ray beam, obtained by the three different methods described above, are presented in Figs. 5(a)–5(c). Fig. 5(a) shows the horizontal and vertical cuts of intensity distributions calculated at the geometrical focal point. One can notice the deviation from the Gaussian shape when simulations are performed with the elliptical mirror shapes, in particular in the horizontal direction (Fig. 5a, left graph). This deviation is confirmed by the ratios of the FWHM to RMS size values: for the distributions obtained using both the wavefront propagation and the geometrical ray-tracing, this ratio is ∼1.8, whereas for the Gaussian distribution one would expect it to be ∼2.35, as approximately obtained by the wavefront propagation when the mirrors are simulated by ideal thin lenses (see Table 1). This deviation of the intensity distribution shape from the Gaussian has a straightforward qualitative interpretation. By definition, an elliptical mirror images a point source located in one of its foci, to a point in the other focus. The demagnification of a source with a finite transverse size is formally varying along the mirror length: for the longitudinal positions between the upstream mirror edge and the mirror center the demagnification is lower, and for the positions between the mirror center and downstream edge it is higher than for the ideal thin lens located at the mirror center. If the mirror length is much smaller than the distances from the mirror center to the source and to the image, this effect is negligible. However, if this is not the case, the demagnification variation along the mirror produces at the final focal point a more `peaky' intensity distribution, with wider `tails' compared with the initial Gaussian distribution. The radiation reflected by the downstream part of the mirror contributes to the narrower peak, whereas the radiation reflected by the upstream part contributes more to the `tails'. Or in other words, the downstream half of the mirror has a higher numerical aperture than the upstream half.

Figure 5 Horizontal (left) and vertical (right) cuts (along the horizontal and vertical plane containing the optical axis) of the intensity distributions of partially coherent X-ray beam at 12 keV photon energy, obtained by three calculation methods as described in the caption to Fig. 2: (a) at the longitudinal position of the focal point; (b) 2 cm upstream from the focal point; (c) 2 cm downstream from the focal point.

Note that the vertical FWHM to RMS size ratio of the intensity distributions at the focus obtained by the wavefront propagation (both using the extended mirror and the ideal thin lens approximations) is ∼1.4, i.e. also less than 2.35, whereas for the ray-tracing simulations it is ∼2.3 (see Table 1). Such small values of the ratio obtained from the wavefront propagation can be explained by the contribution of diffraction of the partially coherent radiation, which is important in the vertical direction and tends to increase the RMS size values of the intensity distributions.

The effects of the radiation diffraction at the aperture (located before the VKB and comparable in size with the projected VKB and HKB mirror dimensions) can also be observed in the intensity distributions in Figs. 5(b) and 5(c), which were calculated for the longitudinal positions 2 cm before and 2 cm after the geometrical focus. Diffraction patterns with multiple fringes are clearly seen in the vertical cuts of the distributions (in the graphs on the right), and even in the horizontal cuts (graphs on the left) at the edges of the distributions, despite the fact that the coherence is very low along the horizontal direction. The distributions obtained for the elliptical mirrors before and after the focal point are also asymmetric and specular to each other, as was the case for the fully coherent radiation beam [see Figs. 2(b) and 2(c) and Fig. 3(c)]. The asymmetrical shapes of the `out of focus' intensity distributions are very similar in the wavefront propagation with the elliptical mirrors and the geometrical ray-tracing simulations, even though the diffraction patterns are of course missing in the case of the ray-tracing. The structures observed in the intensity distribution cuts of the ray-tracing are due to statistical noise.

The comparison of the efficiency of simulations by the different methods in the case of the low-coherence radiation beam is much less favorable for the wavefront propagation method. In order to obtain smooth and stable intensity distributions, especially at the focal point where the contribution of the incoherent source size is quite significant, one had to average intensity distributions of propagated wavefronts from a very large number, up to a million, of coherent Gaussian sources, non-correlated between each other. These sources had Gaussian distributions over the four-dimensional phase space of the incoherent source (with the phase space coordinates being horizontal and vertical positions and angles). The calculation was performed on 64 nodes of a computer cluster, using a parallel version of the SRW code, which makes use of the message passing interface (MPI), and it lasted more than 6 h. The convergence was faster for the `out of focus' intensity distributions, where the effects related to the incoherent source size are smaller than at the focus. On the other hand, the geometrical ray-tracing simulation, which ignores radiation coherence related effects, took only ∼10 min on a single core. This is not surprising, since the numerical complexity of the ray-tracing is O(), i.e. it is linear with respect to the number of rays, as in the case considered in the previous section, whereas the complexity of the partially coherent SRW simulation is , where is the number of non-correlated wavefronts used, and N is the number of data points in each wavefront. In the low-coherence case considered in this section, had to be very large for ensuring numerical convergence, so that . Note that in cases of higher coherence the required can be considerably smaller than this case.

Even though the partially coherent wavefront propagation simulations are currently feasible even for low-coherence radiation beams, their relatively low efficiency is one of the main reasons why this more accurate calculation method is currently less often used than the less accurate geometrical ray-tracing method. Increasing the efficiency of the partially coherent simulations is one of important directions of algorithmic developments in numerical applications of physical optics.

3.3. Partially coherent undulator radiation in a high-resolution microspectroscopy beamline

Until now we were considering the propagation of fully and partially coherent Gaussian beams through simple optical systems containing KB mirrors as the only focusing elements. In this section we increase the level of detail and complexity of both the source and optics, which are now essentially matching those of the Sub-Micron Resolution X-ray Spectroscopy (SRX) beamline at NSLS-II (De Andrade et al., 2011). To facilitate focusing of X-rays down to 10–100 nm spot sizes with a possibility of controlling the resolution and the flux at sample, the SRX beamline uses a popular optical scheme. A horizontal focusing mirror (HFM) creates an intermediate waist in the horizontal plane, where a variable-size secondary source aperture (SSA) is installed. The SSA can be fully open (for high-throughput experiments) or closed to any size down to ∼5 µm (for the experiments requiring high spatial resolution). Different KB mirror systems are considered to be used in the two types of experiments. In all cases, the VKB creates an image of the primary source whereas the HKB re-focuses the SSA at the same longitudinal position where the VKB focus is located, i.e. the sample position. This is illustrated in Fig. 6. The positions of the optical elements relative to the source are as follows: the HFM is located at 32.95 m, the SSA at ∼49.35 m, the VKB at ∼64.808 m, the HKB at ∼64.915 m and the sample position at ∼64.995 m. The VKB is 149 mm long and the HKB is 59 mm long. All mirrors are installed at a grazing angle of 3.5 mrad. For simplicity, the monochromator, an important element of that microspectroscopy beamline, was assumed to be ideal and was not included in the simulations.

Figure 6 Simplified optical scheme of the SRX beamline of NSLS-II (the monochromator is not shown).

The radiation source used by this beamline is a 21 mm-period 1.5 m-long in-vacuum undulator U21 installed in a low-β straight section of NSLS-II. In the center of this section, when the long-term target horizontal emittance of 0.55 nm will be reached, the horizontal (vertical) RMS size of the electron beam will be ∼33.3 µm (∼2.9 µm) and the horizontal (vertical) RMS divergence ∼16.5 µrad (∼2.7 µrad). The simulations were performed for 7.18 keV photon energy, attainable at the fifth harmonic of the U21 radiation spectrum. The undulator radiation for all types of simulations described in this section, wavefront propagation and geometrical ray-tracing, was calculated using the SRW code (Canestrari et al., 2013).

The results of the simulations with the three methods for four different SSA sizes are presented in Fig. 7 and in Table 1. We note that in all cases the radiation beam is limited along both transverse directions by an aperture located in front of the VKB. The vertical and horizontal sizes of the aperture were chosen to approximately match the VKB and HKB dimensions projected on the transverse plane: ∼0.52 mm (vertical) and ∼0.21 mm (horizontal); whereas the corresponding vertical and horizontal FWHM sizes of the beam before the VKB are ∼1.8 mm and ∼1.6 mm, respectively. The beam truncation by the aperture results in diffraction effects appearing in the wavefront propagation simulations but not in the ray-tracing. As one can see from Fig. 7(a), this difference is quite noticeable even for a relatively large 50 µm SSA size, which is larger than the ∼39 µm FWHM size of the horizontal waist at the SSA. We note, however, that in this partially coherent radiation beam case the effect of diffraction on the KB aperture is less significant than in the case of the fully coherent beam (see Fig. 3 in §3.1).

Figure 7 Horizontal and vertical cuts of intensity distributions of the partially coherent X-rays at 7.18 keV photon energy at the position of geometrical focus of KB mirrors of the SRX beamline for four SSA sizes: (a) 50 µm, (b) 20 µm, (c) 10 µm and (d) 5 µm.

With the reduction of the SSA size from 50 µm to 5 µm, the horizontal size of the focused spot resulting from the wavefront propagation simulations gradually approaches the diffraction limit, whereas the spot size predicted by the geometrical ray-tracing keeps on reducing down to non-physical small values way beyond the diffraction limit. Besides, the peak of intensity in the focal spot predicted by the ray-tracing is practically not dependent on the SSA size, which is not the case in the results of the wavefront propagation simulations, predicting a (physically correct) reduction of the peak intensity because of the diffraction at the SSA. As one can see from Fig. 7, the difference between the wavefront propagation simulation results for the elliptical mirrors and the ideal thin lenses is not very large in this case. The use of the wavefront propagation method, even with the KB mirrors simulated by ideal thin lenses, is more appropriate than the use of the pure geometrical ray-tracing, which does not provide physically consistent results in this case (unless a special post-processing/correction is made).

The numerical convergence of the partially coherent calculations with the SRW code described in this section was somewhat faster than in the low-coherence case described in the previous section, especially for the focal-spot intensity distributions approaching the diffraction limit, which happened when the horizontal SSA size was 5 µm. It is noteworthy that the latter is exactly the case posing certain difficulties for the ray-tracing simulations, because the portion of rays passing through the tiny aperture is small (∼14% of the total number of rays). The wavefront propagation simulation results illustrated in Fig. 7 were obtained by averaging ∼10⁵ wavefronts from individual `macro-electrons', distributed according to a Gaussian over the five-dimensional phase space of the electron beam. The phase space `coordinates' were the horizontal and vertical positions and angles, and the energy of the electrons. One partially coherent simulation typically lasted several hours on 64 nodes of a cluster. In each geometrical ray-tracing simulation, performed using the SHADOW code, ∼10⁸ rays were used. Each simulation typically lasted ∼20 min on a single core.

4. Summary

A new method for simulation of fully and partially coherent wavefront propagation through extended grazing-incidence mirrors has been recently implemented in the SRW computer code. The method makes use of a local stationary phase expansion to propagate an electric field from a plane before such a mirror to a plane after it. By using a two-dimensional interpolation, the electric field can be re-sampled after such local propagation on a regular mesh. The operation can be combined with the `free-space' and `thin optical element' propagators and can significantly extend the applicability of numerical Fourier optics based simulations to complex optical systems without compromising the overall performance. Thanks to such seamless integration, the new method does not complicate the partially coherent radiation propagation calculations that are possible with the SRW code.

The new method has been extensively benchmarked against two other methods which are widely used for simulating propagation of optical beams: wavefront propagation using ideal thin lenses instead of grazing-incidence mirrors (with the SRW code) and the geometrical ray-tracing taking into account accurate mirror surface shapes (with the SHADOW code). The benchmarking simulations were performed for realistic X-ray optical systems including elliptical mirrors in the popular Kirkpatrick–Baez geometry at high demagnification, with fully and partially coherent Gaussian beams and an undulator radiation source. Good agreement between the results of the simulations using the new method and the geometrical ray-tracing is observed in the cases without diffraction at mirror apertures, even for fully coherent radiation beams (given the appropriate choice of the input beam sizes and divergences for the ray-tracing simulations). In the cases dominated by diffraction, the differences between the predictions of the wavefront propagation simulations using the new method and the geometrical ray-tracing were quite significant. The results obtained with the new method in these cases were closer to those obtained using the wavefront propagation with the mirrors approximated by ideal thin lenses. All results obtained have a clear physical interpretation, which suggests that the new method can be used for miscellaneous optical simulations, and in particular for those targeting design and optimization of X-ray beamlines and instruments for modern third- and fourth-generation sources.

The benchmarking of the new method being discussed, as well other methods that were recently implemented (or are currently under implementation) in the SRW code, is also planned to be compared with the experimental results obtained at new low-emittance synchrotron sources and X-ray FELs.