1. Introduction

High accuracy calculations of the partially coherent Synchrotron Radiation (SR) propagation through optical beamlines and experimental end-stations of light sources play an important role in a large variety of applications of these scientific facilities. Optimal design of X-ray beamlines for different types of experiments, simulation of experiments to estimate their feasibility on existing beamlines, finding the most appropriate beamline settings for given experiments without spending expensive beam time, and experimental data processing are the examples of such applications. With the continuing reduction of electron beam emittance in modern synchrotron light sources, targeting to increase brightness and coherence of emitted X-rays, the use of the wave optics methods for such calculations becomes increasingly important. Many of the popular experimental techniques, such as coherent X-ray scattering, high-resolution X-ray microscopy, coherent diffraction imaging, ptychography entirely rely on wave-optics features of X-rays. Therefore, for the efficient development and improvement of these techniques, the use of wave-optics calculations, with accurate tracking of radiation degree of coherence, is mandatory. Such partially coherent calculations are, on the other hand, known to be numerically intensive, and the improvement of their efficiency is an important topic of developments in this area.

The processes of single-electron SR emission and propagation can be accurately described in the framework of classical electrodynamics [1], with the retarded potentials based approach being used for the emission [2], and numerically-efficient Fourier optics [3] (extended with compatible methods) for the propagation calculations. Accurate treatment of the finite electron beam emittance, having major impacts on the SR brightness and coherence in storage ring SR sources, can be done by the averaging of SR characteristics, e.g. intensity or cross-spectral density (CSD), also referenced as mutual optical intensity (eventually after a propagation through an optical system) over the phase space occupied by the emitting electron beam [4–6]. This “direct” method, however, has a relatively low numerical efficiency: even if the calculation of SR electric field from an individual electron and its subsequent propagation through beamline optics can be done quickly (within seconds at a sequential execution), accurate treatment of electron beam emittance effects requires performing such single-electron calculations at different initial conditions of electrons a large number of times (up to hundreds of thousands or millions for hard X-rays in 3rd generation SR sources). Though fully feasible for the 3rd and more recent generations of SR sources [7–16] and easily parallelizable, such calculation may require quite extensive computational resources – many hours of computation even at a parallel execution on a multi-core computer server or cluster.

A general method for a compact representation of spatially partially coherent radiation based on coherent modes, that can be determined based on a coherent mode decomposition (CMD) procedure was introduced in [17,18]. This representation was applied to partially coherent undulator radiation (UR) in a 3rd generation storage ring [18], with an experimental method being used for determining the coherent modes (without the calculation of the UR CSD). The CMD based calculation method was proposed to be used for SR sources in [20], as it promised a considerable reduction in the number of fully coherent wavefront (electric field) propagations, compared to the direct method [7], provided that the SR coherent modes are known. An implementation of the CMD of partially coherent UR, and the coherent mode propagation tests were described in [21]. For obtaining 2D coherent modes of UR, authors calculated CSD on a dense 4D mesh, required for numerically “resolving” the oscillatory behavior of the CSD as a function of horizontal and vertical mutually conjugate coordinates. As a result, the 4D CSD of UR in a relatively low-emittance 3rd generation SR source occupied terabyte-size memory, and for the CMD, a parallel software method, executing on a large computer cluster, had to be used. The overall numerical efficiency of partially coherent UR calculations including the CMD in this implementation was comparable to that of the direct calculation [7].

In this paper, we describe a CMD method, which has much lower memory requirements for the initial 4D CSD (from several hundred MB to ∼10 GB, instead of TB) and can be executed on a computer with one multi-core CPU. The reduction of the memory requirements becomes possible thanks to analytical treatment of the common quadratic phase terms, that develop in the CSD (the same way as in electric fields of single-electron emission) at some distance from the radiation source. Along with the direct method [7] and the method of performing the CMD at a radiation beam waist [22], the new CDM based method described in this paper is implemented in the “Synchrotron Radiation Workshop” (SRW) software package [23]. The calculation examples described in this paper refer to the Hard X-ray Nanoprobe (HXN) beamline operating at the National Synchrotron Light Source II (NSLS-II) at the Brookhaven National Laboratory. Even though this work is presented mostly in the context of applications of synchrotron light sources, we believe that the calculation method discussed is general and can be used for different sources of partially coherent radiation.

2. Analytical treatment of quadratic phase terms in the electric field of fully coherent radiation

One of the known difficulties of performing numerical wave-optics calculations is related to the necessity of dealing with (strongly) oscillating functions describing different polarization components of the radiation complex electric field. To “resolve” numerically the fast oscillations of these electric fields, e.g. in a transverse plane perpendicular to the optical axis, one may need to use dense grids with large numbers of data points in the two transverse directions, largely increasing the computer memory and CPU resource consumptions, even if the relatively efficient numerical apparatus of Fourier optics is used.

Since the leading phase terms of electric field components of fully coherent radiation in the largely applicable paraxial approximation are usually known a priori (even if approximately), their analytical treatment can be applied to the numerical calculation of the wavefront propagation [24,25]. In this approach, a component of a complex frequency-domain electric field $U(x,y)$ is represented in a transverse plane as:

