The fractional Fourier transform as a simulation tool for lens-based X-ray microscopy

The fractional Fourier transform (FrFT) is introduced as a tool for numerical simulations of X-ray wavefront propagation. By removing the strict sampling requirements encountered in typical Fourier optics, simulations using the FrFT can be carried out with much decreased detail, allowing, for example, on-line simulation during experiments. Moreover, the additive index property of the FrFT allows the propagation through multiple optical components to be simulated in a single step, which is particularly useful for compound refractive lenses (CRLs). It is shown that it is possible to model the attenuation from the entire CRL using one or two effective apertures without loss of accuracy, greatly accelerating simulations involving CRLs. To demonstrate the applicability and accuracy of the FrFT, the imaging resolution of a CRL-based imaging system is estimated, and the FrFT approach is shown to be significantly more precise than comparable approaches using geometrical optics. Secondly, it is shown that extensive FrFT simulations of complex systems involving coherence and/or non-monochromatic sources can be carried out in minutes. Specifically, the chromatic aberrations as a function of source bandwidth are estimated, and it is found that the geometric optics greatly overestimates the aberration for energy bandwidths of around 1%.

1. Introduction

The introduction of compound refractive lenses (CRLs) (Snigirev et al., 1996; Vaughan et al., 2011) has extended full-field X-ray microscopy to X-ray energies above 15 keV (Lengeler et al., 1999). With a numerical aperture of order 1 mrad, CRL-based objectives are well matched to the high brilliance of synchrotron beams. A range of methodologies have been developed: magnified bright-field imaging (Lengeler et al., 1999), Zernike contrast microscopy (Falch et al., 2018), high-resolution microscopy for imaging colloidal aggregates (Bosak et al., 2010) and dark-field microscopy, where orientation and strains of deeply embedded grains or domains are mapped in three dimensions (Simons et al., 2015, 2016). At the same time, direct space imaging can be complemented by diffraction in the back focal plane (Bosak et al., 2010; Ershov et al., 2013).

As for any microscope, it is important to be able to specify optical properties with high accuracy. Owing to a favorable ratio between wavelength and the aperture of the CRL, a combination of geometrical optics and the Abbe diffraction limit is believed to be adequate for describing simple imaging set-ups. As such, ray transfer matrix (RTM) formalisms have been developed to describe CRLs (Protopopov & Valiev, 1998; Pantell et al., 2003; Poulsen & Poulsen, 2014), leading to exact analytical solutions for the numerical aperture, vignetting, chromatic aberration and resolution for the general thick-lens case (Simons et al., 2017). This work has since been extended to describe the reciprocal-space resolution for dark-field microscopy and to identify schemes for sampling six-dimensional direct space-momentum space (Poulsen et al., 2017).

Nevertheless, there is a need for a corresponding tool based on wave propagation: to validate the RTM work, to be used in more complex set-ups where the effect of diffraction is more pronounced, and in particular to enable forward propagation of the beam in cases where phase contrast and/or coherent scattering are relevant. Numerical simulations are scarce (Schroer & Lengeler, 2005; Osterhoff et al., 2013), and their use is typically constrained to one-dimensional (1D) or very small two-dimensional (2D) areas due to limited computational resources and time. Analytical approaches are more computationally efficient, but exist only for certain optical configurations (Kohn, 2003).

This paper presents the fractional Fourier transform (FrFT) (Namias, 1980; Almeida, 1994; Ozaktas et al., 2001) as a general wave-propagation tool to be used with X-ray microscopy. The advantages of using FrFT for cascades of axially centered lenses are well known in the context of visible-light optics (Ozaktas et al., 2001; Ozaktas & Mendlovic, 1995). More recently, its use for simulating free space propagation of X-rays has been presented by Le Bolloc'h et al. (Le Bolloc'h et al., 2012; Mas et al., 1999). Performing an FrFT is equivalent to solving the Fresnel diffraction integral, but the FrFT is associated with favorable mathematical properties, such as being additive in index. This implies that the complete FrFT of a microscope setup involving numerous lenses and slits can be described in terms of FrFTs of the individual components.

We begin by summarizing the properties of the FrFT and providing a robust and fast algorithm for simulating CRL-based microscopy. Next, we demonstrate the applicability by addressing two topics of interest to X-ray microscopy, with parameters from the dedicated microscopy instrument at beamline ID06 at ESRF. We assess the resolution of an imaging system with a CRL objective combined with a slit in the back focal plane, and estimate the chromatic aberrations when using a large-bandwidth X-ray beam (pink beam).

