 1. Introduction
 2. Methods for describing partial coherent beams from undulators in a storage ring
 3. Propagation of wavefronts along the beamline
 4. Description of the optical system
 5. Results and discussion of multioptics simulations
 6. Summary and conclusions
 7. Related literature
 Supporting information
 References
 1. Introduction
 2. Methods for describing partial coherent beams from undulators in a storage ring
 3. Propagation of wavefronts along the beamline
 4. Description of the optical system
 5. Results and discussion of multioptics simulations
 6. Summary and conclusions
 7. Related literature
 Supporting information
 References
research papers
A fast and lightweight tool for partially coherent beamline simulations in fourthgeneration storage rings based on coherent mode decomposition
^{a}ESRF – The European Synchrotron, 71 Avenue des Martyrs, 38000 Grenoble, France
^{*}Correspondence email: srio@esrf.eu
A new algorithm to perform coherent mode decomposition of undulator radiation is proposed. It is based on separating the horizontal and vertical directions, reducing the problem by working with onedimension wavefronts. The validity conditions of this approximation are discussed. Simulations require low computer resources and run interactively on a laptop. The focusing with lenses of the radiation emitted by an undulator in a fourthgeneration storage ring (EBSESRF) is studied. Results are compared against multiple optics packages implementing a variety of methods for dealing with partial coherence: full twodimension coherent mode decomposition, Monte Carlo combination of wavefronts from electrons entering the undulator with different initial conditions, and hybrid raytracing correcting geometrical optics with wave optics.
Keywords: partial coherence; undulator radiation; wave optics; beamline simulation; OASYS.
1. Introduction
The migration to fourthgeneration storage rings has significantly improved the ^{−3} to 10^{−2} at 10 keV). This has a beneficial impact^{1} for many applications requiring coherent beams, such as Xray photon correlation spectroscopy, coherent diffraction imaging, propagationbased phasecontrast imaging, and ptychography (Paganin, 2006). However, the diffraction effects produced by the interaction of the beam with the edges and distortions of the optical elements strongly affect the quality of the beam. Diffraction patterns show higher visibility due to the increased coherent fraction in new sources, and its accurate modelling is fundamental for the design and optimization of beamlines. The physical models for the limiting cases of full incoherence (usually simulated by geometrical raytracing) or by propagating a single wavefront (valid for fully coherent radiation) are not sufficient for a complete understanding of the beam transport (Sanchez del Rio et al., 2019). The coherent fraction of the radiation emitted by new generation storage rings, although much improved with respect to previous generations, is still of the order of a few percent at hard Xrays, which means that it is mandatory to account for partial coherence. In recent years, several modelling approaches have been demonstrated to work for beamlines using undulator radiation. Starting from incoherent beams, Shi et al. (2014) proposed some correction algorithms to include diffraction effects that occur with coherent radiation. More accurate methodologies exploit the well known propagation of coherent wavefronts. The partial coherence is treated by propagating a set of wavefronts that all together describe the undulator radiation. Two approaches are possible. One consists of calculating the wavefronts emitted by electrons entering the undulator with different initial conditions, sampled by Monte Carlo from the electron beam emittance (multielectron in SRW) (Chubar et al., 2011). A second method, the coherent mode decomposition (CMD), assigns these wavefronts to the coherent modes, which are the eigenfunctions of the crossspectral density (CSD), and can be numerically calculated for the undulator source (Glass & Sanchez del Rio, 2017). A full treatment of CMD with twodimensional (2D) wavefronts was implemented a few years ago in the COMSYL package (Glass, 2017a). Both methods require the use of highperformance computer (HPC) resources that are not always at hand. The problem in CMD is to manipulate and diagonalize a huge stack representing the CSD with enough precision, which is a fourdimensional (4D) entity when using 2D wavefronts. After the release of COMSYL, different techniques have been proposed to deal with the magnitude of the problem. The singlevaluedecomposition method presents some advantages when used for the diagonalization of the CSD (Xu et al., 2022). When the wavefronts are highly convergent or divergent, sufficient sampling of the electric field phase requires a very fine grid. In these cases, the sampling is dictated by a quadratic phase. A method developed by Li & Chubar (2022) consists of subtracting the quadratic phase which is analytically propagated, thus reducing the computational effort to limits acceptable by an average laptop.
and coherence of Xray synchrotron sources. The transverse coherent fraction of the new sources is increased by at least one order of magnitude with respect to the thirdgeneration sources (typically from 10In this paper, we propose a new method for dealing with partial coherence of undulator beams. The key point is to reduce the dimensionality of the problem to deal with onedimensional (1D) wavefront cuts (i.e. separating horizontal and vertical directions).
We demonstrate that whenever this 1D approximation can be used, like in many cases of practical interest, the results are comparable with the other 2D methods but require much fewer resources, thus allowing simulations using a standard laptop.
The new code, referred to here as WOFRY1D, is benchmarked against other existing codes that are available in the OASYS simulation ecosystem (Rebuffi & Sanchez del Rio, 2017). The optical system studied here derives from the project for the new EBSL1 beamline at the upgraded EBSESRF storage ring. We have compared results for different setups implementing two refractive systems (transfocators), plus a slit placed upstream of the transfocators. The beam properties simulated by four different transfocator configurations are studied in detail using four packages available in OASYS: (i) the novel 1D CMD, implemented in the code WOFRY1D, (ii) full CMD in 2D with COMSYL (Glass, 2017a), (iii) SRWME: multielectron simulations in SRW (Chubar & Elleaume, 1998), and (iv) HYBRID raytracing as described by Shi et al. (2014) and implemented in ShadowOUI (Rebuffi & Sanchez del Rio, 2016).
2. Methods for describing partial coherent beams from undulators in a storage ring
In this section, we summarize the basic theory underneath partially coherent emission from electrons in storage rings. We start by showing that a relativistic single electron emits fully coherent radiation when passing through an undulator magnetic field. We then move to the emission from relativistic electron bunches showing that an electron beam with nonnegligible emittance will produce a partially coherent emission and that a higher coherent fraction is associated with a lower electronbeam emittance. Finally, we present the basic principles underlining the numeric calculations within the packages used.
2.1. Description of undulator emission
We quickly remind that an ultrarelativistic charged particle following a curved trajectory [usually wiggly as produced by alternated magnetic fields in insertion devices (IDs)] emits radiation. This electric field can be calculated in the framework of classical electrodynamics [see, for example, equation (14.14) of Jackson (1999)]. In the frequency domain, the electric field at an observation point r = (x, y, z) can be written as
where the subscript ω indicates the photon frequency, e is the c the velocity of light, ε_{0} the electric constant, γ ≃ [GeV] is the Lorentz factor with being the electron energy in practical units, β = is the electron relative velocity and the dot denotes the time derivative. Also n(t) = r − r_{e}(t)/r − r_{e}(t) is the unit vector pointing from the particle to the observation point r; the electron trajectory is represented by r_{e}(t), which is completely determined by the 3D distribution of the magnetic field inside the ID and the electron initial conditions prior to entering it. The origin of the vector r is usually at the centre of the insertiondevice/straightsection. Figure 1 serves as a visual aid to equation (1) and its parameters.
Equation (1) describes a fully spatially coherent field and has been conceptualized for a single electron. A common abstraction that derives from it is the `filament beam', where N_{e} electrons overlap in space following the same trajectory r_{e}(t), which is useful to represent an idealized zeroemittance storage ring. In this case, a multiplicative factor N_{e} is applied to equation (1). Much like the single electron emission, the filament beam also radiates a fully transverse coherent wavefront.
Several codes are available in the synchrotron community to calculate the undulator emission characteristics in different cases. The codes URGENT (Walker & Diviacco, 1992) and US (Dejus & Luccio, 1994) compute undulator emission in the farfield for undulators with a sinusoidal magnetic field. The codes SPECTRA (Tanaka & Kitamura, 2001) and SRW (Chubar & Elleaume, 1998) are more generic as they calculate emission in the near and farfield for any electron trajectory (with different initial conditions) and submitted to an arbitrary magnetic field.
2.2. Electron beam distribution in storage rings
At any position s in the storage ring, an electron can be described by five coordinates, = (x_{e}, , y_{e}, , ), representing the phase space coordinates and a term expressing the relative deviation of the electron energy from the main storage ring energy (also known as the energy spread). It follows that at any given s the many electrons in a bunch follow a 5D Gaussian distribution,
with M as the inverse of the generalized variance 5 × 5 matrix. A common assumption is that the variables are correlated only if they are in the same plane (x or y). In some particular points s where the covariance between spatial and angle terms is zero, only the diagonal terms in M are nonzero: (, , , , ). This is usually the case at the centre of the straight sections, where the undulators are often placed. We also assume that the electron bunch has a Gaussian distribution with σ_{z} along the longitudinal direction, as most particle beams do [cf. §8 in Wiedemann (2019)].
2.3. Emission from electron bunches
Having summarized the coherent emission from a single electron in Section 2.1 and how the electrons are statistically distributed in a bunch (Section 2.2), we now turn our attention to the emission from the electron bunch with finite emittance.
The total electric field emitted from all N_{e} electrons in a bunch circulating in a storage ring is given by
In terms of intensity,
The electric field is the emission by a single electron with a trajectory defined by the undulator magnetic field and electron initial conditions at the observation point r [see equation (1)]. The first term in equation (4) is a sum at r of the intensity of every electron emission weighted by its probability f, which describes temporally incoherent synchrotron radiation (SR). The second term stands for the enhancement of the intensity due to coherent superposition of the emission of the N_{e} electrons, modelling temporally coherent synchrotron radiation. It follows that I_{bunch} = I_{iSR} + I_{cSR}. For emitted wavelengths λ shorter than the electron bunch length (σ_{z} > λ), the power associated with the term I_{cSR} vanishes quickly (Hirschmugl et al., 1991; Wiedemann, 2019). Considering typical undulator radiation emission, i.e. Xray energy ranges from a few hundred electronvolts to a few hundred keV, and typical electron bunch lengths in storage rings (σ_{T} > 30 ps), I_{cSR} can be neglected when considering standard monochromatization schemes in beamlines.^{2}
Similarly, the mutual correlation of the electric field between two observation points r_{1} and r_{2} is
where the superscript asterisk ( ^{*}) indicates the complex conjugate, the angular brackets 〈…〉 indicate the sum over the bunches, r_{1} = (x_{1}, y_{1}, z_{1}) and r_{2} = (x_{2}, y_{2}, z_{2}) (see Fig. 1). This equation is the crossspectral density, that will be discussed later and rewritten in a more manageable form.
2.4. Multielectron Monte Carlo (SRWME)
Synchrotron radiation emitted by undulators in storage rings is a fundamentally random process and should be treated probabilistically. The SRWME algorithm used to account for partial coherence implements equations (4) and (5) by individually calculating the synchrotron emission (electric field) of several electrons subjected to the initial conditions sampled from assuming these are uncorrelated and passing through an arbitrary magnetic field describing the Xray source. Each resulting electric field from this Monte Carlo sampling is then propagated through the beamline until the observation point, where the contributions from different electrons are added in intensity (Chubar et al., 2011). It is impractical (and unnecessary) to account for the emission of every single electron in a beam that very often has a current of a few hundred mA. Electrons are then divided into socalled macroelectrons (me), which is an abstraction that allows the emission of several individual electrons to be grouped into one `superparticle' emitting a fully coherent wavefront but with resulting intensity given by the total intensity I_{bunch} divided by the number of macroelectrons used in the simulation,
An advantage of the SRWME approach is that, since the electric fields of the macroelectrons propagate independently from each other, a convenient parallelization of the wavefront propagations among many processors is possible, requiring the use of HPC in most cases.
2.5. Coherent mode decomposition of undulator radiation
The crossspectral density, generally represented as
expresses the correlation of the emitted radiation between any two spatial points r_{1} and r_{2}. It is the fundamental object that we will use to describe all partially coherent properties of the synchrotron beams. We justify first, in the context of the existing literature, the conditions of its usage for synchrotron light. Then we present the CMD and its practical implementation in 2D (with COMSYL) and the proposed 1D algorithm (with WOFRY1D).
2.5.1. Validity of CSD usage for emission in storage rings
Geloni et al. (2008) show^{3} that, although synchrotron radiation emission (a random process obeying Gaussian statistics) is not intrinsically stationary nor homogeneous, the secondorder coherence theory of scalar fields as presented by Mandel & Wolf (1995) can be applied when the following conditions are observed:
(1) Radiation frequency ω is `sufficiently high'.
(2) ebunch timelength σ_{T} is `sufficiently large' so that ωσ_{T} ≫ 1.
(3) Radiation bandwidth Δ_{ω} is not `too narrow' (Δ_{ω} ≫ 1/σ_{T}).
Excluding infrared frequencies and below, condition (1) holds for soft and hard Xrays; condition (2) is satisfied by storage rings, where typically σ_{T} > 30 ps, but not at freeelectron lasers, where σ_{T} < 0.1 ps due to microbunching effects; and, finally, condition (3) is generally met by standard monochromatization schemes. This set of conditions, related to the longitudinal electronbeam direction, ensures that synchrotron radiation emission is a quasistationary (or a widesense stationary) process.
In the (2D) transverse direction, following Geloni et al. (2008), the source is quasihomogeneous when (i) the intensity slowly varies compared with the coherence length, and (ii) the complex degree of coherence depends only on Δr = r_{2} − r_{1}. The following additional conditions delimit the applicability of quasihomogeneity:
(4) N_{x} ≫ 1 and/or D_{x} ≫ 1.
(5) N_{y} ≫ 1 and/or D_{y} ≫ 1, with
where σ_{x,y} and represent the electron beam transverse sizes and divergences, L_{u} is the undulator length (number of periods N_{u} times the magnetic period λ_{u}) and = λ/2π. If conditions (4) and (5) are met with `and', we have a Gaussian quasihomogeneous source. For thirdgeneration synchrotron light sources, condition (4) is met with `and' and condition (5) with `or'. Here the source is quasihomogeneous and can be separated in H and V, as shown by Geloni et al. (2008) [equation (56)]. As the horizontal emittance reduces, as it is the case of fourthgeneration light sources (for the cases studied here N_{x} = 13, D_{x} = 2.6, N_{y} = 0.4, D_{y} = 0.3), we slowly approach a region where quasihomogeneity breaks down.
2.5.2. Coherent modes, coherent fraction and coherent length
The CSD in equation (7) is used to define the spectral density^{4} as
for the case where r = r_{1} = r_{2}. We also define the normalized crossspectral density function or spectral degree of coherence^{5}, hereafter referred to as DoC, as
The modulus value of equation (10) is limited to 0 ≤ 1, where μ = 0 means total uncorrelation and μ = 1 denotes full correlation of the fluctuations at positions r_{1} and r_{2}.
A well known result from coherence theory is the coherent mode representation of partially coherent fields in free space [see §4.7.1 in Mandel & Wolf (1995)]. It is possible to decompose the CSD into an infinite sum of orthonormal coherent modes,
where Λ_{n} (eigenvalues) are the intensity weights and Φ_{n} are the coherent modes (eigenfunctions). Some important characteristics of this coherent mode decomposition are:
(1) The modes Φ_{n} are orthonormal in the integral sense.
(2) The decomposition maximizes the CSD making the truncation optimal,
(3) The eigenvalues Λ_{n} are a measure of the intensity of the corresponding mode Φ_{n}.
(4) We define the occupation η of the ith mode as its normalized intensity,
(5) Radiation is considered fully coherent if and only if there is only a single mode.
From these arguments, it is now natural to rigorously define coherent fraction (CF) as the occupation of the first coherent mode,
The transverse coherence length (CL) is related to the width of a cut of the modulus of the DoC. There is no unanimous accepted way of defining the CL, a parameter of practical importance for daily beamline operations. Quantitative values of CL are discussed in §5.1. Importantly, the blind application of the van Cittert–Zernike theorem may lead to errors, as discussed in §4.1 of Geloni et al. (2008), because this theorem was originally derived for incoherent sources.
The interest in the coherent mode decomposition method applied to the optical design of Xray beamlines is manyfold: (i) the possibility of propagating a partially coherent beam along a beamline by just propagating coherent modes; (ii) the ability to compute the CSD and therefore all the related coherence parameters from the propagated modes; (iii) the use of the coherent fraction (a scalar parameter) as a well defined measure of how coherent is the beam at any position of the beamline; (iv) the numerical storage of the N_{m} modes that depend on two spatial variables is more economic than the storage of the CSD, a complex function of many variables; and (v) the infinite series converges smoothly to the CSD, therefore the truncation at a limited number of modes N_{m} always guarantees that it is the best possible reduced representation of the CSD and can be quantified.
2.5.3. Coherent mode decomposition with COMSYL
COMSYL (COherent Modes for SYnchrotron Light) is a software package to perform the CMD of undulator radiation in a storage ring (Glass, 2017a). A complete description of the code is given by Glass (2017b) and here we summarize the principles underlying it. Applications of this software for beamline modelling at the ESRFEBS are presented by Glass & Sanchez del Rio (2017) and Sanchez del Rio et al. (2019). COMSYL was used to study the specked structure of the CSD and the presence of Xray coherence vortices and domain walls (Paganin & Sanchez del Rio, 2019).
Coherent mode decomposition consists of calculating Λ_{i} and Φ_{i} in equation (11). This operation can be seen as the `diagonalization' of W where the eigenfunctions are the coherent modes (Φ_{i}) and the eigenvalues their intensity weights (Λ_{i}). These are the solution of the homogeneous Fredholm integral equation of the second kind,
The eigenvalue in equation (14) can be written
where we define the Hermitian operator for a generic function g as
A first look at the CSD expression [see equation (5)] is enough to get an impression of how it is resourceintensive calculating and storing this 6D function. For synchrotron beams it is useful to decouple the longitudinal coordinate – along the optical axis in a beamline (see Fig. 1) – so that the CSD reduces to 4D, and wavefronts are 2D. We use W_{2D} notation for this CSD. Knowing the CSD at a given position z, that is, , it is possible to propagate it to another position z′ and also backpropagate to the source position.
Kim (1986) developed a propagation theory of synchrotron radiation using Wigner distribution. He introduced the `brightness convolution theorem' stating that the source due to a collection of electrons randomly distributed in their phase space is calculated by a convolution of the source due to a reference electron with the electron distribution function. COMSYL applies Kim's convolution theorem in a plane (s_{x}, s_{y}) at s_{0} = −L_{u}/2, which is where the electrons enter the undulator. This distance s_{0} is taken from the centre of the ID (origin of the optical axis) – see Fig. 1. We have [cf. equation (3.37) in Glass (2017b)]
where r_{e} = (x_{e}, y_{e}), θ_{e} = , ) and Δr = r_{2} − r_{1}. The electric field E_{ω} is calculated using SRW in an arbitrary plane at a distance z and then backpropagated to s_{0} using the standard Fresnel freespace propagator (see Section S2 of the supporting information).
COMSYL starts with a matrix method that discretizes the crossspectral density operator in a step function basis set [see equation (4.4) of Glass (2017b)]. The discretization is followed by an iterative diagonalization using SLEPc (Hernandez et al., 2005). The implementation uses parallel computation and requires the use of an HPC. The key point of COMSYL is to avoid storing the full representing matrix because it requires significant memory and computational resources. It scales essentially with N_{x}^{2}×N_{y}^{2} where N_{x} and N_{y} are the numbers of grid points in the x and y dimension, respectively. Typical sizes for N_{x} and N_{y} can easily reach a few hundred up to a few thousand. In the latter case, the memory requirements would reach several thousand terabytes. To reduce the memory requirements of the matrix method, COMSYL uses a twostep method that first performs a coherent mode decomposition for a zero divergence electron beam and based on this decomposition performs a second decomposition that takes the divergence into account. The memory requirement for our undulator applications is drastically reduced to about 4 × N_{x} × N_{y} × N_{m} where N_{m} is the number of requested coherent modes. This allows the calculation of higher harmonics or higher emittance rings where N_{x} × N_{y} ≫ N_{m}.
The modes calculated with the justdescribed method can be propagated along the beamline and used to compute the spectral density and crossspectral density at any point of the beamline. These modes can be conveniently propagated downstream of a beamline with SRW or WOFRY using the OASYS environment (see §3). However, due to modifications by optical elements (cropping/truncation and/or absorption), the propagated wavefronts generally lose their orthonormality. Thus, for computing the CF at a given point of the beamline, it is necessary to perform a new CMD with the propagated CSD. It is possible to apply the very same method used to compute the initial coherent modes, but COMSYL implements the action of the integral operator [equation (14)] directly to the coherent modes, which is much more economic in terms of memory usage.
2.5.4. Coherent mode decomposition with separate horizontal and vertical directions (WOFRY1D)
Resuming the discussion in §2.5.1, we now assume to be in a quasihomogeneous regime. This allows us to decompose the CSD as a product of horizontal and vertical crossspectral densities,
The CSD for the horizontal (x) direction is
and similarly for the vertical direction (replacing x by y), with λ_{n} the eigenvalues (scalars) and ϕ_{n} the 1D eigenfunctions. The CSD described in equation (19) is now a 2D function. This dimension reduction is very welcome as it becomes very easy to store, manipulate and diagonalize the CSD using common tools (e.g. the Python numpy and SciPy libraries). Much like that presented in §2.5.3 for the COMSYL software, we calculate the 1D undulator emission at an arbitrary plane located at z from the origin using pySRU (Thery et al., 2016), and backpropagate this field to the centre of the undulator using WOFRY's zoom propagator (see Section S3 of the online supporting information). W_{1D} is then obtained using Kim's convolution theorem at z = 0, which is effectively constructed using equation (3.51) of Glass (2017b).
This recipe is implemented in the WOFRY wavefront propagation addon in OASYS. For a typical beamline, two calculations are done: one for the horizontal direction and another for the vertical. These two results can be combined to obtain the CSD = × , implying numerical tensor operations. Similarly, the 2D spectral density is retrieved as = . This will be extensively used in later sections. Note that if the intensity is stored in arrays, the product is indeed an outer product. Finally, this decomposition gives rise to two values of coherence fraction: CF_{x} for the horizontal direction and CF_{y} for the vertical direction. The 2D CF can be retrieved as CF = CF_{x} × CF_{y}. Like in COMSYL, the ϕ_{n} eigenfunctions can be propagated to any position of the beamline like standard wavefronts. The propagated W_{1D} can, again, be decomposed in coherent modes to obtain the CF after propagation.
The factorization in equation (18) has been discussed by Geloni et al. (2008) [see their §3.1, equation (56)] where it is stated that this rearranging of equation (7) is only valid for quasihomogeneous sources (see §2.5.1). In the context of the Wigner distribution, this separation in horizontal and vertical components has been discussed by various authors and is believed to work well for undulator radiation with photon energies near the resonance of odd harmonics (Bazarov, 2012; Tanaka, 2014; Nash et al., 2021). The Cartesian factorization in equation (18) is supposed to work well when the electron beam parameters in the horizontal are different from the vertical ones. The Cartesian separation does not work for rotationally symmetric sources (round beams). However, the Wigner function in this case (and therefore the CSD, being related by a Fourier transform) can be treated as a 1D problem, as suggested by Agarwal & Simon (2000) and developed by Gasbarro & Bazarov (2014). It would be interesting for future fourthgeneration sources that will create round beams to develop a CMD tool using polar coordinates, also including wave propagators like those discussed by Li & Jacobsen (2015).
2.6. Hybrid raytracing (HYBRID)
Simulations methods using raytracing are simpler and faster than those using waveoptics, and they offer useful insight when studying and designing a synchrotron beamline (Sanchez del Rio et al., 2019). However, pure raytracing methods (based on geometrical optics) are not adequate when analyzing optical systems dealing with (partially or fully) coherent beams subjected to strong diffraction effects (e.g. beam clipping by either physical acceptance of an optical element or by slits and optical errors). It is possible to add diffraction effects to the geometric methods, by convolution of the beam divergence profiles calculated by raytracing with those resulting from the diffraction integrals, then proceeding with classical raytracing methods. This hybrid concept (Shi et al., 2014) is implemented as an extension to the raytracing code SHADOW (Sanchez del Rio et al., 2011) available in the ShadowOUI (Rebuffi & Sanchez del Rio, 2016) addon in OASYS. Since its release, this HYBRID raytracing implementation has been successfully used as an efficient and fast tool to design beamlines also including coherence effects (Shi et al., 2017; Rebuffi & Shi, 2020; Lordano, 2022).
3. Propagation of wavefronts along the beamline
In terms of design using physical optics, an Xray beamline is composed of two main elements: drift spaces and optical elements. The first category is handled with diffraction integrals, and a brief explanation of the wavefront propagators used by the software in Section 2 is presented here. Optical elements can usually be represented by transmission elements and their treatment is also covered in this section.
3.1. Drift spaces
Under the scalar theory, a generic wavefront obeying the wave equation and completely described at a position z, that is, E_{ω}(r), known for all the xy plane, will propagate (evolve) between two parallel planes separated by a distance L = z′ − z as
Equation (20) is the first Rayleigh–Sommerfeld diffraction equation (Huygens–Fresnel principle) and is valid for the case where r′ − r ≫ λ, with λ = 2π/k. We define a normal vector parallel to the optical axis (zaxis) ℓ so that θ is the angle between −ℓ and the vector r′ − r; Σ is the xy plane in z where the integration takes place with = (see Fig. 1). Further simplification to equation (20) can be done using the paraxial approximation. In this case, it is assumed that ≃ 1; and that the term r′ − r = [(x′ − x)^{2} + (y′ − y)^{2} + L^{2}]^{1/2} can be expanded in a Taylor series with L^{2} ≫ (x′ − x)^{2} and L^{2} ≫ (y′ − y)^{2}. Retaining the quadratic term in the exponential function, but dropping it for the denominator,
This approximation is known as the Fresnel diffraction integral and strategies for its numerical calculation are plenty (Kelly, 2014; Goodman, 2017); see also Rees (1987), Stern & Javidi (2004) and Zhang et al. (2020) for other practical considerations.
The different waveoptics packages in use have different implementations of these propagators. The selection of the propagator and its parametrization (sampling, padding, interpolation, etc.) is made depending on the particular characteristics of the optical element and propagation distances. SRW uses four different propagators, summarized in Section S2 of the supporting information. WOFRY1D basically uses two different propagators (see Section S3).
3.2. Optical elements
Optical elements will interact with the wavefront by manipulating its amplitude and phase. A very wide range of optical elements can be represented by the complex transmission element,
which is applied to the electric field E_{ω}(r) as a multiplicative factor (Cloetens et al., 1996). This thinelement approximation can be expanded to represent thick optical elements by applying the multislicing method, which represents an optical element as intercalation of several thin elements (slices) and freespace propagation between them (Paganin, 2006; Li et al., 2017; Munro, 2019) (see Section S4 of the supporting information).
A generic aperture (slit, pinhole and beamstop) is represented by a binary opaque object: it masks part of the wavefront,
When the element is a slit, then T = 1. If it is a beamstop, then T = 0. A is the region defined by the object boundary.
Absorption filters, test patterns, Xray refractive lenses, refractive correctors (freeform refractive optics), zone plates and transmission gratings are usually well represented by this thinelement approximation in projection approximation with
Δ_{z} is the projected thickness of an optical element along the zaxis. We define the index of refraction as n_{ω} = 1 − δ_{ω} + iβ_{ω}. Optical errors can be simulated in a similar way (Laundy et al., 2014; Celestre et al., 2020; Sanchez del Rio et al., 2020).
Refractive systems, like the lenses used in this work, are calculated using the thinobject approximation with the lens described by its profile and refraction index (material). Usually, a single lens has a parabolic interface z = x^{2}/(4R) with R the radius at the apex. A lens has usually two parabolic interfaces (front and back) separated by a lens thickness d_{L}. The interfaces are flat outside the lens aperture A. Therefore, the lens thickness profile is
used to compute the complex transmission with equation (22).
Reflective optics can also be reasonably well approximated by transmission elements with more complex calculations, e.g. stationaryphase methods [application of equation (20) over the mirror surface], geometric raytracing calculation of the optical path difference (Canestrari et al., 2014) or even multislicing (Li et al., 2017). Alternatively, 1D wavefronts can be easily propagated in grazingincidence mirrors and gratings calculating directly the numeric integral in equation (20) (Raimondi & Spiga, 2015; Sanchez del Rio et al., 2020).
4. Description of the optical system
The optical system under consideration is based on the future `EBSL1' beamline at ESRF.^{6} This will be a long beamline (200 m) for applications exploiting the beam coherence. The source considered is an undulator with 138 periods of 18 mm (length close to 2.5 m). We analyze a focusing system with two transfocators, at 65 and 170 m from the source. They contain sets of 2D and 1D lenses that will permit modifying independently the focal lengths for the horizontal and vertical directions. The use of two transfocators allows great flexibility in beam transport (Vaughan et al., 2011). The first one can be used to modify the divergence of the beam, even to collimate it, to guarantee a full illumination at the second transfocator.
Each of the two transfocators in use (T1 and T2 in Fig. 2) are idealized as two crossed 1D Be lenses. For each plane (H and V), lens1 and lens2 have variable curvature radii R_{1} and R_{2} that match the focal distances f_{1} and f_{2} (f = R/2δ), with δ = 6.96 × 10^{−6} for Be at 7 keV. The focal lengths for the lenses are different for the horizontal and vertical directions to adapt to the beam characteristics. A slit (CS in Fig. 2) of aperture a_{x} in the horizontal and a_{y} in the vertical is placed upstream of the lens1. We set the distances matching the requirements of the EBSL1 beamline (see Fig. 2), and we analyzed the system at a photon energy of 7 keV for different focal distances of lens1 and lens2.
We are interested in the beam properties (intensity distribution, size, flux) at the sample plane for four cases. The first case is selected to obtain a small spot (about 5 µm) and the second one a large spot (more than 30 µm). For these cases, the slits are selected to match CF_{x} = CF_{y} = 90% for a photon energy of 7 keV. The values are shown in Table 1. Cases 3 and 4 follow the same logic but the slits are opened to increase intensity at the expense of reducing coherence (CF_{x} = CF_{y} = 70%).

5. Results and discussion of multioptics simulations
Calculations are done using the four different methodologies discussed previously, implemented in four different addons of the OASYS ecosystem.
5.1. Source characteristics and its propagation to the entrance slit
We first calculated the source and the illumination at the entrance slit plane (z = 36 m). The spectral density calculated with the different codes is shown in Fig. 3. The distributions calculated by the different methods are close with differences in full width at halfmaximum (FWHM) of less than 10%.
The CSD can be calculated using wave optics codes. Figure 4 shows the horizontal and vertical W_{1D} at the source plane calculated with SRWME and WOFRY1D. Again, an excellent agreement is found, with similar FWHM of the profiles crossing (0, 0) [9 µm for both codes in the horizontal (H) and 12 µm in the vertical (V)].
Figure 5 shows the horizontal and vertical DoC at the slit plane expressed in the new coordinates () = [(x_{1} + x_{2})/2, x_{2} − x_{1}] (for the horizontal, similarly with y for the vertical). The interest in using these coordinates is to redress the plot of the CSD that lies on a diagonal (as shown in Fig. 4) and obtain the `coherent length' (CL) as a `width' of the modulus of the DoC versus Δ_{x} at = 0. If using the FWHM, we obtain a horizontal CL of 76 µm with WOFRY1D and 80 µm with SRWME and a vertical CL of 444 µm with WOFRY1D and 402 µm with SRWME. A naive application of the van Cittert–Zernike theorem for a source with a Gaussian intensity profile permits the calculation of the coherence length (the FWHM of the Fourier transform of the source intensity profile) as CL = 0.88λz/S, with z the source–observation plane distance, and S the source FWHM. In our case [z = 36 m, source FWHM 70.6 µm (H) and 15.0 µm (V), and λ = 1.77 × 10^{−10} m] it gives CL = 79 µm (H) and CL = 374 µm (V). The (rough) agreement of the values from the van Cittert–Zernike theorem (corresponding to a fully incoherent source) with our numerical values (for the partially coherent beam) is justified by the fact that z is large enough to lie in the zrange where the CL is linear [see discussion and Fig. 17 of Geloni et al. (2008)].
It is worth mentioning that the modulus of the DoC (and therefore the CL) is obtained experimentally by measuring the interference of the two beams originated by a doubleslit (Thompson & Wolf, 1957). Examples of this type of experiment with synchrotron radiation are presented elsewhere (Chang et al., 2000; Paterson et al., 2001; Leitenberger et al., 2003; Tran et al., 2005). We performed simulations with WOFRY1D placing two slits of 2.5 µm with horizontal separation s_{A} in the plane at z = 36 m, and propagating the resulting two beams to z = 46 m. The results are shown in Fig. 6(a), where it can be appreciated that the visibility of the fringes decreases when increasing s_{A}. For a given s_{A}, the intensity profile [e.g. Fig. 6(b)] is used to compute the visibility = which is equal to the modulus of the DoC. We then obtained the visibility versus the slit separation s_{A} that gives (as expected) the same values as the modulus of DoC versus x_{2} − x_{1} in Fig. 5 [see Fig. 6(c)].
The principal role of the slit is to control the beam coherent fraction. Closing the slit will increase the CF, with an obvious decrease in integrated intensity. The slit aperture necessary to obtain a `good' coherence is somehow related to the CL, but quantitative values are better calculated using the CF versus the slit aperture. Within the CMD theory, this requires calculating coherent modes after the slit, which can be easily done with WOFRY1D (see Fig. 7). From the CF versus slit aperture plot, one can pick the aperture values to match the desired CF (we selected 90% or 70% values of CF to select slit apertures in Table 1), and at the same time estimate the losses in flux.
5.2. Image at sample position
The intensity distribution at the sample plane is displayed in Figs. 8(a)–8(d), for results using the WOFRY1D, COMSYL, SRWME and HYBRID codes. Beam dimensions are obtained by calculating the FWHM from the intensity distribution in one direction, resulting from integration along the other direction. They are displayed in the plots and summarized in Table 2. The results for case 1 show a horizontal profile mostly triangular with shoulders that evidence small diffraction fringes. The fringes are more resolved in the vertical direction. Case 2 presents in the horizontal a soft Gaussianlike profile, but in the vertical important symmetric shoulders are visible. Case 3 shows a smooth Gaussian profile in the horizontal and a small shoulder with fringes in the vertical. Case 4 shows a conventional smooth profile in the horizontal but an original threelobe plateau in the vertical. This variety of profile distribution demonstrates how relevant the diffraction effects are, which modulate the beam shape in a nontrivial way.

The WOFRY1D results are shown in Fig. 8(a). The 1D intensity profile for each direction is obtained from the summation of several modes. High modes have very low eigenvalues. It is sufficient to consider only ten modes for accounting for more than 99% of the spectral density. The 2D intensity distribution shown in Fig. 8(a) is obtained by combining the calculated horizontal and vertical 1D profiles via the outer product. We can observe in the intensity distributions the same structures due to the diffraction effects as those observed with the other calculation methods. The beam profiles calculated, obtained with COMSYL and propagated with WOFRY, are shown in Fig. 8(b). SRWME results are given in Fig. 8(c). The good convergence of the values displayed is guaranteed by a convergence analysis described in Section S5. Hybrid raytracing for the four cases defined in Table 1 are shown in Fig. 8(d). The obtained intensity distributions are less structured than those calculated with the other waveoptics methods (e.g. the threelobe plateau in case 4v is not reproduced). However, the FWHM values agree with full waveoptics calculations for most cases [except for two particular cases: 1h (HYBRID 17.3 µm, WOFRY1D 8.5 µm) and 2v (HYBRID 74.0 µm, WOFRY1D 32.4 µm)]. They will be discussed in the next section.
The agreement between the results of WOFRY1D in Fig. 8(a) and SRWME [Fig. 8(c)] is striking. All intensity distributions reproduce the same features, and the differences in FWHM values are less than 12%, a value that is compatible with the errors of the simulations. This result validates the 1D CMD method proposed here, whose requirements in computer power are extremely low (it runs very fast on an average laptop).
The numeric value for sizes calculated with the different methods (Table 2) depends not only on the code itself but also on the particular specific parameters in each method (number of pixels for sampling wavefronts, propagation parameters, etc.). To estimate the calculation error in the final size numbers, we vary randomly these specific parameters over a reasonable range (e.g. 10%). The dispersion (standard deviation) of the sizes obtained is a good estimator of the error in this parameter. This exercise would take considerable computational effort using 2D methods, but it can be easily done with WOFRY1D. We run 200 cases with 10% random variation in the number of pixels and zoom factor for drift spaces. The obtained sizes (horizontal × vertical) are 8.49 ± 0.60 µm × 4.97 ± 0.37 µm (case 1), 39.94 ± 2.98 µm × 32.77 ± 2.69 µm (case 2), 36.39 ± 2.89 µm × 6.12 ± 0.51 µm (case 3), and 24.18 ± 1.80 µm × 133.44 ± 10.06 µm (case 4). We confirmed that the values given in Fig. 8(a) are within these error intervals.
The calculated beam sizes should be completed with ^{15} photons s^{−1} (0.1% bandwidth)^{−1}. Each of the three elements studied (slit, lens1 and lens2) absorbs part of the The estimation of the absorption by the slit can be done using simple geometrical arguments, and the absorption by the lenses depends on the average Be thickness presented to the beam. The of Be at 7 keV is μ = 3 cm^{−1}, giving 1.45% attenuation for a 50 µmthick layer (like the lens thickness used in the simulations^{7}). From the simulated data we extracted the absorption for the different absorbing elements (slit, lens1 and lens2) (see Table 3). We note a high absorption in lens2 in cases 1 and 3; this is because lens2 is overilluminated, therefore the 1 mm physical aperture absorbs the beam considerably.
At 7 keV, the undulator in the selected configuration emits a of 1.5 × 10

5.3. Computer resources
COMSYL requires highperformance computing to perform full CMD of the undulator beam, by solving the Friedholm problem and obtaining the full 2D eigenfunctions (coherent modes) and eigenvalues. The simulation of the source with COMSYL used to calculate Fig. 8(b) took 55 min using 28 × 3.30 GHz CPUs of 251.82 GB RAM, for getting 174 modes of 1691 × 563 pixels. The modes calculated by COMSYL were propagated with WOFRY in the OASYS environment (Rebuffi & Sanchez del Rio, 2017). The propagation used the 2D zoom propagator (see Section S3) and the optical elements described in Section 3.
The good convergence of the SRWME results in Fig. 8(c) is guaranteed by a convergence analysis described in Section S5 of the supporting information. It was used to determine the minimum number of electrons that produce accurate results. The SRWME simulations for the cases analyzed converge with only a few thousand electrons in a node with 28 CPUs totalizing 256 GB. This is because the beam after the slit has a relatively high CF.
Concerning running times, the simulations with the new WOFRY1D code run in a few seconds on a laptop. HYBRID also requires light computer resources and also runs on a laptop. The simulation of the full CMD with COMSYL required about 1 h with 28 cores. The source is then reused for propagating the different configurations. SRWME required a full source simulation for each configuration that also runs in about 1 h with 5000 electrons in a similar cluster.
5.4. Further simulations
The simulations presented, motivated by the EBSL1 project, use a simplified optical layout. They are used to validate the WOFRY1D tool before proceeding with further analyses.
A systematic and exhaustive study has been carried out using WOFRY1D for the new EBSL1 beamline also including other optical elements not considered here and multiple transfocator configurations. It is important to match correctly the two transfocators to guarantee that the focal point is located at the sample plane. Heat load deformation must be controlled at the whitebeam mirrors and also at the monochromator. For that, the deformations calculated by finiteelement methods are used in WOFRY1D for assessing the optical impact, as we did previously (Sanchez del Rio et al., 2020). It is imperative to study the effect of mirror slope errors and surface errors at the lenses [as described by Celestre et al. (2020)]. These results will be described elsewhere.
6. Summary and conclusions
We presented in Section 2 the theory of partially coherence optics applied to undulator radiation and its implementation in two different solutions: (i) Monte Carlo multielectron simulations as implemented in SRWME and (ii) coherent mode decomposition as implemented in COMSYL and in the new package WOFRY1D. The key point of WOFRY1D is that it needs very scarce computer resources. The factorization of the CSD in two directions (horizontal and vertical) is possible in many cases of major interest when the undulator is tuned close to oddharmonic resonances, and when the horizontal and vertical emittances of the storage ring are not the same (i.e. nonround beams, as for the lowemittance storage ring EBSESRF). Section 3 summarizes the propagation of wavefronts along empty drift spaces and thin objects, which include the elements used in our simulations: slits and Xray lenses.
We studied a particular case of focusing a partially coherent beam produced by an undulator in EBSESRF by a system of two transfocators (implemented as single parabolic lenses). This case is of interest to the new EBSL1 beamline being constructed at the ESRF. The combined effect of beam diffraction at the slit and global focusing by the two lenses produces images with a variety of intensity profiles (see Fig. 8). We included HYBRID raytracing results as this method can be used in the first simulation phase and produces approximated values of beam sizes and flux.
We have verified that simulations with the new WOFRY1D code are consistent with the other three simulation codes typically used to simulate synchrotron beamlines. A remarkable agreement is found between WOFRY1D and SRWME for the functions that describe partial coherence (CDS, DoC and CL) (see Figs. 4 and 5). We further used WOFRY1D to discuss the coherent fraction versus slit aperture (Fig. 7) and to verify that the modulus of the DoC is consistent with the results of a (simulated) experiment of twobeam interference (Fig. 6).
Partial coherence calculations using 2D wavefronts are expensive from the computation point of view, either because many thousands of wavefronts are propagated (like in SRWME) or because of the need to diagonalize an extremely large 4D crossspectral function (like in COMSYL). The newly developed coherent mode decomposition uses 1D wavefronts and is very rapid and light. It can run interactively on a laptop. This opens new paths for intensive simulations of experiments using partially coherent beams, in particular for imaging applications or beamline optimization, where thousands of runs are necessary. It can also serve as a simulation engine for systems exploiting machine learning or for digital twins of beamline instruments.
The software tools developed here are available in the WOFRY addon of the OASYS suite (Rebuffi & Sanchez del Rio, 2017). The OASYS workspaces and scripts for the simulations performed in this work are also available from https://github.com/oasysesrfkit/papermultiopticsresources.
7. Related literature
The following references, not cited in the main body of the paper, have been cited in the supporting information: Chubar (2021); Chubar & Celestre (2019); Cowley & Moodie (1957); He et al. (2020); Pirro (2017); Schmidt (2010).
Supporting information
Extension and detailed information of some sections. DOI: https://doi.org/10.1107/S1600577522008736/ay5600sup1.pdf
Footnotes
^{1}At thirdgeneration light sources, very restrictive pinholes and slits were used for spatial filtering with a dramatic loss of when improving and tuning the coherent fraction of the beam. Due to the increased coherent fraction at fourthgeneration sources, much less restrictive filtering is necessary, with higher photon transmission as a consequence.
^{2}Geloni et al. (2008) provide a counterexample discussing the highresolution monochromator of their reference 15.
^{3}During the review of this manuscript, the authors came across the work of Trebushinin et al. (2022), which also covers the theoretical background of undulator statistical properties.
^{4}Spectral density is sometimes loosely called intensity. The spectral density and the intensity functions are equivalent for the stationary case. The same holds for the CSD and the mutual optical intensity (MOI) (Mandel & Wolf, 1995).
^{5}It is sometimes called complex spatial (or spectral) degree of coherence, see §4.3.2 in Mandel & Wolf (1995).
^{6}Throughout this work we have used the following values of electron beam sizes at the centre of the straight section: σ_{x} = 29.7 µm, = 4.37 µrad, σ_{y} = 5.29 µm, = 1.89 µrad, corresponding to beam emittances: ε_{x} = 130 pm rad, ε_{y} = 10 pm rad, and beta functions β_{x} = 6.8 m, β_{y} = 2.8 m.
^{7}In the simulations, the horizontal and vertical focusing are separated in two lenses, with accumulated thickness 100 µm thus absorption 3%.
Acknowledgements
G. Geloni (EU XFEL) is acknowledged for discussion on statistical optics applied to undulator radiation in lowemittance storage rings; O. Chubar (NSLSII/BNL) for clarifications on selected wavefront propagators; S. Lordano for sharing information on beamline design for the Sirius Light Source; and J. P. Guigay for discussions about DoC and CL.
Funding information
This project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreements No. 823852 (Photon and Neutron Open Science Cloud – PaNOSC) and No. 101007417 (NFFAEurope Pilot Joint Activities – NEP).
References
Agarwal, G. S. & Simon, R. (2000). Opt. Lett. 25, 1379–1381. Web of Science CrossRef PubMed CAS Google Scholar
Bazarov, I. V. (2012). Phys. Rev. ST Accel. Beams, 15, 050703. Web of Science CrossRef Google Scholar
Canestrari, N., Chubar, O. & Reininger, R. (2014). J. Synchrotron Rad. 21, 1110–1121. Web of Science CrossRef IUCr Journals Google Scholar
Celestre, R., Berujon, S., Roth, T., Sanchez del Rio, M. & Barrett, R. (2020). J. Synchrotron Rad. 27, 305–318. Web of Science CrossRef IUCr Journals Google Scholar
Chang, C., Naulleau, P., Anderson, E. & Attwood, D. (2000). Opt. Commun. 182, 25–34. Web of Science CrossRef CAS Google Scholar
Chubar, O. (2021). SRW, https://www.github.com/ochubar/SRW. Google Scholar
Chubar, O., Berman, L., Chu, Y. S., Fluerasu, A., Hulbert, S., Idir, M., Kaznatcheev, K., Shapiro, D., Shen, Q. & Baltser, J. (2011). Proc. SPIE, 8141, 814107. CrossRef Google Scholar
Chubar, O. & Celestre, R. (2019). Opt. Express, 27, 28750–28759. Web of Science CrossRef PubMed Google Scholar
Chubar, O. & Elleaume, P. (1998). Proceedings of the Sixth European Particle Accelerator Conference (EPAC98), 22–26 June 1998, Stockholm, Sweden, pp. 1177–1179. THP01G. Google Scholar
Cloetens, P., Barrett, R., Baruchel, J., Guigay, J.P. & Schlenker, M. (1996). J. Phys. D Appl. Phys. 29, 133–146. CrossRef CAS Web of Science Google Scholar
Cowley, J. M. & Moodie, A. F. (1957). Acta Cryst. 10, 609–619. CrossRef IUCr Journals Web of Science Google Scholar
Dejus, R. J. & Luccio, A. (1994). Nucl. Instrum. Methods Phys. Res. A, 347, 61–66. CrossRef CAS Web of Science Google Scholar
Gasbarro, A. & Bazarov, I. (2014). J. Synchrotron Rad. 21, 289–299. CrossRef IUCr Journals Google Scholar
Geloni, G., Saldin, E., Schneidmiller, E. & Yurkov, M. (2008). Nucl. Instrum. Methods Phys. Res. A, 588, 463–493. Web of Science CrossRef CAS Google Scholar
Glass, M. (2017a). COMSYL GitHub repository, https://www.github.com/oasyskit/comsyl. Google Scholar
Glass, M. (2017b). Statistical optics for synchrotron emission: numerical calculation of coherent modes. Phd thesis, Université Grenoble Alpes, France (https://tel.archivesouvertes.fr/tel01664052). Google Scholar
Glass, M. & Sanchez del Rio, M. (2017). Europhys. Lett. 119, 34004. Web of Science CrossRef Google Scholar
Goodman, J. W. (2017). Introduction to Fourier Optics, 4th ed. W. H. Freeman and Company. Google Scholar
He, A., Chubar, O., Rakitin, M., Samoylova, L., FortmannGrote, C., Yakubov, S. & Buzmakov, A. (2020). Proc. SPIE, 11493, 114930H. Google Scholar
Hernandez, V., Roman, J. E. & Vidal, V. (2005). ACM Trans. Math. Softw. 31, 351–362. CrossRef Google Scholar
Hirschmugl, C. J., Sagurton, M. & Williams, G. P. (1991). Phys. Rev. A, 44, 1316–1320. CrossRef CAS PubMed Web of Science Google Scholar
Jackson, J. D. (1999). Classical Electrodynamics. New York: Wiley. Google Scholar
Kelly, D. P. (2014). J. Opt. Soc. Am. A, 31, 755. CrossRef Google Scholar
Kim, K.J. (1986). Proc. SPIE, 0582, 2–9. CrossRef CAS Google Scholar
Laundy, D., Alcock, S. G., Alianelli, L., Sutter, J. P., Sawhney, K. J. S. & Chubar, O. (2014). Proc. SPIE, 9209, 920903. Google Scholar
Leitenberger, W., Wendrock, H., Bischoff, L., Panzner, T., Pietsch, U., Grenzer, J. & Pucher, A. (2003). Physica B, 336, 63–67. Web of Science CrossRef CAS Google Scholar
Li, K. & Jacobsen, C. (2015). J. Opt. Soc. Am. A, 32, 2074–2081. CrossRef Google Scholar
Li, K., Wojcik, M. & Jacobsen, C. (2017). Opt. Express, 25, 1831–1846. Web of Science CrossRef PubMed Google Scholar
Li, R. & Chubar, O. (2022). Opt. Express, 30, 5896–5915. CrossRef PubMed Google Scholar
Lordano, S. (2022). Private communication. Google Scholar
Mandel, L. & Wolf, E. (1995). Optical Coherence and Quantum Optics. Cambridge University Press. Google Scholar
Munro, P. R. T. (2019). J. Opt. Soc. Am. A, 36, 1197–1208. CrossRef Google Scholar
Nash, B., Goldring, N., Edelen, J., Webb, S. & Celestre, R. (2021). Phys. Rev. Accel. Beams, 24, 010702. CrossRef Google Scholar
Paganin, D. M. (2006). Coherent Xray Optics. Oxford University Press. Google Scholar
Paganin, D. M. & Sanchez del Río, M. (2019). Phys. Rev. A, 100, 043813. CrossRef Google Scholar
Paterson, D., Allman, B. E., McMahon, P. J., Lin, J., Moldovan, N., Nugent, K. A., McNulty, I., Chantler, C. T., Retsch, C. C., Irving, T. H. K. & Mancini, D. C. (2001). Opt. Commun. 195, 79–84. Web of Science CrossRef CAS Google Scholar
Pirro, G. (2017). Application of Scaled Wave Optics Propagator to Model Synchrotron Beamlines. MSc thesis, Politecnico di Milano, Italy. Google Scholar
Raimondi, L. & Spiga, D. (2015). Astron. Astrophys. 573, A22. Web of Science CrossRef Google Scholar
Rebuffi, L. & Sánchez del Río, M. (2016). J. Synchrotron Rad. 23, 1357–1367. Web of Science CrossRef IUCr Journals Google Scholar
Rebuffi, L. & Sanchez del Rio, M. (2017). Proc. SPIE, 10388, 103880S. Google Scholar
Rebuffi, L. & Shi, X. (2020). Proc. SPIE, 11493, 1149303. Google Scholar
Rees, W. G. (1987). Eur. J. Phys. 8, 44–48. CrossRef CAS Google Scholar
Sanchez del Rio, M., Canestrari, N., Jiang, F. & Cerrina, F. (2011). J. Synchrotron Rad. 18, 708–716. Web of Science CrossRef CAS IUCr Journals Google Scholar
Sanchez del Rio, M., Celestre, R., Glass, M., Pirro, G., Herrera, J. R., Barrett, R., da Silva, J. C., Cloetens, P., Shi, X. & Rebuffi, L. (2019). J. Synchrotron Rad. 26, 1887–1901. Web of Science CrossRef IUCr Journals Google Scholar
Sanchez del Rio, M., Wojdyla, A., Goldberg, K. A., Cutler, G. D., Cocco, D. & Padmore, H. A. (2020). J. Synchrotron Rad. 27, 1141–1152. Web of Science CrossRef IUCr Journals Google Scholar
Schmidt, J. D. (2010). Numerical Simulation of Optical Wave Propagation. SPIE Press. Google Scholar
Shi, X., Reininger, R., Harder, R. & Haeffner, D. (2017). Proc. SPIE, 10388, 10388. Google Scholar
Shi, X., Reininger, R., Sanchez del Rio, M. & Assoufid, L. (2014). J. Synchrotron Rad. 21, 669–678. Web of Science CrossRef IUCr Journals Google Scholar
Stern, A. & Javidi, B. (2004). Opt. Eng. 43, 239–250. CrossRef Google Scholar
Tanaka, T. (2014). Phys. Rev. ST Accel. Beams, 17, 060702. Web of Science CrossRef Google Scholar
Tanaka, T. & Kitamura, H. (2001). J. Synchrotron Rad. 8, 1221–1228. Web of Science CrossRef CAS IUCr Journals Google Scholar
Thery, S., Glass, M. & Sanchez del Rio, M. (2016). PySRU GitHub repository, https://www.github.com/oasyskit/pySRU, https://github.com/oasyskit/pySRU. Google Scholar
Thompson, B. J. & Wolf, E. (1957). J. Opt. Soc. Am. 47, 895–902. CrossRef Google Scholar
Tran, C. Q., Peele, A. G., Roberts, A., Nugent, K. A., Paterson, D. & McNulty, I. (2005). Opt. Lett. 30, 204–206. CrossRef PubMed CAS Google Scholar
Trebushinin, A., Geloni, G., Rakshun, Y. & Serkez, S. (2022). Optica, 9, 842–852. CrossRef Google Scholar
Vaughan, G. B. M., Wright, J. P., Bytchkov, A., Rossat, M., Gleyzolle, H., Snigireva, I. & Snigirev, A. (2011). J. Synchrotron Rad. 18, 125–133. Web of Science CrossRef IUCr Journals Google Scholar
Walker, R. P. & Diviacco, B. (1992). Rev. Sci. Instrum. 63, 392–395. CrossRef Web of Science Google Scholar
Wiedemann, H. (2019). Particle Accelerator Physics, 4th ed. Springer International Publishing. Google Scholar
Xu, H., Zhu, Z., Li, X., Liu, P., Dong, Y. & Zhou, L. (2022). Opt. Express, 30, 7625–7635. CrossRef PubMed Google Scholar
Zhang, W., Zhang, H. & Jin, G. (2020). Opt. Express, 28, 39916–39932. CrossRef PubMed Google Scholar
This is an openaccess article distributed under the terms of the Creative Commons Attribution (CCBY) Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the original authors and source are cited.