As shown in the next section, a similar approach can be used for reducing the problem size of performing CMD of 4D CSD of partially coherent radiation at a distance from source (or from radiation beam waist).

3. Coherent mode decomposition of partially coherent radiation

3.1 Conventional method

The CSD of one polarization component of temporally incoherent SR at a fixed frequency / photon energy, emitted by a relativistic electron beam moving in magnetic field of an insertion device or a bending magnet of a storage ring, considered as a function of horizontal and vertical transverse coordinates $W({x_1},{y_1},{x_2},{y_2})$ in a given observation plane, can be calculated by integrating contributions from the CSDs corresponding to individual electrons distributed over the phase space volume of the emitting electron beam (see e.g. [4], though the “CSD” term / abbreviation was not used there):

We note that Eq. (2) can be readily used for calculating different polarization components of the CSD (with the corresponding electric field components in place of U), and it can be easily modified for the calculation of generalized Stokes parameters [26]. The relations similar to Eq. (2) can be used to calculate other characteristics of partially coherent SR, such as the intensity distribution, i.e. the “diagonal” term of the CSD when ${x_1} = {x_2},\textrm{ }{y_1} = {y_2}$, and Wigner distribution (i.e. 2D Fourier transform of the CSD, see [4]). The CSD allows to estimate the degree of coherence and coherence lengths in the two transverse directions – the very important characteristics of the partially coherent radiation (see section 4).

The electric field U and the CSD distribution W in Eq. (2) can be considered at any location of optical scheme, e.g. before first optical element of a beamline, or after propagation through the entire beamline, at an experimental sample. The numerical method for calculating the intensity of partially coherent SR after propagation through a beamline has been implemented in SRW code based on Eq. (2) for more than a decade [6,7]. Its numerical efficiency, however, is not very high, especially for relatively low-coherence SR sources in the hard X-ray spectral range (when hundreds of thousands or millions of propagations of electric fields emitted by electrons with different initial conditions distributed over the phase space of the electron beam may be required for achieving a desired accuracy of the propagated radiation intensity).

It should be mentioned, that in some cases, like propagation of radiation from a source in free space or to a waist after the focusing by a lens, or after diffraction on one aperture, convolution-type relations with respect to the single-electron intensity may hold, allowing for a dramatic acceleration of numerical evaluation of Eq. (2). However, this calculation method is not general (as the convolution relation doesn’t hold for large numbers of types of sources and optical systems).

Several approaches to improving the numerical efficiency of partially coherent SR calculations were suggested, based on the use of the decomposition of CSD over pre-defined orthogonal sets of functions, such as Gauss-Hermite and / or Gauss-Laguerre modes [20]. Since these functions are not necessarily the eigenmodes of the SR CSD, a large number of modes / coherent wavefront propagations may be required for reaching a sufficient accuracy in a general case of 4D CSD of SR. The general CMD method [18] was applied to the CSD of the UR calculated at a distance from the undulator in [21]. Following this general method, the CMD is formulated as:

Since the electric fields emitted by different electrons $U(x,y,\Omega )$ in the monochromatic case, considered in this paper, are functions of the horizontal and vertical coordinates, they need be sufficiently sampled to resolve all their oscillations along these dimensions. Their CSD distribution, which is a function of four coordinates, may require a huge amount of the computer memory: in cases where the oscillation rates of the electric field U are dominated by their quadratic phase terms (see Eq. (1)), the number of CSD data points describing the partially-coherent radiation may need to be as large as ${\sim} {({{\Delta \theta _x^2\Delta \theta _y^2{R_x}{R_y}} / {{\lambda ^2}}})^2}$, where $\Delta {\theta _x},\Delta {\theta _y}$ are angular dimensions of the radiation wavefronts, ${R_x},{R_y}$ their corresponding radii of curvature, λ the radiation wavelength (see Eq. (5) in [25]). This means ∼6 × 10¹² single-precision complex numbers, or ∼5 TB of data in the case of a ∼full fan of useful partially coherent UR at 9 keV photon energy considered at a position of first optical element of a hard X-ray beamline at a storage ring light source like NSLS-II. Such a memory requirement considerably complicates the numerical procedure of the eigenvalue problem, requiring special software packages and a lot of CPU time even when using parallelization, which compromises the benefits from the CMD strategy.

3.2 Method with the subtraction of common quadratic phase terms

Equation (1) can be formally applied to the electric fields emitted by individual electrons contributing to the partially coherent radiation beam, where ${R_x}$ and ${R_y}$ are approximate wavefront radii of curvature in the horizontal and vertical planes, which are common for all electrons, and $({x_0},{y_0})$ are the coordinates of the source point averaged over the entire beam (that are also common for all electrons). With this representation, Eq. (2) for the CSD can be re-written as:

By substituting Eq. (6) to Eq. (4), we readily obtain:

Substituting Eq. (5) and (8) to Eq. (3) further gives:

