Modelling undulators in ray tracing simulations

We introduce a model that can accurately simulate radiation from undulator sources for ray tracing applications. It incorporates several key effects relevant to fourth-generation synchrotron sources, such as electron emittance, energy spread, and diffraction-limited beam size. This code has been developed as part of SHADOW4, the latest version of the widely used SHADOW X-ray optics ray tracing program. The approach relies on calculating the field distribution in the far field, which determines the ray divergences. The integration of existing models for electron energy spread is also addressed. Rays sampled at the source follow a size distribution derived by backpropagating the far-field radiation. These models are detailed, and several examples are provided.

1. Introduction

Undulators are the most popular magnetic structures for producing synchrotron radiation in third- and fourth-generation sources. The radiated beam by an undulator is usually more brilliant than the other sources: it is much more collimated than in wigglers and bending magnets, and it is as intense or more than the wiggler at certain photon energies. The theory of undulator radiation (UR) is well understood and several comprehensive texts are available (Duke, 2000; Onuki & Elleaume, 2003; Clarke, 2004). The radiation emitted by the undulator exhibits distinct structures both in its spectrum, presenting peaks at some photon energies (resonances), and in its geometry (wavefront size and derived divergences). Several software tools are available to compute the characteristics of the UR. Among them, SRW (Chubar & Elleaume, 1998) and SPECTRA (Tanaka & Kitamura, 2001) are the most advanced.

Ray tracing packages create undulator sources by sampling rays according to the distributions given by the undulator theory. SHADOW, since its first version SHADOW1 (Cerrina, 1984), included an undulator model (Chapman et al., 1989). In SHADOW3 (Sanchez del Rio et al., 2011) and its ShadowOui interface (Rebuffi & Sanchez del Rio, 2016), the undulator calculations were refactored, and partially replaced by new Python code. Moreover, ShadowOui also provides an `Undulator Gaussian' application, that creates a source with rays that follow a Gaussian distribution that approximate the undulator distributions. This has been found very useful when in a first phase or prototyping beamlines using undulators as sources. In SHADOW4 (Sanchez del Rio & Rebuffi, 2023), the newly refactored and enhanced version of this popular ray tracing code, we have reimplemented, improved, corrected, and upgraded the undulator algorithms, significantly improving the performance and accuracy of ray tracing simulations. Special attention has been given to accurately implement certain features that have become crucial with the emergence of fourth-generation synchrotron sources. Notably, this includes the impact of electron energy spread, which is particularly important when utilizing radiation at high harmonics, as well as the accurate description of diffraction-limited beam size.

This paper aims to describe the methods and algorithms used in SHADOW4 for simulating undulator sources. In a first section we summarize the most important results of the undulator theory used, with a detailed discussion of the Gaussian approximations for beam sizes and divergences. Next, the algorithms and methods for sampling rays from undulator sources, both in the Gaussian approximation and in the full model, are described. Finally, several examples are provided, followed by a discussion.

2. Summary of the theory of undulator emission

The undulator magnets induce in the electrons a periodic (mostly sinusoidal) trajectory. The small deflection of each electron at each oscillation of the magnetic field makes it possible that the photons produced in a crest of the electron periodic trajectory interfere with the photons originating from the next oscillation crest, thus producing radiation with non-smooth characteristics (spectral and spatial). The spectrum contains peaks at photon energies proportional to the so-called resonance. It depends on the deflection parameter K. The K value for an electron traveling in an oscillating magnetic field (with z the spatial coordinate along the undulator, B the maximum magnetic field, and λ_u the undulator period) is

with m and e the mass and charge, respectively, of the electron, and c the velocity of the light. The resonance is found at the photon wavelength

with θ the observation angle (θ = 0 on-axis) and γ ≃ 1957 [GeV] the Lorentz factor with the electron energy.

2.1. Emission from a relativistic electron

An ultrarelativistic charged particle traveling along a curved, often wiggly, trajectory, typically generated by alternating magnetic fields in insertion devices (IDs), emits radiation. The electric field can be calculated in the framework of classical electrodynamics [see, for example, equation (14.14) of Jackson (1999)]. The electric field at an observation point r = (x, y, z) is proportional to the following integral along the electron trajectory,

where c is the velocity of light, ω is the radiated frequency, β = is the electron relative velocity, and the dot denotes the time derivative. Also n(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 center of the ID/straight section.

Equation (3) describes a fully spatially coherent field emitted by a single electron. In an idealized zero-emittance storage ring, the electrons follow a `filament beam'. It contains N_e electrons that follow exactly the same trajectory r_e(t), therefore the radiation intensity, calculated from the square of equation (3), will be affected by a multiplicative factor N_e, or for practical effects will be expressed as a function of the electron current.