2. The fractional Fourier transform

The FrFT is a generalization of the conventional Fourier transform. Whereas the conventional FT transforms between real and momentum space, the fractional Fourier transform can be interpreted as a continuous rotation in the real-momentum phase space.

2.1. Definition of the FrFT

The FrFT of order a, , of the function f(x) is defined as follows (Ozaktas & Mendlovic, 1995),

where is the sign of and

The transformation kernel K_a is defined for . For a = 0 the kernel is simply and at a = it is . When a = 1, the usual Fourier transform is recovered and a = −1 is the inverse Fourier transform. The kernel is periodic in a with a period of 4, so any transform order can be mapped onto the definition interval. Intermediate values of a are interpreted as rotations in phase space. One of its most important properties is the additive nature of the transform order, = , i.e. that carrying out two consecutive transforms of order a₁ and a₂ is equal to a single transform of order a₁ + a₂. We will revisit this property later.

2.2. Relation to the Fresnel diffraction integral

The Fresnel diffraction integral (Goodman, 2005) can be described in terms of the fractional Fourier formulation. Let p_1,2(x_1,2) describe the wavefields in two planes perpendicular to the optical axis, with free space propagation over the distance d between these planes. With reference to Fig. 1, consider the free space propagation at a distance d from plane p₁ to p₂. Within the paraxial approximation we have (with λ being the wavelength) (Goodman, 2005)

Mapping this onto the FrFT requires quadratic phase factors at the initial and final planes. Intuitively, this construction may be interpreted as propagation between two parabolic surfaces, q₁ and q₂, with radius of curvature at the apex of R₁ and R₂, respectively (cf. Fig. 1). The parabolic and flat simulation surfaces relate to each other as follows,

As shown by Ozaktas & Mendlovic (1995), free space propagation between two surfaces can be calculated by an FrFT. To describe this propagation we introduce the scaled variables

By inspection of equations (1) and (3) it appears that the FrFT describes the propagation if and only if the following three equations are fulfilled (Ozaktas & Mendlovic, 1995),

From these equations we can see that g₁g₂ = , from which it follows that . We can then formulate the following equation for calculating the propagation of the electric field,

To recap, we now have two parameters describing the source plane, R₁ and s₁, two parameters describing the detector plane, R₂ and s₂, as well as the transform order a, all of which in principle can be chosen arbitrarily. This gives a total of five parameters with only three equations to satisfy. Hence, infinitely many combinations may be used to describe the exact same propagation. None of these parameters therefore have a direct physical meaning, except for integer values of a in some cases (see §3.1). As such the FrFT method is completely equivalent to normal Fourier optics, and the process of calculating the free space propagation is identical in both cases: multiply by a quadratic phase, perform a (fractional) Fourier transform, and multiply by a second quadratic phase (Goodman, 2005). By setting the transform parameter a = 1 and the object scaling s₁ = 1, and solving equation (7) one will reach the normal Fourier optics parameters.

Figure 1 The fractional Fourier transform can be seen as a free space propagation between two planes p1 and p2 by using quadratic phase factors, as indicated by the two parabolic surfaces q1 and q2. See the main text. In this case the sign of the two curvatures are and .

2.3. Implementation for direct space propagation

In the following, we shall resolve the ambiguity by defining the source plane parameters R₁ and s₁. In this case, the image plane parameters and the transform order are given by

The choice of scale parameter s₁ affects the sampling in both real and momentum space. The numerical implementation used in this work is from Ozaktas et al. (1996), and assumes that the scaled real space and scaled momentum space are equally long, i.e. real-momentum phase space is square. This makes the computation of the FrFT easier, but restricts the value of s₁ to be [see Ozaktas et al. (1996) for details]

where Δx₁ is the field of view of the simulation and N is the number of pixels. Using the scaled variables has the benefit that the detector space will automatically be scaled by a factor s₂/s₁, removing any aliasing as the detector plane will be large enough to accommodate all scattered intensity. Having a fixed scaling parameter for the numerical simulations only leaves R₁ as a free parameter.