I.e., ${\tilde{\psi }_n}(x,y)$ are eigenfunctions of Eq. (7), corresponding to the same eigenvalues ${\lambda _n}$ as the eigenfunctions ${\psi _n}(x,y)$ of Eq. (4). Solving the eigenvalue problem of Eq. (7) is therefore equivalent to solving the eigenvalue problem of the original Eq. (4). The numerical sampling density of the function $\tilde{W}({x_1},{y_1},{x_2},{y_2})$, due to the absence of the leading phase terms resulting in oscillations, on the other hand, can be much lower than that of the original CSD,, where these terms are present. From the point of view of numerical calculations, it can therefore be advantageous to solve Eq. (7) instead of Eq. (4). In the next section, we illustrate this by simulations performed for partially coherent X-ray beams of UR produced at the NSLS-II storage ring.

4. Calculation example of coherent mode decomposition of undulator radiation at NSLS-II

In this section, we discuss parameters of the NSLS-II electron beam for which all numerical simulations were performed, illustrate the CSD of partially coherent UR calculated for these parameters with and without the subtraction of common quadratic phase terms, and present results of CMD performed on the CSD calculated with the subtraction of common quadratic terms on meshes with the same ranges of transverse position, but different (reduced) numbers of points.

4.1 Parameters of NSLS-II electron beam

The main parameters of electron beam in a Low Beta straight section, that is used for brightness and coherence demanding undulator based X-ray beamlines at NSLS-II, are listed in Table 1. As in all 3rd generation synchrotron light sources, the emittance, and separately size and divergence of the electron beam in NSLS-II in the horizontal direction are much larger than the corresponding parameters in the vertical direction. Because of this, the radiation degree of coherence and transverse coherence length in the hard X-ray spectral range are also considerably higher in the vertical direction than in the horizontal one.

Table 1. NSLS-II electron beam parameters: general and related to low-beta straight section.

View Table

To produce high brightness radiation in the hard X-ray range, In-Vacuum Undulators (IVUs) are used in modern medium energy storage rings such as NSLS-II. For our simulations, we used parameters of 20 mm period 3 m long IVU20 undulator, the radiation source of the Hard X-ray Nanoprobe (HXN) high-resolution microscopy beamline at NSLS-II. We will first discuss the UR CSD calculations for the HXN beamline, and in section 5 we will discuss in detail the partially coherent UR propagation simulations performed for this beamline, without and with the use of CMD. We did our simulations for 9 keV photon energy, that is frequently used at HXN. This photon energy can be reached at 7th harmonic of UR from the IVU20.

4.2 Partially coherent undulator radiation cross-spectral density calculations with and without the subtraction of common quadratic phase terms

In this work, we used the CSD calculation method based on Eq. (2). The calculations were done using a Message Passing Interface (MPI) parallel computing method, implemented in the Python interface in SRW. To efficiently improve the code performance, the computation was split among several MPI “groups”, each consisting of several “worker” processes performing UR (and, if necessary, fully coherent wavefront propagation) calculations on a subset of electrons, and one “master” process calculating the corresponding “partial” 4D CSD according to Eq. (2). The CSD data generated by the groups were averaged, producing the final 4D CSD.

The calculation based on this algorithm was taking up to several tens of hours on one multi-core computer server for ∼3 × 10⁵ of “macro-electrons” contributing to the CSD (with variations depending on the mesh chosen for the final 4D CSD). In an alternative implementation, the calculation of 4D CSD from electric fields according to Eq. (2) was done on GPU by one “master” process; in that case, it was possible to run the entire calculation with only one MPI “group” without creating performance “bottleneck”. These general methods are applicable for calculating CSD at any location of a beamline (at a source, at some distance from it, after a number of optical elements, at a sample, etc.). In some special cases, e.g. at the observation plane located at a distance from undulator, without any optical elements in optical path, more numerically-efficient methods can be developed, making use of convolution-type relations [21]. A similar method is also under development in SRW.

Examples of cross-sections / cuts of the UR 4D CSD calculated for the electron beam and undulator parameters described above are presented in Fig. 1. The calculations were performed using Eqs. (2) and (6) for the 20.5 m observation distance from the undulator. The calculations were done on the transverse meshes with 512 points vs. each of the horizontal mutually conjugate coordinates $({x_1},{x_2})$ and 60 points vs. each of the vertical mutually conjugate coordinates $({y_1},{y_2})$. This is much less than the number of points required for fully “resolving” the single-electron SR electric fields and CSD without subtraction of the common quadratic phase terms: according to an estimate based on Eq. (5) of [25], 2220 (h) x 356 (v) points would be required in the latter case. Image plots in Figs. 1(a) and 1(c) show 2D cuts of the 4D CSD vs. the horizontal coordinates $({x_1},{x_2})$ at the vertical coordinates close to 0 (Fig. 1(a)) and vs. the vertical coordinates $({y_1},{y_2})$ at the horizontal coordinates close to 0 (Fig. 1(c)). The image plots on the left in Figs. 1(a) and 1(c) depict 2D cuts of the CSD after the subtraction of the common quadratic phase terms, whereas the image plots in the center show the cuts of the CSD made on the same mesh without the subtraction of these terms, and the plots of the right show the cuts of the CSD still without the subtraction of the common quadratic terms, but on a much denser mesh. The comparison of the plots on the right with those in the center clearly demonstrate that the CSD calculated on the mesh with (512 × 60)² points without the subtraction of the common quadratic terms (plots in the center) is not sufficiently well “resolved” (as the image patterns are strongly different and not smooth). On the other hand, the cuts of the CSD after the subtraction of the common quadratic terms (plots on the left) are very smooth and well “resolved”. The 1D cuts of the 4D CSD vs. the horizontal position x₁ (Fig. 1(b)) and vs. vertical position y₁ (Fig. 1(d)) confirm this observation: the number of points (512 × 60)² is clearly insufficient for resolving all oscillations of the CSD without the subtraction of the common quadratic phase terms, as this can be particularly well seen in the cuts made at x₂ and y₂ values equal to the halves of the corresponding maximum positions (x₂ ≈ 0.5 x_max and y₂ ≈ 0.5 y_max, see graphs on the right in Figs. 1(b) and 1(d)).

 