Several codes are available in the synchrotron community to calculate the undulator emission characteristics in different cases. The codes URGENT (Walker & Diviacco, 1992) and UR or US (Dejus & Luccio, 1994) compute undulator emission in the far field for undulators with a sinusoidal magnetic field. The codes SRW (Chubar & Elleaume, 1998) and SPECTRA (Tanaka & Kitamura, 2001) are more generic as they calculate emission in the near and far field for any electron trajectory (with different initial conditions) and submitted to an arbitrary magnetic field. We use pySRU (Thery et al., 2016), an open source code developed in Python, that calculates the emission using equation (3). It is well integrated in Python ecosystems, such as OASYS (Rebuffi & Sanchez del Rio, 2017), which includes the SHADOW4 user interface. Portions of pySRU have been incorporated into the internal code of SHADOW4.

The flux spectrum F(E), with F the flux in photons s⁻¹ (0.1% bandwidth)⁻¹ and E = ℏω the photon energy (in eV), is obtained by fixing a coordinate z (the distance from the center of the undulator to an observation plane) and integrating over the x, y variables (x is in the horizontal plane and y in the vertical plane¹). The spectrum `on-axis' [i.e. integrated over an infinitesimal interval of (x, y)] presents peaks at energies corresponding to the values in equation (2) (with θ = 0). These peaks have the form of a function with x = (Elleaume, 2003) (N is the number of undulator periods, E₀ is the resonance energy, and n is the harmonic number). As far as one opens the integration range in θ or in (x, y) (or in other words, we open an acceptance slit) the peaks become wider because the resonance shifts with the conservation angle θ in equation (2).

The intensity map of the radiation for z sufficiently large (far field) does not change its shape but it only expands with z. We can then speak about `divergence', in terms of the radial angle θ = (x² + y²)^1/2/z (the horizontal angle is θ_x = x/z and the vertical one θ_y = y/z). Under some approximations (far field, sinusoidal field, on-resonance) the intensity reduces to [equations (23)–(24) of Elleaume (2003)]

with = θ(L/λ)^1/2, L is the undulator length and λ the photon wavelength at a given harmonic.

To obtain the source size (intensity map at z = 0), equation (3) cannot be directly applied to points within the electron trajectory. Therefore, it is necessary to compute the electric field at a z point external to the undulator and then backpropagate it (using Fresnel or Fraunhofer propagators) to the plane at z = 0. Using, as before, some approximations, the source size can be expressed as a Hankel transform (the Fraunhofer propagator in radial coordinates) of equation (4) which gives [equation (29) of Elleaume (2003)]

where J₀(x) is the Bessel function of the first kind and zeroth order, and is the `reduced' radial coordinate, = 2πr(2λL)^−1/2.

2.2. Gaussian approximation of undulator size and divergence at resonance

The divergence or angular distribution of the UR can be calculated by representing the flux (F) as a function of the horizontal (θ_x) and vertical (θ_y) angles, or the radial angle = ]. Near resonance (and its odd harmonics), the distribution displays a pronounced peak, known as the `central cone', along with some surrounding rings. At exact resonance, the distribution is described by equation (4). As far as the photon energy is reduced (red-shifted), the peak opens a hole in the middle that separates into two or more peaks to eventually disappear.

The width of the intensity profile of the radiation cone is a fundamental parameter for researchers and engineers working with synchrotrons. Different works found in the literature use different approximations under multiple hypotheses with some discrepancy in the results. We summarize them from a historic perspective and identify the equations implemented in SHADOW4.

The natural divergence of synchrotron light for all sources (bending magnets, wigglers and undulators) is approximately proportional to γ⁻¹. For undulator sources, it is smaller by a factor that depends on the number of undulator periods. Using simple arguments (Krinsky, 1983) affirms that the angular broadening of the radiation is defined as

The same expression is used by Kim (1986a; Kim, 1986b) and is also given in the X-ray Data Booklet (Thompson, 2001) or by equation (14.21) of Duke (2000). Notice that in the original texts the width is said to be a `half width', which in principle is different from a `sigma' (σ) in a Gaussian distribution, which has a full width at half-maximum of ≃ 2.355σ. In Kim's papers from 1989 (Kim, 1989), the Gaussian approximation is obtained by matching its integral with the angular distribution of intensity [the sinc function in equation (4)] and obtaining a smaller divergence,

Elleaume (2003) performed a Gaussian fit on the intensity versus emission angle at the resonance [equation (4)] and obtained

Note that in this equation the numeric factor 0.69 is very close to = 0.707 in equation (7). Therefore, if not identical, they are in close agreement (within 2%). We repeated this fit with the same result [see Fig. 1(a)]. However, one can remark by a simple visual inspection that the fit is not good: the intensity profile is far from being Gaussian. Another practical method to obtain the σ value is to compute the root mean square (r.m.s.) of the intensity distribution, identical to the standard deviation because the mean is zero. This is practical for numeric calculations, but may end in infinite values using the theoretical equations, as discussed by Elleaume (2003). Therefore, care must be taken when using these approximated values. In recent publications, most authors agree with equation (7), for example Tanaka (2014) and Walker (2019).