From equation (4a) it is seen that the quadratic phase term is inversely proportional to the wavelength. For the short wavelength of X-rays the phase will change rapidly, even for small objects, which requires 10⁵–10⁶ samples per dimension for tens to hundreds of micrometer-sized objects (Schroer & Lengeler, 2005; Osterhoff et al., 2013). Here we show that the FrFT approach allows the sampling requirements to be relaxed such that numerical simulations of large objects in 2D become feasible with economical computer hardware. This is achieved by choosing R₁ = for an incoming plane wave, which eliminates the quadratic phase. Thus, all sampling requirements are removed, except that the object must be sampled properly (which is only a concern if the object has a rapidly changing phase). If instead the X-ray source is converging or diverging, for example by placing the object before or after the focal point of a condenser lens, we can choose R₁ to be either negative or positive to match the curvature of the source wavefront. This again removes the quadratic phase factor in the object plane.

With s₁ = and choosing R₁ = , the transform order is

At large propagation distances, a is approaching a value of 1, corresponding to a usual Fourier transform as we go into the Fraunhofer diffraction regime. A note of caution here is that the transform order is implicitly a function of the number of pixels in the source plane through s₁, and a comparison of the a parameter along the propagation path is only valid if R₁ and s₁ are not changed.

This approach does not seek to remove the phase curvature in the detection plane, and R₂ may very well be finite for a given propagation. Therefore, one should consider the sampling in the detection plane only if the phase is required (the intensity is not affected by the phase).

3. Application to microscopy based on CRLs

For a simple free space propagation, the FrFT approach provides two major benefits: first, by removing the strict sampling constraints of standard Fourier optics, and, second, by avoiding aliasing through the automatic scaling of s₁ and s₂. Furthermore, when analyzing systems containing numerous lenses, we can capitalize on the additive nature of the transform order to speed up calculations drastically.

3.1. Lenses and the FrFT

In the approach below, we first of all take advantage of the fact that with classical refractive optics the numerical aperture (NA) is limited to approximately NA = , where δ is the refractive index decrement (Als-Nielsen & McMorrow, 2011). With δ of the order of 10⁻⁶, the NA is of the order of 3 mrad, i.e. the paraxial approximation is fulfilled. Next, we shall assume that the lens shape is parabolic, and will treat each lens as an infinite and perfect thin lens, therefore implying that its effect on the wave propagation is only a shift in the phase (Goodman, 2005). Analysis using geometrical optics shows that for the X-ray energies under consideration here the focal length, f, of a single refractive lens element (lenslet) is much larger than its thickness, T (Simons et al., 2017). The resulting complex transmission function, t_lens, depends on the distance to the optical axis as follows [Goodman (2005); equations (5)–(10)],

where f is the focal length of the lens. The validity of the latter approximation may be experiment-specific, and we shall evaluate it for a particular setting below.

Fig. 2 depicts a typical setup with propagation from a flat object plane to a lens, phase shift within the lens, and propagation to a flat detector plane. The relevant FrFT parameters are indicated.

Figure 2 Schematic showing an optical system with a single thin lens. The parameters associated with the FrFT have been defined. See the main text.

The wavefront in the detector plane can be described as follows using the FrFT,

For the first propagation we set R₁^- and s₁^-, and so R₁⁺ will be given. For the second propagation, we specify R₂^- such that the exponential terms vanish,

As a part of the thin lens approximation we use the paraxial approximation, meaning that s₂⁺ = s₁^-. With values for R₂^- and s₂^- defined from these two relations, the remaining parameters are given. The propagation now takes the following form (using the additive transform order property and R₁^- = ),

This means that we can perform the entire wavefront propagation using a single transformation. Furthermore, this holds for arbitrary positions of the object and detector planes, i.e. the detector can be placed at any point, not only in focal planes or in imaging planes.

This FrFT approach does not take any absorption into account, so the lenses have an infinitely large pupil. Therefore, to take absorption into account the simulation has to be stopped at each plane where absorption occurs. We will discuss this in terms of CRLs in §3.2.

The concept of cascading the FrFTs can be straightforwardly extended to a system with any number of lenses (which can have varying distances and varying focal lengths, and negative focal lengths as well). In such a case the FrFT parameters (R_i^-, s_i^-, a_i, R_i⁺ and s_i⁺) may be calculated iteratively, based on initial choices of R₁^- and s₁^-,