Fig. 1. SD distribution of partially-coherent UR at 9 keV, the resonant photon energy of 7th harmonic of 20 mm period 3 m long in-vacuum undulator IVU20, calculated for 20.5 m observation distance from the undulator center, which is also the center of the Low Beta straight section of NSLS-II: (a) vs. mutually conjugate horizontal coordinates x₁, x₂ at fixed conjugate pair of vertical coordinates y₁ and y₂ close to 0, with the average quadratic phase terms subtracted (left), with all phase terms present over the mesh with 512 points vs. x₁ and x₂ used for the initial 4D CSD computation (center), and over a more dense mesh with the CSD values obtained by interpolation (right); (b) vs. x₁ at x₂ value of the mesh close to 0 (left), and at x₂ ≈ 0.5 x_max (right), both at y₁, y₂ close to 0; (c) vs. mutually conjugate vertical coordinates y₁, y₂ at x₁, x₂ values of the mesh close to 0, with the average quadratic phase terms subtracted (left), with all phase terms present over the mesh with 60 points vs. y₁ and y₂ used for the initial 4D CSD computation (center), and over a more dense mesh with the CSD values obtained by interpolation (right); (d) vs. y₁ at y₂ value of the mesh close to 0 (left), and at y₂ ≈ 0.5 y_max (right), both at x₁, x₂ close to 0. The legend in the left graph (b) applies to all graphs in (b) and (d).

Download Full Size | PDF

Figures 2 and 3 show important UR characteristics obtained directly from the 4D CSD illustrated in Fig. 1: the UR intensity distribution vs. horizontal and vertical coordinates $(x,y),$ i.e. the diagonal term of the 4D CSD at x₁ = x₂ = x and y₁ = y₂ = y, and the UR degree of coherence in the horizontal and vertical mid-planes, calculated from the CSD according to:

 

Fig. 2. UR intensity (diagonal term of the 4D CSD illustrated in Fig. 1 at x₁ = x₂ = x and y₁ = y₂ = y): the 2D distribution (image plot on the left) and its cuts by the horizontal (graph in the center) and vertical (graph on the right) mid-planes.

Download Full Size | PDF

 

Fig. 3. UR degree of coherence obtained from the 4D CSD illustrated in Fig. 1: (a) vs. (x₁+ x₂)/2 and (x₁-x₂)/2 (left), and vs. (x₁-x₂)/2 at x₁ = -x₂ (right), both at y₁ = y₂ ≈ 0; (b) vs. (y₁+ y₂)/2 and (y₁-y₂)/2 (left), and vs. (y₁-y₂)/2 at y₁ = -y₂ (right), both at x₁ = x₂ ≈ 0. Box-shaped markers correspond to the points of the mesh used for the CSD and CMD calculation after the subtraction of the common quadratic phase terms (512 mesh-points vs. horizontal and 60 vs. vertical coordinates).

Download Full Size | PDF

For convenience, the degree of coherence is plotted as function of combinations of the mutually conjugate variables: ${{({x_1} + {x_2})} / 2},$ ${{({x_1} - {x_2})} / 2},$ ${{({y_1} + {y_2})} / 2},$ ${{({y_1} - {y_2})} / 2}.$ The FWHM sizes of the central “lobes” of the degree of coherence cuts vs. ${{({x_1} - {x_2})} / 2}$ and ${{({y_1} - {y_2})} / 2}$ can be interpreted as the horizontal and vertical transverse coherence lengths. In the case under consideration, these lengths are equal to ∼14 µm and ∼146 µm. Note that these values are smaller than the FWHM sizes of the UR intensity distribution shown in Fig. 2, that are ∼1.04 mm and ∼448 µm, respectively. The difference between the UR beam size and its transverse coherence length is particularly large in the horizontal plane. This is related to the fact that the electron beam emittance and size in the horizontal direction are much larger than in the vertical one (see Table 1), which results in the considerably lower radiation coherence in the horizontal plane. This is typical for 3rd generation SR sources, and should be taken into account in optical designs of X-ray beamlines at these sources. Figures 2 and 3 clearly show that the (relatively sparse) mesh with (512 × 60)² points that was used for the 4D CSD calculation is perfectly sufficient for analyzing both the UR intensity distribution and its degree of coherence in the case under consideration.

4.3 Coherent mode decomposition of the total undulator radiation beam