Figure 1 Gaussian fits of the intensity versus reduced emission angle (a) and reduced size (b) at the center of the undulator plane, after Onuki & Elleaume (2003).

The situation is more controversial when discussing the radiation size σ_r, or the spatial width of the radiation at the center of the undulator. Before summarizing the bibliographic results, we remind that a Gaussian beam (the first mode of a Gaussian–Shell model beam) verifies

Kim [equation (21) of Kim (1986b) or equation (6.37) of Kim (1989)] supposed that the UR verifies equation (9). Therefore, they implicitly assumed the validity of approximating UR by a Gaussian Schell-model, or, in other words, accepting that the emission at the resonance is Gaussian. In consequence, depending on the divergence value used, two results are found: combining equations (9) and (6), Kim obtains [equation (28) of Kim (1986b)]

or, combining equations (9) and (7) [equation (6.37) of Kim (1989)],

Elleaume (Onuki & Elleaume, 2003) followed a different direction. He did not suppose that the observed approximately Gaussian divergence comes from the Fraunhofer propagation of a Gaussian beam as hypothesized by Kim, but obtains a numerical fit of the calculated radiation expressed as a function of the real space at the source position [equation (5)] and obtained

The fitted Gaussian for spatial [equation (12)] and angular [equation (8)] representations of the UR are not related via Fourier transform [equation (9)], or, in other words, their product is not λ/(4π) but

which can be interpreted as the phase space volume of UR being approximately twice the phase space volume of the first coherent mode of a Gaussian beam.

In the literature, we can find papers that follow Elleaume's [e.g. Borland (2012); Hettel (2014)] and Kim's [e.g. Huang (2013)] model. To conclude this section, it is worth mentioning, citing Onuki & Elleaume (2003), that `these are approximations and should not be considered as fundamental results'. Moreover, Walker (2019) affirms that `the reason why different models for the source size and divergence have been put forward is that the radiation phase space is not at all Gaussian in nature'. Recently, several papers discussed the undulator's phase space and obtain expressions of brighness, coherence, etc. adapted to new generations of sources (Geloni et al., 2008; Tanaka, 2014; Lindberg & Kim, 2015; Walker, 2019).

In the SHADOW4 code, and for the following discussion, we adopt Elleaume's approach [equations (8) and (12)].

2.3. Description of electron sizes and emittance

At any position s in the storage ring, an electron can be described by five coordinates: = 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 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), thus defining 2 × 2 matrices. For x (horizontal plane),