Equation (8) assumes positive s_i ^+/-, which is not guaranteed from these equations. Changing the sign of s is equivalent to inverting the image, so a negative s_i⁺ may have its sign changing by adding 2 to the corresponding a_i, which also performs an inversion. So mathematically : (s_i⁺ = = a_i+2).

Fig. 3 shows a numerical example of the curvature of the intermediate planes involved in a CRL simulation using FrFT. For this example, we used a small number of lenses with a short focal length to make the figure easier to interpret. We used N = 7 lenses each with a focal length of f = 2.12 m, giving an effective focal length of the CRL of f_N = 26.5 cm. The blue and red surfaces correspond to R_i ^- and R_i ⁺, respectively. In this example, the object and detector have been placed such that an M = 2 times magnified image is depicted on the detector. The object field-of-view in this example is = 100 µm, and the size of the surfaces are proportional to s_i⁺/s_i^-, the width of the simulated area (automatically adjusted so that all diffracted intensity stays within the field-of-view). As we have used an imaging geometry with a magnification of 2, the detector plane field-of-view is twice as large as the object plane field-of-view. Also note there is a quadratic phase in the image plane, which has no impact on the measured intensity but has to be taken into account for further propagations after the image plane.

Figure 3 Illustration of the curved surfaces appearing in an FrFT simulation of CRLs, in this case for a CRL with seven lenslets and a setting with a magnification of M = 2. The blue line to the left is the sample plane and the red curve to the right represents the image plane. The blue and red surfaces correspond to Ri^- and Ri^+, respectively. The sizes of the planes are 1000 times larger than the simulated area to show the curvatures more clearly, and they are proportional to si^+/si^- to show the extent of the diffracted intensity. The focal points of the CRL are indicated as well.

For a given optical setup, the transform order a can be calculated at any plane, where a is an odd or even integer corresponding to that plane being a Fourier or imaging plane, respectively. This is true for any choice of scaling parameters and sampling, as long as R₁ ^- = . In the example in Fig. 3, = 2, and the detector therefore records an inverted image of the object. The magnification is given by the ratio s_detector⁺/s_object^- = 2.

3.2. Effective pupil function for a CRL

While the FrFT provides an elegant way of propagating the wavefront through a cascade of non-absorbing lenses, the question remains of how to handle the attenuation. Can we find an effective pupil function to apply at a single plane, or is it necessary to treat the absorption at each physical lens? In the latter case, the simulation requires propagation of the electric field from lens to lens, followed by a manual reduction of the amplitude (the square-root of the intensity). This raises two issues for practical simulations: the numerical stability of the propagation method and a significant increase in execution time.

In the following, we shall investigate (a) the numerical stability of the slow approach (one lens at a time), (b) the possibility of defining an effective aperture that represents the wavefield propagation coming from a point in the sample plane and on the optical axis (simulation of point spread function) and (c) the possibility of defining another effective aperture that represents the vignetting.

For these numerical tests, we will use a MatLab implementation (Ozaktas et al., 1996; Bultheel & Martíez Sulbaran, 2004) on a workstation PC and use double-precision complex numbers (128-bit floating point). The experimental parameters are imported from an existing experimental setup at ID06 at ESRF. The CRL objective lens consists of N = 69 identical Be lenses, each with a curvature of R = 50 µm and a center-to-center distance of T = 1.6 mm. Using a photon energy of E = 17 keV the refractive index decrement is = 2.36×10^-6 and the linear attenuation coefficient is = 40.7 m⁻¹, giving an individual lens focal length of f = 21.2 m, a combined focal length of f_N = 0.270 m and a numerical aperture of NA = 3.61×10^-3 (Simons et al., 2017). The CRL objective lens and detector positions have been chosen to form a real image with a magnification of M = 10, corresponding to an object to entry-of-lens distance d₁ = 0.302 m and an exit-of-lens to detector distance d₂ = 3.54 m [equations (18) and (19) of Simons et al. (2017)]. Using equation (16) and performing the FrFT propagation we find the image plane to be 0.46 mm closer to the CRL. For the physical propagation this small difference is insignificant, but the origin of the difference is unknown. The underlying assumptions for the RTM formalism and the FrFT propagation are the same, and the error seems be too large for numerical errors. For the simulations we shall assume that the lenses are infinitely large, not taking the physical aperture into account.

3.2.1. Numerical stability