The CMD was applied to the 4D CSD calculated on the sparse mesh consisting of (512 × 60)² points, with the subtraction of the common quadratic terms, that was illustrated in section 4.2. We note that in our calculations, the ranges of the horizontal and vertical positions were 2.5 mm and 1.0 mm respectively, i.e. large enough to cover nearly the total “footprint” made by partially coherent UR at 9 keV photon energy at 20.5 m from the undulator in the storage ring with electron beam parameters listed in Table 1, see Fig. 2. Reduction of the transverse position ranges, if compatible with X-ray beamline acceptances, would reduce the requirement to the numbers of mesh points on which the CSD has to be calculated, and would enable the CMD to be performed effectively without the subtraction of quadratic phase terms (see e.g. [19]).

In our work, the CMD was done by solving Eq. (7) in its matrix form for the eigenmodes. Three numerical algorithms were successfully applied for this: two from the SciPy package, that use the ARPACK Fortran package [27]: “scipy.linalg.eigh” and “scipy.sparse.linalg.eigsh” methods of the Linear Algebra and Sparse Linear Algebra modules, and the “primme.eigsh” method of the PRIMME library [28], via a wrapping module implemented in Python in SRW code. The time to solution was typically under half an hour on one multi-core server. In the above case of a medium-sized mesh, the fastest results were delivered by the “primme.eigsh” method. In general, differences among the computational efficiencies of these methods are present, depending on the problem size and the CSD matrix sparsity.

The intensity distributions of a subset of obtained coherent modes, sorted (with their corresponding eigenvalues) according to the flux “carried” by them, using 1-based numbering, are presented in Fig. 4. With such sorting, modes contributing to the radiation beam intensity in the horizontal plane appear to have first numbers (see the upper row of image plots in Fig. 4), whereas the modes contributing to the beam intensity in the vertical plane have larger numbers (see lower row of plots in Fig. 4). In this calculation, 1000 coherent modes were considered and used in subsequent analysis.

 

Fig. 4. Examples of intensity distributions of coherent modes of UR obtained by decomposition of 4D CSD illustrated in Fig. 1, with the subtraction of the common quadratic phase terms on the mesh with 512 (h) x 60 (v) points.

Download Full Size | PDF

The flux, carried by these modes individually, and the integrated (summed-up) flux, are shown in Fig. 5 (in graphs on the left and on the right respectively). The graph on the left exhibits an exponential decay of the high order modes contribution to the total flux. As one can see from the graph on the right in Fig. 5, the total flux carried by the first 1000 coherent modes of the UR at 9 keV at the NSLS-II storage ring, is still by ∼5% smaller than the total flux according to the direct calculation of the UR intensity distribution or 4D CSD (the latter value is shown by horizontal line in the right graph).

 

Fig. 5. Spectral flux of coherent modes of UR obtained by decomposition of 4D CSD illustrated in Fig. 1 with the subtraction of the common quadratic phase terms on the mesh with 512 (h) x 60 (v) points (left), and the corresponding integrated spectral flux (right) in comparison with that of the total UR beam, according to the direct calculation of CSD and intensity distributions illustrated in Figs. 1 and 2.

Download Full Size | PDF

Figure 6 shows the integrated (summed-up) 2D intensity distributions (Fig. 6(a)), their 1D cuts vs. the horizontal and vertical positions (Fig. 6(b)), and the degree of coherence (Fig. 6(c)) compared among subsets of major modes as well as against the direct calculations. As the number of major modes increases, both the total intensity and degree of coherence recover those from the direct calculation, with the sizes of intensity distributions gradually increasing (see Fig. 6(b)), and the sizes of the central “lobes” of the degree of coherence curves (i.e. the coherence lengths, see Fig. 6(c)) decreasing with the number of contributing modes. As one can see from the graph on the right in Fig. 6(b), the intensity corresponding to the sum of 1000 modes differs from the direct calculation essentially only at large vertical coordinate values (|y| > 350 µm), which are larger than the vertical acceptance of the HXN beamline, for which we did the analysis in this paper. This suggested that limiting the number of coherent modes in this analysis to 1000 probably resulted in smaller than 5% loss of accuracy of wavefront propagation calculations.

 

Fig. 6. UR characteristics determined as sums of contributions from different numbers of coherent modes, in comparison with direct calculations (without the CMD): (a) image plots of 2D intensity distributions vs. horizontal and vertical positions; (b) cuts of the intensity distributions by the horizontal (left) and vertical (right) mid-planes; (c) degree of coherence in the horizontal (left) and vertical (right) mid-planes.

Download Full Size | PDF

4.4 CMD performed on CSD calculated on sparse meshes

The method of performing CMD on 4D CSD with the subtraction of common quadratic phase terms was tested with the CSD calculated on transverse meshes with different, gradually reduced, numbers of points, but with the same ranges vs. the horizontal and vertical positions, 2.5 mm (h) x 1.0 mm (v). The numbers of points vs. the horizontal positions in these calculations varied from 512 down to 128, and vs. the vertical positions from 60 down to 40. The 4D CSD calculated on the mesh with (512 × 60)² points occupied ∼7.5 GB of memory, whereas with the (128 × 40)² point mesh, it occupied only ∼210 MB. We note that all these numbers are dramatically smaller than the TB level memory requirement for the 4D CSD that would need to be applied if the CMD would be done for the UR parameters considered in this paper without the subtraction of the common quadratic phase terms. The small numbers of mesh points also lead to fast CMD time to solutions (within a magnitude of seconds in the case of the (128 × 40)² point mesh).