and similarly for the y coordinate (vertical plane). We also indicate the expression as a function of the Twiss functions (α, β and γ) and emittance [ε_x = (〈xx〉〈x′x′〉 − 2〈xx′〉²)^1/2]. In some particular points or the storage ring, the covariance between spatial and angle terms is zero (ρ = α = 0), thus only the diagonal terms ( are sufficient to define the electron beam. This is the case at the center of the straight sections, where the undulators are usually placed. When the undulator is in another position at a distance s from the center of the straight section (with Twiss parameters α₀, β₀ and γ₀), the new parameters are

2.4. Divergence and size of the photon source at resonance in Gaussian approximation

Consider a filament beam that emits a radiation wavefront. At resonance, the size and divergence distributions are supposed Gaussians given by equations (12) and (8), respectively. Consider now that the emission is not given by a filament beam, but instead by a bunch of electrons distributed, as discussed, with values , , , .

Suppose first that all the electrons have exactly the same energy ( = 0). We can assume that the emission of the photon source is the convolution of the photon source (filament beam) with the electron beam. Therefore, the sizes and divergences of the photon source are

Tanaka & Kitamura (2009) have studied the effect of the electron energy dispersion in the UR. They found that the approximated photon source size and divergences are

where corrective terms for sizes Q_s and angles Q_a have been introduced. They depend on the electron `normalized energy spread' σ_ε = , with the electron energy dispersion, N the number of periods of the undulator and n the harmonic number in use. The Q_a functions is

with erf(x) the error function (Wikipedia contributors, 2024); and the Q_s function² is

Geloni et al. (2018) treat how the brightness is influenced by the energy spread in a context not restricted to Gaussian approximations. They obtain the approximated expressions compatible with those of Tanaka & Kitamura (2009) and argue that they `may constitute a good approximation in some region of the parameter space, when it comes to the limit for a diffraction-limited beam with non-negligible energy spread, a more detailed study is needed'.

For the calculations in this paper we use the parameters of the ESRF U18 undulator installed in the ID06 beamline at EBS-ESRF. It has a period of λ_u = 18 mm, N = 111, and therefore a length of 2 m. The K value ranges from 0.2 to 1.479. The resonance is set to E = 10 keV with K = 1.3411. The EBS storage ring has electron beam energy = 6 GeV, electron energy spread = 9.334 × 10⁻⁴ and electron sizes at the center of the straight section where the undulator is installed σ_x = 30.18 µm, σ_y = 3.64 µm, = 4.37 µrad, = 1.37 µrad. They give beam emittances: ɛ_x = = 132 pm rad, and ɛ_y = = 5 pm rad.

As an example, the undulator sizes and divergences for the ESRF ID06 U18 undulator calculated using equation (18) are shown by the solid lines of Fig. 2. From Figs. 2(c) and 2(d) it is evident that for high harmonics the divergence is influenced by the electron energy spread, thus the necessity of including it in the ray tracing simulations. The effect in beam size is moderate [Figs. 2(a) and 2(b)]. Walker (2019) discussed the sizes and divergences (ibid., Figs. 8 and 9) in the context of numeric calculations of brightness also including the energy spread. The results for divergence r.m.s. agree well with the model used here. For the size r.m.s. they observe a discrepancy mainly due to changes in the shape (narrowing the cone and introducing wide tails) of the distribution that becomes less and less Gaussian when increasing .

Figure 2 Full width at half-maximum values of size [in the horizontal (a) and vertical (b) directions] and divergence [in the horizontal (c) and vertical (d) directions] of the photon source as a function of the resonance photon energy for the ID06-U18 at EBS/ESRF. The continuous lines are calculated using equation (18). The points (dots) correspond to ray tracing using the `Undulator Gaussian' widget. Every point corresponds to the average value after 50 SHADOW4 simulations, using a random Monte Carlo seed each. The corresponding standard deviation is shown as error bars [note the small range of horizontal size on the left (a) resulting in large error bars].

We verified numerically the suitability of Q_a to correct the angular width of the undulator emission for the previously discussed ESRF ID06 U18 undulator tuned at 10 keV (first harmonic). We supposed here zero emittance and a variable energy spread from zero to 0.005. We calculated numerically, using the WOFRY wave optics package (Sanchez del Rio et al., 2024), the intensity distribution at 100 m (far field). Without changing the undulator configuration, we repeated the calculation for different values of electron energy around = 6 GeV. We then constructed the pattern for each energy spread by summing the patterns for each electron energy weighted by a Gaussian with the corresponding . We finally calculated the FWHM and the r.m.s. values of each intensity pattern and normalized them to the value obtained for = 0. Fig. 3 shows the results calculated for the far field [Figs. 3(a) and 3(b)] and backpropagated to the center of the ID [Figs. 3(c) and 3(d)]. In Fig. 3(b) a comparison of Q_a with the numerical values at the far field of the FWHM and SD (a width calculated from the r.m.s. as if it was Gaussian, i.e. SD = 2.355 × r.m.s.) is shown. We observe a good agreement of Q_a with the numerical values (both FWHM or SD) for n = 1. The agreement is less good for higher n and . For the backpropagated radiation we see, as previously noticed, that the effect of the electron energy spread is moderate. Indeed, it seems from Fig. 3(c) that there is a shrink in the width when increasing . When examining the numerical values in Fig. 3(d) we observe a discrepancy between the FWHM and SD, indicating a non-Gaussian behavior. While SD increases (as predicted by Q_s), the observed FWHM slightly decreases. This is attributed to a narrowing of the peak accompanied by an expansion of the tails, as noted by Walker (2019). In summary, the positive takeaway is that, in all cases, Q_a and Q_s values fall between the numerical values of FWHM and SD, highlighting the difficulty of selecting a single parameter to describe a non-Gaussian distribution. Keeping in mind that the value of is close to 0.001 for most synchrotron sources, we conclude that the use of Q_a and Q_s is a reasonable choice for incorporating the electron energy dispersion in ray tracing simulations when working at resonance.

Figure 3 Results of the wave optics simulation of the intensity versus electron energy spread calculated for the ID06 U18 undulator. The radiation is calculated at the far field [(a) and (b)] and backpropagated to the center of the undulator [(c) and (d)]. (a) Intensity map at the far field as a function of calculated at the first harmonic (10 keV). (b) Values of the angular-width correction (normalized to the width at zero energy dispersion) as a function of the electron energy spread. The widths are measured by the FWHM and the SD value (2.355 × r.m.s.) of the far-field intensity profiles [as shown in (a) for n = 1]. They are compared with Qa. (c) The same as (a) but for the backpropagated radiation. (d) The same as (b) but for the size corrections compared with Qs.

2.5. Divergence and size of the photon source off-resonance

The values of beam size and divergence obtained in the previous section are valid only when working with photons at resonance (or at a particular odd harmonic). It is common that the experimentalist set the monochromator close to, but not exactly at, resonance. For example, the photon energy corresponding to the maximum intensity integrated over a finite θ interval (e.g. using a slit) is not exactly at resonance, but red-shifted by an amount that depends on the aperture. Moreover, going out of resonance, the intensity distribution changes from a well defined peak (Fig. 1) to other shapes, also presenting a double-peak. This is illustrated in Fig. 4 where numeric values of FWHM and SD are computed for photon energies scanning the first harmonic peak. It can be observed [see Fig. 4(b)] that the minimum of the divergence is obtained at a position blue-shifted with respect to the resonance, but [see Fig. 4(d)] the minimum FWHM of the size tends to a red-shifted position.

Figure 4 Results of wave optics simulations of the intensity versus photon energy calculated for the ID06 U18 undulator. The radiation is calculated at the far field [(a) and (b)] and backpropagated to the center of the undulator [(c) and (d)]. (a) Intensity map at the far field (100 m) as a function of photon energy calculated for the first harmonic (E0 = 10 keV). (b) Values of the width of the intensity distribution at the far field as a function of the photon energy. The widths are obtained from the FWHM and the SD value (2.355 × r.m.s.). The intensity (in arbitrary units) is also shown (green curve). (c) The same as (a) for the backpropagated radiation at the center of the undulator. (d) The same as (b) for the backpropagated radiation at the center of the undulator.

The effect of detuning of the electron energy has a similar effect as detuning the photon energy from resonance. Indeed, from equation (2), = (λ_u/2)(1 + K²/2) when θ = 0 (on-axis). Therefore is a constant for a particular tuned undulator, or, in other words, a wavelength (or energy) shift is equivalent to a corresponding shift of the electron energy. One can compensate the other as far as is constant. This is illustrated in Figs. 5(a) and 5(b).

Figure 5 Maps of the intensity at the far field (at 100 m from the source, where x is the vertical spatial coordinate). (a) Intensity versus photon energy shift from resonance (E0 = 10 keV) and x ( = 6 GeV). (b) Intensity versus electron energy (at resonance E0).

Nash et al. (2019) present functions to correct the flux, size and divergence considering the electron energy spread and the detuning from the resonance. They built two-dimensional maps of functions that correct the flux, size and divergence standard deviation, versus both the energy spread and the spectral detuning. In this way, the mentioned correction functions Q_a,s(σ_ε) are replaced by functions F_a,s(E − E₀, σ_ε) also including the detuning from resonance E − E₀. Following the same idea, we calculated numerically using WOFRY the maps of the r.m.s. [Fig. 6(a)] and the FWHM [Fig. 6(b)] of the far-field intensity for the ESRF ID06 undulator around the resonance. Because we want to observe the changes versus , we normalize each value to the corresponding value at = 0. The differences observed between these two maps indicate, again, that the distributions are not Gaussian, therefore FWHM is not related to r.m.s. by the constant 2.355. Another observation is that the most changes (values that separate from one) are observed for small values of E − E₀. We compared these results with the values from Nash et al. (2019) applied to our particular undulator. The map for angles [Fig. 6(c)] agrees well with our corresponding r.m.s. map [Fig. 6(a)]. For completeness we also evaluated the size correction F_s(E − E₀, σ_ε) [see Fig. 6(d)]. Our results, confirmed by those of Nash et al. (2019), conclude that the corrections by energy spread are important at the resonance, but are not so important far from resonance. This remark will be used in the undulator model presented in Section 3.2.

Figure 6 (a) Calculation of the r.m.s. value of the intensity distribution at 100 m downstream of the ID06 U18 undulator. (b) The same as (a) but the map represents the FWHM. (c) The correction function Fa calculated using data from Nash et al. (2019) applied to the same undulator. (d) The correction functions Fs calculated using data from Nash et al. (2019) applied to the same undulator. Data in each map are normalized the values at zero electron energy spread.

3. Algorithms used in ray tracing undulator sources

We present in this section the ideas behind the models for creating an undulator source with SHADOW4, in its two different applications: the simplified and approximated `Undulator Gaussian', and the more accurate `Undulator Light Source'.

3.1. The `Undulator Gaussian'

It is generally a good idea to start ray tracing simulations for prototyping an undulator beamline with a simplified and quick model. This is done with the `Undulator Gaussian'. It supposes that we work at the resonance energy, or at an energy very close to it to suppose the results at resonance are applicable. The rays are sampled following Gaussian distributions in both size and divergence. We use equation (17) to calculate the sigmas of the photon source. For this, the only requirement is to know the electron beam sizes, the working photon energy, and the undulator length. In SHADOW4 we give the option to consider the effect of the electron energy spread. If we activate this option, the user is requested to enter the values of and also N and n (the undulator period and the harmonic in use) to correct the sizes and divergences using equation (18).

Due to the assumptions made, the source is considered monochromatic. However, when modeling crystal monochromators, which typically have a very narrow energy bandwidth, it is beneficial to generate a polychromatic source over an energy range slightly broader than the monochromator's acceptance. To achieve this, we introduce an energy interval ΔE around the resonance, within which ray energies are sampled according to a simplified flat distribution.

A key enhancement in SHADOW4 is that the sources now include information on the number of photons, allowing the application to directly provide data on absorbed and transmitted intensity and power. This eliminates the need to manually rescale the number of rays to represent the number of photons. In the simplified model of the `Undulator Gaussian', the user can either enter this value manually or allow it to be calculated using the equation of the flux in the central cone [equation (17) of Kim (2001)],

where α is the fine-structure constant, I_SR is the electron beam current, Δω/ω is the photon energy bandwith, typically 10⁻³; and Q_n(K) = F_n(1 + K²/2)(1/n), with F_n a universal function defined by Kim (2001).

3.2. The `Undulator Light Source'

This is the primary application to simulate undulator sources. The user selects input parameters for the electron beam (sizes, divergences, energy and current), undulator parameters (N, period and K) and sampling parameters (number of rays N_rays and the photon energy interval). An important parameter is θ_max, the maximum radial angle to be considered in the calculations. It affects the total number of photons and the sampling of the rays. Other `advanced' parameters permit the number of sampling points to be defined, and flags to activate some modeling options described below.

In SHADOW4, like in the original SHADOW1, the sampling of points that follow a given (1D, 2D or 3D) distribution is done using the `inverse method', an optimized algorithm proposed by John von Neumann in a famous letter to Stan Ulam [Fig. 3 of Eckhardt (1987)]. The application of this method to sample undulator ray energies and divergences is discussed in detail in Section 6 of Chapman et al. (1989).

The steps to create the SHADOW4 source are the following:

(i) Construct a 3D stack of the electric fields (for σ and π polarizations) of the radiation emitted by the undulator in the far field as a function of the photon energy E, the radial angle θ and the azimuthal angle φ for the filament beam. This requires a previous calculation of the electron trajectory. The 3D stack has N_E × N_θ × N_φ points, which is set by the user. In the case that the user wants a `monochromatic' source, N<sub>E</sub> = 1. This main step uses the theory of undulators based on equation (3). Like in Chapman et al. (1989), the stack is computed in polar coordinates, that is more efficient than using Cartesian coordinates, because the `almost' axial symmetry permits to limit N_φ to low values.

(ii) From these stacks of electric field for the two polarizations A_E,σ and A_E,π, compute the stacks of intensity = |A_E,σ|² + |A_E,π|² and polarization = |A_E,σ|/(|A_E,σ| + |A_E,π|).

(iii) When the electron energy spread is taken into account, the array containing the angular array θ is multiplied by Tanaka's function Q_a calculated at the resonance. This option is only allowed when simulating monochromatic sources. As discussed before, the user should only activate this correction when working at the resonance or very close to it.

(iv) is used as a 3D probability density function, thus integrate it to obtain the 2D and 1D cumulative distribution functions. Then obtain the sampled arrays of the photon energy and divergences (directions) of the rays. Up to here, the rays are directed as if there were emitted by a filament beam. Then, correct these directions for electron beam emittance by adding sampled values that follow Gaussian distributions with and . This is equivalent to performing the mathematical convolution. In case of finite Twiss α, use a 2D Gaussian as defined in equation (15).

(v) Calculate the polarization for each ray by interpolating with the sampled E_i,θ_i,φ_i values. With it, construct the SHADOW electric field vectors A_σ and A_π. The first is directed along the horizontal axis and the second along the vertical axis. The intensity for each ray is normalized to unity: |A_σ|² + |A_π|² = 1.

(vi) Sample the ray positions. For the undulator, the photon source has `no depth', i.e. z = 0 for all rays. Contrary to other sources (bending magnets and wigglers) where the rays start from different positions along the electron path, the undulator rays start from a volume corresponding to the backpropagation of the far field to the plane z = 0. The horizontal and vertical distribution are due to (i) the backpropagated far emission or the filament beam (a kind of `diffraction limit' size), and (ii) the electron sizes σ_x and σ_y. They are calculated and combined as follows:

 (a) The backpropagated field for the filament field can be selected among three options: (i) to neglect it, setting a point source x = y = z = 0 (before convolution with electron beam sizes), as done in the first model of SHADOW1 (a good approximation for `high-emittance' storage rings like those of the 20th century because the electron sizes are much larger than the `diffraction-limited' sizes); (ii) sample rays following the Gaussian approximation for the emission at resonance [equations (12)] (this is a good approximation when simulating monochromatic sources at the resonance); and (iii) use a most accurate method backpropagating the far-field radiation to the plane z = 0 and sampling rays accordingly. This is a costly and delicate operation as it involves a careful sampling of the radiation field, and implementing and setting the wave propagation for the N_E photon energies. A compromise has been found to obtain a reasonable solution without exploding the calculation time (see Appendix A).

 (b) If the option of considering the electron energy spread is on, the array with the sizes r is multiplied by the Tanaka's correcting function Q_s. Again, this option is only available for monochromatic sources.

 (c) Once the size distribution is found, the rays are sampled and corrected by adding the Gaussian sampling of the electron source with σ_x and , or from a 2D Gaussian [equation (15)] in the case of finite Twiss α.

(vii) The collection of rays with sampled photon energies, directions (divergences), sizes, and electric field (polarization) constitute the ray tracing source.

4. Examples and discussion

Here we present calculations using the two SHADOW4 undulator sources, the `Undulator Gaussian' and the `Undulator Light Source' that implement the methods described in Sections 3.1 and 3.2, respectively. The aim is to confirm that the rays generated at the source to represent the undulator accurately match the expected intensity distributions based on theoretical predictions. This serves as a benchmarking process to validate the reliability of the new code.

We use for the calculation the ESRF U18 undulator described before. For testing the `Undulator Gaussian', several simulations have been performed in monochromatic mode at different photon energies. From the intensity distributions as a function of spatial or angular coordinates the FWHM has been evaluated. The results are shown in Fig. 2, overplotted with the calculations using the analytical expressions of the energy spread effect [equation (18)]. The plot includes error bars in the SHADOW simulations, derived from the standard deviation of multiple runs. It is observed that the theoretical value falls within the error bars of the ray tracing results, confirming that the parameters extracted from the rays do align with the underlying theoretical model.

Another interesting feature of SHADOW4 undulators is the possibility of getting the source photon flux and use it for simulations. In Fig. 7 we compare the flux calculated with the `Undulator Gaussian' [that using equation (21)] with the ones obtained from OASYS add-on XOPPY that uses SRW.

Figure 7 Flux comparison from harmonic 1 to 9 of SHADOW4 Gaussian undulator and XOPPY/SRW.

Using the full undulator application `Undulator Light Source', we first calculate the intensity distribution as a function of the vertical angle at three different photon energies (exactly on resonance, slightly red-shifted and a bit blue-shifted). The results are given in Fig. 8, showing that the three methods implemented in SHADOW4 produce similar results (the default `internal' method, the use of the WOFRY library, and the use of SRW to compute the radiation in the far field).

Figure 8 Divergences for the ID06 U18 undulator, set close to the first harmonic resonance. (a) 2D plot for photon energy red-shifted E = E0[1 − (Nn)^−1] = 9910 eV, (b) 2D plot at resonance E = E0 = 10000 eV, (c) 2D plot at blue-shifted E = E0[1 + 0.4(Nn)^−1] = 10036 eV, all calculated with the `internal' algorithm. On the right, we compare the histograms of the vertical divergences given by the three calculation modes SHADOW4 internal algorithm, pySRU + WOFRY, and SRW for the three photon energies: (d) red-shifted from resonance, (e) at resonance, and (f) blue-shifted from resonance.

In Fig. 9 we study the effect of the electron energy spread  = 0.001 on the vertical divergence at the n = 5 harmonic. The same intensity profiles are obtained using the three modes of calculation in SHADOW4, namely the `internal' method, the one using `pySRU + WOFRY', and the one using SRW.

Figure 9 Left column: vertical divergences for the ID06 U18 undulator at 50 keV (n = 5, in green) and two off-resonance energies (red- and blue-shifted). Right column: including energy spread ( = 0.001). (a) SHADOW4 internal algorithm, (b) pySRU + WOFRY, and (c) SRW. Electron beam emittance has been considered.

It is often useful to simulate a polychromatic source, for example, when covering the energy range of a given undulator harmonic. In Fig. 10, we display the U18 first and third harmonic histograms using the undulator in polychromatic mode (with the resonance energy at 10 keV). Both with 101 energy points, the first harmonic covers a range from 9.6 keV to 10.2 keV in an emission cone of 16.6 µrad, and the third harmonic from 29.4 keV to 30.4 keV in 9.6 µrad. We compared the intensity distribution from the rays with the theoretical one.

Figure 10 Energy histogram from the rays of the U18 first (a) and third (b) harmonic using the undulator in polychromatic mode. Both are compared with the theoretical flux distribution that has been used in the sampling process.

To verify the proper functioning of the polychromatic option, we compared the rays generated by it with those produced by multiple monochromatic sources, scanned over a range of photon energies. This was done in 101 steps, ranging from 9.6 keV to 10.2 keV, with each energy step intensity weighted by its flux. We then compared these results with the size and divergence obtained from the polychromatic mode, the same simulation from which the first harmonic in Fig. 10 was obtained. In Fig. 11, we show the comparison results of both planes sizes and divergences profiles. We note that the correction for electron energy spread is not available when using the polychromatic option. This limitation arises from the challenge of extending the correction functions Q_a,s, derived for resonance energy (monochromatic case), to the polychromatic scenario. It is recommended to begin simulations with the simplest configuration (monochromatic setup) to assess whether electron energy spread corrections are necessary. If these corrections are significant, the polychromatic option should be avoided, necessitating an iteration over monochromatic sources. However, as noted earlier, these corrections are typically negligible outside resonance, reducing the impact of their unavailability in the polychromatic case.

Figure 11 Comparing the polychromatic source with the accumulation of monochromatic steps. (a) Horizontal size, (b) vertical size, (c) horizontal divergence and (d) vertical divergence. Intensities are normalized to the intensity peak.

Finally, we evaluated the diffraction-limited size. Fig. 12 displays the patterns at three photon energies near resonance, similar to Fig. 8. Slight differences are evident among the various backpropagation methods implemented in SHADOW4. Fine-tuning the backpropagation parameters for the selected method is necessary due to their significant impact on the source size.

Figure 12 Sizes at the center of the ID for the ID06 U18 undulator, set close to the first harmonic resonance. (a) 2D plot for photon energy red-shifted E = E0[1 − (Nn)^−1] = 9910 eV, (b) 2D plot at resonance E = E0 = 10000 eV, (c) 2D plot at blue-shifted E = E0[1 + 0.4(Nn)^−1] = 10036 eV, all calculated with the `internal' algorithm. On the right we compare the histograms of the vertical size given by the three calculation modes SHADOW4 internal algorithm, pySRU + WOFRY, and SRW for the three photon energies: (d) red-shifted from resonance, (e) at resonance, and (f) blue-shifted from resonance.

5. Summary and conclusions

This work details the development and improvements of undulator sources in SHADOW4, a ray tracing code designed for synchrotron beamline modeling. Key enhancements address critical elements for the fourth-generation synchrotron sources, including electron energy spread and diffraction-limited beam size.

A description and analysis of the existing models for the effects of the electron energy spread demonstrate the need to correct radiation divergence when using high harmonics. This is effectively modeled by the Q_a function from Tanaka & Kitamura (2009), as confirmed by our wave-optics numerical simulations. The size correction is considerably less significant, primarily affecting the shape of the intensity distribution rather than its overall width. Additionally, we found that the corrections in angle and size are not significant for photon energies outside of resonance. These insights are valuable when utilizing the new undulator features in SHADOW4: `Undulator Gaussian' and `Undulator Light Source.'

The software tools developed here are part of the SHADOW4 add-on available in the OASYS suite (Rebuffi & Sanchez del Rio, 2017). The OASYS workspaces and scripts for the simulations performed in this work are also available (Sanchez del Rio & Reyes-Herrera, 2024).

APPENDIX A

SHADOW4 methods for backpropagating the far-field radiation

SHADOW4 stores the radiation field in a 3D stack in polar coordinates of dimension (N_E, N_θ, N_φ). It is calculated at a large distance or far field (typically D_far = 100 m). The radiation produced by a filament electron beam constitutes a single (coherent) wavefront. We want to calculate the intensity map of this wavefront at the source position. Thus we need to backpropagate the wavefront from D = D_far to D = 0. Backpropagation requires the use of a propagator. This propagator can be implemented using the full numeric integral (which is impractical because of the long calculation time) or using two Fourier transforms (for the Fresnel propagator). The latter requires a Cartesian gridding of the wavefront with equally distributed points. We discarded the solution of interpolating the electric field from the polar grid (N_θ, N_φ) to a Cartesian dense grid (N_x, N_y). The interpolation introduces high errors in the wavefront (particularly in the phase) that strongly affect the propagated wavefront. We left for a future work the implementation of a good Fresnel propagator in polar coordinates, using the Hankel transform (instead of Fourier transform in Cartesian coordinates).

Currently three solutions are available in SHADOW4: an internal solution based on backpropagation of a 1D wavefront, and other solutions based on external libraries (pySRU + WOFRY and SRW). In summary, these three methods proceed as follows:

(i) A simple solution (`internal' method) is based on reusing the original far-field (N_E, N_θ, N_φ) array. Without interpolation, we calculate the list of coordinates (x, y) for the available points. They cover one quadrant. These points are mirrored to cover the four quadrants using the symmetry of the radiated far field. Then the list of points and fields are propagated using a 2D integral propagator into a line on the y axis. We cannot use Fourier transforms as the points are not over a Cartesian homogeneous grid. Limiting the points in the backpropagated plane to a line makes reasonable calculation times. Then, the intensity maps are calculated for each energy and then added. Supposing axial symmetry, we can then create the 2D intensity map from the calculated 1D intensity distribution. It is from this 2D distribution map that the ray's coordinates (x, y) are sampled.

(ii) A second method uses the external packages pySRU and WOFRY. The far field is re-calculated in a Cartesian grid using pySRU and the 2D wavefront is backpropagated using two Fourier transforms with WOFRY. The obtained backpropagated fields (one for each photon energy) are used to compute the intensities. The intensity map for the whole energy interval (obtained by adding the intensities of the individual photon energies) is used to sample the (x, y) coordinates of the rays.

(iii) A third solution, similar to the previous one, uses SRW (Chubar & Elleaume, 1998) for calculating the source and backpropagating it. All calculations are done in Cartesian coordinates.

Footnotes

¹In SHADOW4 the axes naming is different: the propagation is along y, and at the source level x is horizontal and z vertical.

²Note that Tanaka & Kitamura (2009) have an additional factor of two. This is related to the mentioned discussion on calculating the size using equations (10) or (11). In our case, it is included in our equation (12), therefore Q_s → 1 for zero energy spread x → 0.