In terms of stability, we find the maximum intensity errors after 69 individual propagations compared with a single propagation to be 0.1%–0.3%. The accuracy is not good at a few pixels distance from the edges; however, this is irrelevant in practice as some space should always be left blank at the edges (i.e. a support). The phase deviation is less than 1% in the areas with significant intensity, but is unreliable in the dark regions (where the phase is ill-defined anyway). See comparison in Fig. 4.

Figure 4 The point spread function simulated with physical attenuation at each lens (lens-by-lens simulation), compared with using an effective aperture placed at the last lens in the CRL stack. The figure shows the intensity (a), intensity error (b), phase (c) and phase error (d). We see a small phase difference between the two procedures, which is equivalent to the focusing point being moved about 150 nm away from the sample, and is likely due to numerical errors accumulating in the stepwise simulation. Simulation details: 2D wavefront, 1000 × 1000 pixels, 1 nm pixel size; see also `Example_Figure_4.m' of Pedersen (2017). The plots are cross sections perpendicular to the optical axis.

3.2.2. Effective pupil function in the case of negligible vignetting

In the following, we consider propagating the scattering from a 1 nm object on the optical axis, a size that is much smaller than the point spread function (PSF). Fig. 4(a) shows the resulting image intensity on a detector in the image plane, i.e. the PSF, in relation to the object plane coordinates. The figures show the results of two simulation runs, one using the physical attenuation at each lens (lens-by-lens simulation, 70 propagation steps) and the other using an effective pupil at the last lens in the CRL stack (two propagation steps). The RMS width of the Gaussian pupil function is calculated analytically using the RTM formalism,

where is the RMS acceptance angle of the CRL, given by Simons et al. (2017), and = (T/f )^1/2 [see the supporting information (SI) for the derivation]. As the pupil function is Gaussian, the PSF is Gaussian as well. Fig. 4(b) shows the absolute and relative intensity error of the effective aperture compared with the lens-by-lens simulation. The differences are negligible. Similarly, Figs. 4(c) and 4(d) show the phase and phase error of these two simulations, again showing negligible differences. The sampling is not important in calculating the PSF as long as the pixel size is smaller than the PSF. As shown in Fig. S1 in the SI, the intensity is indeed independent of sampling. This is also true for larger and more complex objects (see Fig. S2 in the SI). We can also model the intensity and phase in the vicinity of the image plane (see the intensity and phase maps in Fig. S3).

3.2.3. Effective pupil function in the case of vignetting

To include the effect of vignetting as well, we propose to supplement the effective aperture (positioned in the plane of the last lens) found above with a Gaussian attenuation in the object plane. We do not have a known expression for the width of this attenuation profile a priori, but it can be derived numerically from the rotational symmetry of the CRL within fractions of seconds (see the SI for details).

To test this proposal, we made simulations for a uniform incoming plane wave extending the entire 1 mm field-of-view of the object plane. Here we compare the propagation when using lens-by-lens attenuation (70 propagations) and the aforementioned combination of the effective aperture at the exit of the CRL (fixed, given by geometrical optics) and a Gaussian attenuation in the object plane (variable width) (two propagations). Fig. 5 shows the vignetting function (relative to the object plane) obtained with the two simulations for the optimal width of the aperture in the sample plane of 185 µm. The intensity and intensity error can be seen in Figs. 5(a) and 5(b), respectively, and the phase and phase error are shown in Figs. 5(c) and 5(d), respectively. The difference is seen in all cases to be negligible.

Figure 5 Comparison of the vignetting function when using lens-by-lens attenuation and using a combination of two effective apertures, one in the sample plane and one at the exit plane of the CRL. The intensity relative to the object plane is shown in (a) and the corresponding error in (b). The phase and phase error are shown in (c) and (d), respectively. Simulation details: 2D wavefront, 1000 × 1000 pixels, 1 µm pixel size; see also `Example_Figure_5.m' of Pedersen (2017). The plots are cross sections perpendicular to the optical axis.

As an example of the effect of vignetting, we propagate a 2D image from the object plane to a 10× magnified imaging plane, using both the effective vignetting and pupil as well as lens-by-lens simulations. The results are given in Fig. 6, and as expected from the results above the images with effective attenuation are indistinguishable from the much slower lens-by-lens simulation. Note the object size difference between Figs. 5(a) and 5(d), and also notice that the images are inverted and magnified ten times (as expected).