Horizontal and vertical mid-plane cuts of intensity distributions obtained by summing-up contributions of different numbers of coherent modes, from the 4D CSD calculated on meshes with different numbers of points (with the subtraction of the common quadratic phase terms), are presented in Fig. 7.

 

Fig. 7. Horizontal (left) and vertical (right) mid-plane cuts of UR intensity distributions determined from different numbers of coherent modes obtained from 4D CSD calculated on meshes with different numbers of points, from 512 (h) x 60 (v) down to 128 (h) x 40 (v), for the same ranges of the horizontal and vertical positions.

Download Full Size | PDF

Overall good agreement between the partial intensity distributions created by modes calculated on meshes with different numbers of points can be observed, with an exception for the 100 mode case, when the sum of the modes, calculated on the most sparse mesh with (128 × 40)² points, results in ∼7% higher peak intensity (see blue dashed curves in the graphs) than sums of the 100 modes obtained from CSD on meshes with larger numbers of points. These intermediate results suggest that the coherent modes, obtained with the subtraction of the common quadratic phase terms in all the considered mesh cases, can in principle be used for subsequent wavefront propagation calculations. The results of such tests are described in section 5.4 below.

5. Examples of CMD calculation after X-ray beam propagation through a beamline

In this section, we describe partially coherent UR propagation calculations for the HXN beamline at NSLS-II. After a short description of the HXN optical layout, we present the results of 4D CSD calculation at the beamline sample position (at the end of its optical path), then show and discuss the results of the CMD performed at that location. First, we perform these calculations using the direct method [7], i.e. by calculating single-electron UR at different initial conditions of electrons, propagating these electric fields through the beamline to its sample position, calculating the corresponding 4D CSD, and applying the CMD to it. For comparison, we also calculate the same 4D CSD and CMD at the sample position by propagating the coherent modes, that were derived at the beamline entrance and discussed in the previous section.

5.1 Optical layout of hard X-ray nanoprobe beamline used for the simulations

The pictogram of a simplified optical layout of the HXN beamline, generated by the Sirepo web-browser GUI [29] supporting wave-optics simulations with SRW, is shown in Fig. 8. The S1 slits located at ∼20.5 m distance from the center of undulator have 2.5 mm (h) x 0.7 mm (v) aperture for the X-ray beam. The horizontal collimating mirror (HCM) with a focal length of ∼28.35 m (equal to the distance between the mirror and undulator centers) creates nearly parallel X-ray beam, that passes through a double-crystal monochromator (not shown in Fig. (8)) and is focused (in the horizontal plane) by the horizontal focusing mirror (HFM) on the secondary source aperture (SSA). In the vertical plane, the X-ray beam is also focused on the SSA by a compound refractive lens (CRL). The SSA dimensions in the operation mode, for which the calculations were performed, were equal to 20 µm (h) x 100 µm (v). The X-ray beam focused on the SSA is truncated (especially in the vertical plane) by the SSA and propagated (with a divergence) further to the location of the aperture of the final focusing optics (AFFO) and the FFO itself, which is simulated by an ideal lens. The dimensions of the AFFO in our simulations were 150 × 150 µm². The FFO (which can be a multilayer Laue lens [30] or a zone plate) produces nano-scale size X-ray beam spot at ∼18 mm from it, at the “sample” position.

 

Fig. 8. Pictogram of a simplified optical scheme of HXN beamline with the optical element positions with respect to undulator center shown in the row on the top (see explanations in text).

Download Full Size | PDF

This type of optical scheme is used relatively frequently for high-resolution microscopy beamlines in 3rd generation SR sources. Among its advantages is a possibility to easily switch between the “high resolution” mode (and ultimately even reach the diffraction limited resolution, if the quality of FFO allows for this) and the “high throughput” mode, or / and “trade” between the spatial resolution and flux at sample, as needed by a given experiment, by simply changing the SSA dimensions [8].

5.2 X-ray beam cross-spectral density after propagation to beamline sample position

To calculate 4D CSD at the sample position of the HXN beamline, we applied the wavefront propagation method, first using the direct calculation of UR electric fields from macro-electrons, propagating them through the beamline and calculating their contributions to the CSD [7]. This calculation is similar to the one described in [8], with the difference that instead of intensity, 4D CSD is calculated after the propagation of individual electric fields, as described in section 4.2. The CSD was calculated on the final 4D mesh with (100 × 100)² points (with the data occupying 800 MB of memory) and the same range vs. horizontal and vertical position: −100 nm ≤ x_1,2 ≤ 100 nm, −100 nm ≤ y_1,2 ≤ 100 nm.