Figure 6 Example of vignetting. A small (a) and a large (d) test object is propagated to the imaging plane with a magnification of 10. Attenuation has been implemented lens by lens [(b) and (e), 70 propagations] and with an effective vignetting and pupil function [(c) and (f), two propagations]. Simulation details: 2D wavefront, 1200 × 1500 pixels, 5 nm [(a), (b) and (c)] and 400 nm [(d), (e) and (f)] pixel size; see also `Example_Figure_6abc.m' and `Example_Figure_6def.m' of Pedersen (2017).

4. Examples of use

We anticipate that the FrFT can be a powerful tool for optimizing entire microscopy beamlines comprising CRLs as pre-condensers, main condensers and objectives. For reasons of simplicity, below we present two examples of use, involving only an objective. The examples address issues of current research and refer to the ID06 setup presented in §3.2.

4.1. Aperture in the back focal plane: effect on spatial resolution

In an imaging setup, by using a CRL objective one can access both the real space imaging plane and the Fourier transform in the back focal plane (BFP). This provides a number of benefits: as an example, filtering the signal with an aperture in the BFP allows imaging of regions with different strain separately. However, introducing a small aperture in the Fourier plane will reduce the obtainable real space resolution. As the CRL objective also limits the resolution, the question becomes: what is the effective PSF of this combined system? To answer the question, we simulated the intensity from a single point through the imaging setup described above, with the only addition of a square aperture in the center of the BFP with variable side length D. In these simulations, we obtain a nearly Gaussian PSF, and the resolutions reported in this section are the RMS widths of the Gaussian PSFs. Furthermore, we compare these simulated resolutions with what we would expect from geometrical optics. The resolutions reported here are relative to the object plane, i.e. the magnification has been taken into account.

Within geometrical optics the resolution from the CRL alone can be derived from the PSF [equation (31) of Simons et al. (2017)¹],

The resolution for a slit in the BFP alone is (see derivation in SI)

From geometric optics, we would expect the effective resolution simply to be a geometric sum of the two contributions,

Fig. 7 shows these analytical resolutions and the fitted resolution of the FrFT simulations as a function of the aperture size. As expected, the analytical and simulated resolutions agree very well when the resolution is dominated by either of the two components. Nevertheless, in the intermediate region the simulated resolution deviates from the geometrical optics result with . The error on the fit of the simulated PSF is of the order of the thickness of the line. The PSF shape inherently relies on diffractive effects, and we therefore trust the FrFT wavefront calculations. This example highlights the fact that the wavefront simulations are more accurate than a simple geometrical optics approach, even for rather basic experimental setups.

Figure 7 Resolution of an imaging system with a combined CRL objective and a square aperture in the back focal plane. Shown are results from geometrical optics, , and from FrFT simulations, , as a function of aperture size. In the intermediate region where the two components contribute equally to the resolution, geometrical optics fails at predicting the resolution. Simulation details: 2D wavefront, 1000 × 1000 pixels, 10 nm pixel size, PSF widths extracted from Gaussian fitting to vertical line profiles through the optical axis; see also `Example_Figure_7.m' of Pedersen (2017).

4.2. Chromatic aberration for pink-beam operation

Owing to the relatively high speed at which optical propagations can be simulated using the FrFT, it becomes feasible to perform extensive wavefront calculations of partially coherent beams. The following example focuses on partial longitudinal coherence due to an energy spread of the source. Determining chromatic aberrations due to the incoming beam energy spread is of great interest, in view of plans to operate with so-called pink beams with bandwidth of up to 1% (Falch et al., 2016).

The procedure of calculating propagation with partial longitudinal coherence is detailed by Voelz (2011; ch. 9.1). In short, the idea is to calculate the fully coherent image intensity at different adequately spaced wavelengths, and then to calculate a weighted sum of these images according to the energy distribution in the source. In this section, we will assume the energy distribution in the source to be Gaussian with an RMS width between = 10^-4 (typical value for a crystal monochromator) and = 10^-2 (typical width for a pink beam). We will assume a central energy of E₀ = 17 keV. As we will be testing many different energy distributions, we perform a large amount of fully coherent simulations within [following Voelz (2011)] an energy spread of . The step size was set at 1 eV; inspection showed that further subdivision leads to no changes in the final image. This gives a total of 1021 energy steps, which for a 1000 by 1000 pixel image can be calculated in a few minutes on a modern workstation PC. In this simulation we used an effective aperture to speed up the simulation time.

The test object for this simulation is a horizontal line, which allows us to extract the PSF with a vertical profile through the optical axis and compare it with geometric optics results. The setup is the same as mentioned above, with the distances optimized for a magnification of ten with an energy of E₀ = 17 keV. For these simulations, the geometrical distances are kept constant, whereas the material parameters, and thus the focal length and absorption, are updated for each photon energy step (Henke et al., 1993). All the PSFs are corrected for the magnification, and refer to the position in the object plane.

The fully coherent simulation step at E₀ = 17 keV gives the same Gaussian PSF as seen for the CRL alone in §4.1 (see Fig. 8). The other photon energies also give a Gaussian beam shape, although much larger, as they are not in perfect focus. When adding the fully coherent simulated intensities according to the bandwidth, the resulting PSF shapes become progressively non-Gaussian for larger bandwidths, as seen in Fig. 8. The chromatic aberrations from a crystal monochromator ( = 10^-4) are negligible, but for the larger bandwidths the main peak broadens and develops significant intensity outside the main peak.

Figure 8 The PSF for selected energy bandwidths. The crystal monochromator with has negligible broadening over the ideal Gaussian PSF, whereas the broader bandwidths have broader peaks and become increasingly non-Gaussian. Simulation details: 1D wavefront, 1000 pixels, 10 nm pixel size; see also `Example_Figure_8_9_1D.m' of Pedersen (2017).