Figure 9 shows 2D and 1D cuts of the 4D CSD distribution calculated at 9 keV photon energy at the sample position of the HXN beamline. The cuts vs. horizontal and vertical positions, made close to the center and off the center of the CSD distribution are well resolved on the mesh with relatively small numbers of points that were chosen for this calculation. The same way as the electric field of a fully coherent radiation beam at its waist, the CSD of the partially coherent beam does not have strong quadratic terms in its phase at the waist and doesn’t require large numbers of points for being well resolved (this feature was effectively exploited in [22]).

 

Fig. 9. Cuts of 4D CSD distribution of partially coherent 9 keV X-ray beam at the sample position of HXN beamline: (a) vs. mutually conjugate horizontal coordinates x₁, x₂ (left), and at x₂ ≈ 0 and x₂ ≈ 0.5 x_max (right), at y₁, y₂ values of the mesh close to 0; (b) vs. mutually conjugate vertical coordinates y₁, y₂ (left), and at y₂ ≈ 0 and y₂ ≈ 0.5 y_max (right), at x₁, x₂ values close to 0.

Download Full Size | PDF

The 2D intensity distribution, i.e. the diagonal term of the 4D CSD at sample position, is presented in Fig. 10 on the left, and the horizontal and vertical 1D mid-plane cuts of this distribution are shown in graphs in the center and on the right. The choice of the SSA dimensions and other setting of the HXN beamline allowed to obtain X-ray spot with nearly equal horizontal and vertical FWHM sizes, ∼23 nm (h) x 24 nm (v) in the case under consideration. The degree of coherence plots, presented in Fig. 11, show that the transverse coherence lengths of this radiation (∼16 nm (h) x 17 nm (v)) are comparable to the intensity spot dimensions. This corresponds to the radiation with relatively high degree of coherence.

 

Fig. 10. X-ray intensity at the sample position of HXN beamline obtained directly from the 4D CSD illustrated in Fig. 9 (by considering its diagonal term at x₁ = x₂ = x and y₁ = y₂ = y): the 2D distribution (image plot on the left) and its cuts by the horizontal (graph in the center) and vertical (graph on the right) mid-planes.

Download Full Size | PDF

 

Fig. 11. X-ray degree of coherence obtained directly from the 4D CSD at the sample position of HXN beamline: (a) at y₁ = y₂ ≈ 0 vs. (x₁+ x₂)/2 and (x₁-x₂)/2 (left), and vs. (x₁-x₂)/2 at x₁ = -x₂ (right); (b) at x₁ = x₂ ≈ 0 vs. (y₁+ y₂)/2 and (y₁-y₂)/2 (left), and vs. (y₁-y₂)/2 at y₁ = -y₂ (right).

Download Full Size | PDF

5.3 Coherent mode decomposition of X-ray beam at HXN sample position

The CMD at the sample position of HXN beamline was performed on the 4D CSD illustrated in Fig. 9. Since the CSD in this case was calculated at X-ray beam waist, no analytical treatment of quadratic phase terms was necessary. The software packages used for the CMD are listed in section 4.3. Intensity distributions of several first modes, obtained by the CMD, are shown in Fig. 12. The flux carried by these modes individually, and the integrated (summed-up) flux, are shown in Fig. 13 (in graphs on the left and on the right respectively). As one can see from the graph on the right in Fig. 13, first ten coherent modes contribute most (∼98%) of the radiation flux in this high-coherence case.

 

Fig. 12. Examples of intensity distributions of coherent modes obtained by decomposition of 4D CSD at the sample position of the HXN beamline.

Download Full Size | PDF

 

Fig. 13. Spectral flux of coherent modes obtained by decomposition of 4D CSD illustrated in Fig. 9 (left), and the corresponding integrated spectral flux (right) in comparison with that of the UR beam directly propagated to the HXN beamline sample.

Download Full Size | PDF

Figure 14 shows the X-ray beam characteristics at the sample position, obtained by summing-up contributions from different numbers of major coherent modes: 2D intensity distributions (Fig. 14(a)), their 1D cuts by the horizontal and vertical mid-planes (Fig. 14(b)), and the degree of coherence in these planes (Fig. 14(c)), in comparison with that from the direct calculations. One can see that in this case of the partially coherent radiation at the beamline sample position, the radiation characteristics, calculated by accounting for contributions from only ten first coherent modes are basically undistinguishable from those of the direct calculation (that are derived from the CSD, see Figs. 9–11). This is consistent with the flux-based observations. The small number of coherent modes required for a complete description of the partially coherent radiation at the sample position of HXN beamline is typical for X-ray beamlines exploiting brightness and coherence of modern synchrotron light sources. This makes the coherent modes based representation of the partially coherent radiation at the sample position very convenient for applications related to the simulation of experiments and experimental data processing.

 

Fig. 14. X-ray beam characteristics determined as sums of contributions from different numbers of coherent modes, in comparison with direct calculations (without the CMD) at the sample position of the HXN beamline: (a) 2D intensity distributions vs. horizontal and vertical positions; (b) cuts of the intensity distributions by the horizontal (left) and vertical (right) mid-planes; (c) degree of coherence in the horizontal (left) and vertical (right) mid-planes.

Download Full Size | PDF

5.4 Testing the accuracy of propagation of coherent modes calculated with the subtraction of common quadratic phase terms