A full overview of the PSF as a function of energy bandwidth is seen in Fig. 9, in which vertical slices are PSFs corresponding to the bandwidth on the x-axis. Fig. 9(a) shows the simulated PSFs and Fig. 9(b) shows the PSFs from geometrical optics. The geometrical optics PSF from the chromatic aberration alone is calculated from equation (38) of Simons et al. (2017), and then convoluted with the CRL resolution PSF to give an effective PSF. The energy distributions used to weight the sums of the fully coherent simulations are not normalized, and so the increase in PSF intensity reflects the increase in flux. The black lines show the 25% (full width at quarter-maximum), 50% (full width at half-maximum) and 75% (full width at three-quarters-maximum) of the relative intensity of each PSF, thus giving a sense of the resolution. A striking difference is seen between the FrFT simulation and geometrical optics: according to the geometric optics the resolution should simply keep deteriorating with increasing bandwidth, whereas the FrFT simulations show that the resolution only decays rapidly up to about = 10^-3, after which the resolution only changes slightly.

Figure 9 The point spread function as a function of energy bandwidth resulting from (a) FrFT simulations and (b) closed expressions from geometrical optics. For each bandwidth, the intensity is shown as a function of position in the image plane. The color represents the intensity on an absolute scale. Hence, the increase in the intensity at the central position with larger bandwidth reflects an increase in incoming photon flux. The black lines are the 25%, 50% and 75% relative intensity lines. Simulation details: 1D wavefront, 1000 pixels, 10 nm pixel size; see also `Example_Figure_8_9_1D.m' of Pedersen (2017).

Why are the two cases so different? The FrFT approach includes diffractive phenomena, which causes the PSF to defocus for the non-optimized wavelengths, in effect reducing the influence of the wavelengths far from the central wavelength, as the intensity will be substantially smeared out. This example highlights the benefits of wavefront propagation methods, which include diffractive phenomena in propagation. The PSF discussed here relates only to the on-axis resolution, and further simulations are required for more detailed off-axis resolutions.

5. Discussion

As mentioned earlier, we assume each individual lens to be a thin lens. According to Goodman (2005) the typical thin-lens approximation for optical lenses involves two approximations: that the lens surface is described by a paraboloid and that the rays are parallel to the optical axis inside the lens, i.e. the paraxial approximation. These approximations hold well in the case of CRLs, since the true lens profile is parabolic and the angular deviations within the lens stack are typically very small. In addition, the focal length of the individual lenses is much longer than the lens thickness.

As demonstrated, the FrFT approach has a great benefit over conventional Fourier optics by removing the quadratic phase term in the object plane, which for X-rays requires a very high sampling rate. Furthermore, it speeds up CRL calculations, as all the free space propagations up to the CRL exit plane can be combined into a single transform. The absorption of the CRL can be handled by an effective vignetting function (applied in the object plane) as well as an effective pupil function (applied at the end of the CRL). The effective vignetting and pupil functions have been demonstrated here mostly in imaging geometries, but the approach is valid in general cases as well [see, for example, `Example_7_1D_CRL_Condenser.m' of Pedersen (2017), in which the object–lens distance is 0]. The effective vignetting and pupil functions have been calculated assuming infinitely large lenses, but we anticipate that it is possible to include finite apertures in the same way with an effective pupil and vignetting function by using a combination of Gaussian and box functions.