In the final part of our work, we tested the accuracy of the wave-optics based propagation of the coherent modes, determined at the position at the HXN beamline entrance (at 20.5 m from the center of undulator), using the CMD method with the subtraction of common quadratic phase terms, from the 4D CSD calculated at different numbers of mesh points. The UR intensity distributions, obtained from summing-up contributions of different numbers of these modes, were presented and discussed in section 4.4.

In our tests of propagation accuracy of these modes, we applied re-sampling of the modes’ electric fields by interpolation, using the “resizing” function available in the SRW code, to increase the mesh density to approximately the same level (required for the propagation accuracy) in each case of the sparse mesh used for the initial 4D CSD calculation. We note that the wavefront “resizing” algorithm in SRW performs an electric field interpolation after the substruction of its quadratic phase terms (whenever possible), therefore the electric fields of the modes, calculated from the CSD, even on a very sparse original mesh (such as the one with (128 × 40)² points), could still be re-sampled at an accuracy nominally sufficient for the propagation simulations. Each of the three sets of 1000 coherent modes with the original numbers of mesh points equal to 512 × 60, 256 × 40 and 128 × 40 was then propagated to the beamline sample position, where a corresponding new 4D CSD was calculated (on the mesh with (100 × 100)² points, as illustrated in Figs. 9–11). The CMD was then applied to each of the CSDs obtained this way. The results of these calculations are presented in Fig. 15. The graphs in Fig. 15(a) show the horizontal and vertical mid-plane cuts of intensity distributions of the first mode and the sum from first 3 modes, for all the considered cases of the initial meshes, in comparison with the corresponding mode intensities obtained from the direct calculations; the graphs in Fig. 15(b) show the corresponding intensity distributions obtained by summing-up contributions from all the modes obtained by the CMD at the sample position in each case. The intensity cuts in Fig. 15(b) formally include contributions from 100 coherent modes, though the contributions from only ∼10 first modes were significant in every case.

 

Fig. 15. Horizontal (left) and vertical (right) mid-plane cuts of X-ray intensity distributions at the HXN beamline sample position, determined from different numbers of coherent modes of the 4D CSD calculated from the propagated coherent modes that were obtained by CMD of the initial 4D CSD at the beamline entrance with the subtraction of common quadratic phase terms: (a) contributions from 1 and 3 final coherent modes; (b) contributions from all (100) final coherent modes. Calculations of the initial CSD were done on sparse meshes with numbers of points varying from 512 (h) x 60 (v) down to 128 (h) x 40 (v), see Fig. 7.

Download Full Size | PDF

As one can see from Fig. 15, the X-ray intensity distributions at the beamline sample position, obtained from propagation of the coherent modes which were produced by CMD with the subtraction of common quadratic phase terms from 4D CSD calculated initially on sparse meshes, are almost undistinguishable from the intensity distribution obtained by the direct calculation [7]. The variation of peak intensity among the considered cases is within 2.5%, which is comparable with typical error of wavefront propagation calculations, and is sufficient for optical design of beamlines, simulations of experiments, and many other applications of this type of calculations. The largest deviation from the direct calculation occurs for the set of modes calculated from the CSD on the initial mesh with the smallest, (128 × 40)², number of points. Further reduction of the number of points in the mesh used for the initial 4D CSD calculation results in more significant errors in the first CMD and in the propagation of the obtained coherent modes. However, even the use of (128 × 40)² mesh points means a reduction of the memory occupied by the 4D CSD from several TB down to ∼210 MB, and the CMD time to solution from many hours on a large computer cluster to seconds on one server.

6. Conclusion

We presented the examples of the CMD application towards the 4D CSD of the partially coherent hard X-ray undulator radiation at the HXN beamline at NSLS-II, with moderate computational resource requirements. The method that we used for the reduction of computational resources consists in the subtraction of leading common quadratic phase terms from the 4D CSD, performing CMD on the modified CSD, and restoring these phase terms in the obtained coherent modes. Since the 4D CSD after the subtraction of these phase terms doesn’t have fast oscillations vs. mutually conjugated horizontal and vertical coordinates, it doesn’t require a dense mesh for its accurate numerical representation. We have shown that with this method the memory requirements for the 4D CSD drop from the level of terabytes down to the level of hundreds of megabytes, and the CMD solving time reduces from many hours down to minutes or even seconds. This dramatically increases the practical feasibility of the CMD based approach to wave-optics simulations with partially coherent synchrotron (undulator) radiation in modern light sources, where the number of fully coherent wavefront propagations can be reduced thanks to the CMD by factor of 100 or 1000. The CMD method with the subtraction of common quadratic phase terms is implemented in SRW code, that was used for all simulations presented in this paper. The new method greatly complements two other methods of partially coherent SR calculations that were implemented in this code earlier – the direct simulation method for propagation of fully coherent SR wavefronts emitted by electrons with different initial conditions [7] and the method of performing CMD on 4D CSD at or near a radiation beam waist [22]. Together with the extensive library of fully coherent numerical “propagators” for different types of optical elements and some types of experimental samples, the updated methods for partially coherent simulations further advance SRW as a powerful wave-optics tool for various applications at modern light sources, including developments of new X-ray beamlines and simulations of experiments.