The simulations presented here have been based on perfect lenses. The effect of imperfect lenses can also be modeled if the effects can be expressed in terms of a complex transmission function (attenuation and phase shift), relative to the ideal lens. If the imperfections of the lenses can be modeled cumulatively, as discussed by, for example, Seiboth et al. (2017), the complex transmission function may be applied at the CRL exit plane along with the effective pupil function. Alternatively, one can apply transmission functions at each lens, but then having to perform N + 1 propagation steps.

The FrFT approach is a powerful simulation tool in experiments involving coherent effects. One example is ptychography, where the probe intensity, shape and phase may be simulated. If the sample interaction can be modeled, the diffraction can be simulated as well, so one can perform the entire experiment in simulations. One can, for example, test the effect of different noise levels and probe displacement on the reconstruction quality. Another example is phase contrast microscopy, in which the diverging beam may be taken into account by setting R₁ ^- so that the phase curvature matches that of the incoming beam. Again, the quadratic phase factor is canceled out, so large 2D areas may be simulated on modest computer hardware. A last example involving CRLs is Zernike phase contrast, in which one can experiment with putting the phase shifter in different locations, either inside the CRL lens stack or in the back focal plane.

Finally, FrFT provides a seamless and computationally fast way to tie together diffraction patterns from various values a, representing, for example, a Fourier plane (far-field Fraunhofer limit or back focal plane in an imaging system), direct space (imaging plane) or any other rotation in direct-Fourier space. The FrFT approach is useful in any wavefront simulation situation, as long as the optical components can be modeled by a complex transmission function. Other types of lenses, such as zone plates and multi-layer Laue lenses, effectively add a quadratic phase, and can therefore be modeled as a simple lens in the FrFT framework. More complex examples are flat mirrors, focusing KB mirror systems and crystal optics. Even free space propagation may benefit from the FrFT approach if the field of view is large.

6. Conclusion

In this paper we have adapted the general FrFT wavefront propagation method to easy-to-use Matlab code for wavefront propagation of X-rays. We have leveraged the extra degrees of freedom associated with the FrFT propagation method, as compared with traditional Fourier optics, to remove the quadratic phase factor in the object plane. This removes the strict sampling requirements typically encountered for X-rays due to the short wavelength, and opens up for large field-of-view simulations in both 1D and 2D.

The FrFT method intrinsically describes fully coherent propagation, but can also be used for partial transverse and longitudinal coherence. Partial coherence is implemented by averaging repeated coherent intensity simulations, and here the removal of the sampling requirements is very helpful, as the otherwise time-consuming simulation can be greatly speeded up and be performed on modest computer hardware.

The FrFT propagation technique is very general, and can describe the propagation through highly complex optical setups. The only requirement is that the effect of the optical element can be described by a complex transmission function. Lenses can be treated implicitly in the FrFT propagation, and a single transform can be used to propagate through any number of lenses, given that they have an infinite pupil. To incorporate any type of aperture, either the pupil function of a lens or a standalone aperture, the attenuation is applied at the plane of the optical element. Here, we have shown that in the case of CRLs the attenuation in each lens can be described by an effective pupil function applied at the end of the CRL and an effective vignetting function applied at the object plane. In this way, propagation simulation through a CRL is no more computationally intense than a single lens.

As numerical examples, we, first, calculate the PSF of an imaging system with a CRL objective lens and a variable-sized square aperture in the back focal plane, and, second, the PSF of an imaging system with a CRL objective lens using pink-beam illumination with variable bandwidth. In both cases the FrFT wavefront propagation provides a more accurate result compared with analytical expressions derived from geometrical optics.

Footnotes

¹The expression given by Simons et al. (2017) should be multiplied by a factor .