Memory and CPU efficient computation of the Fresnel free-space propagator in Fourier optics simulations

Oleg Chubar; Rafael Celestre

doi:10.1364/OE.27.028750

1. Introduction

Wave-optics calculation methods are becoming increasingly important for a large variety of applications in many branches of modern science and technology requiring diffraction-limited spatial resolution or/and dealing with miscellaneous wave-optics phenomena. Modern high-brightness and high-coherence light source facilities − low-emittance synchrotrons and linear accelerator based X-ray Free-Electron Lasers (XFELs) − are among the radiation sources requiring this type of calculations for their design, commissioning, and experimental applications. Special features of these sources are: small radiation wavelength (typically belonging to X-ray spectral range), partial transverse coherence, considerable deviation of radiation beam shapes from Gaussian for storage rings, short (femtosecond-range) duration and often complicated time structure or radiation pulses, requiring time-/frequency-dependent calculations for XFELs. In this context, wave-optics based simulations of radiation propagation through beamlines are often memory-hungry, and in some cases impossible to accomplish by using the conventional Fourier optics based methods (cf. §5 in [1]) on personal computers or even on servers.

Several software packages are available and routinely used throughout different light source facilities, exploring different approaches to numerical wave-optics calculations for beamline design [2–7]. In this paper, we describe some features and applications of “Synchrotron Radiation Workshop” (SRW) [3,4], the code allowing for both the accurate calculation of different types of spontaneous Synchrotron Radiation (SR) [8] and the wave-optics based simulation of radiation propagation through optical systems of beamlines in the infrared, ultra-violet, soft and hard X-ray spectral range at modern synchrotron light source facilities [9–13].

The method described in Sect. 3 of this paper was implemented in the SRW code in 2007 and allowed at that time to perform 3D time-/frequency-dependent wave propagation simulations for an XFEL pulse on a regular personal computer [14]. After this first application, it was successfully used for many other simulations for light sources, and it is now the most frequently used method in SRW for the simulation of electric field propagation through drift spaces [15]. Its application typically results in considerable memory savings and CPU performance increases compared to calculations based on the standard Fourier optics Fresnel propagator (see e.g. §5 of [1] or [16,17]). Since we never described this method (as implemented in SRW) in detail, we decided to give its formal description, together with representative application examples from synchrotron light sources. During the review process of this paper, we were informed that the method in question was published in [18] (though it still remains to be not very well known). Before describing this method, that we reference as the method “with the analytical treatment of the quadratic phase terms”, we remind the basics of the conventional Fresnel free-space propagator for comparison. Besides the applications for synchrotron light sources, the method can be used for other radiation sources and optical systems as long as the assumptions leading to the paraxial approximation are met (see e.g. §4 in [1] or [19,20]). We note however that special precautions should be taken when using this method for a beam considered at its waist before or after the propagation, as discussed at the end of Sect. 3.

2. Fresnel free-space propagator

Following the standard scalar diffraction theory (§3 in [1] or §8 in [21]), the free-space propagation of a transverse component of the frequency-domain complex electric field U can be described by the Rayleigh-Sommerfeld diffraction formula (cf. Eq. (3–42) in [1]) which is, in the small-angle approximation:

(1)$$U({{\bf r}_2}) \approx \frac{k}{{2\pi i}}\int_{{\Sigma _1}} {U({{\bf r}_1})|{{\bf r}_2} - {{\bf r}_1}{|^{ - 1}}\exp (ik|{{\bf r}_2} - {{\bf r}_1}|)d{\Sigma _1}} ,$$

where the integration is performed over the surface ${\Sigma _1}$, ${{\bf r}_1} = ({x_1},{y_1},{z_1})$ is a point on that surface, ${{\bf r}_2} = ({x_2},{y_2},{z_2})$ is an observation point, and k is the radiation wavenumber. If ${\Sigma _1}$ is a plane (e.g. normal to optical axis Z), and ${{\bf r}_2}$ belongs to another plane located at a distance ${z_2} - {z_1} = L$ from ${\Sigma _1}$, then $|{{\bf r}_2} - {{\bf r}_1}|= {[{L^2} + {({x_2} - {x_1})^2} + {({y_2} - {y_1})^2}]^{1/2}},$ $d{\Sigma _1} = d{x_1}d{y_1},$ and Eq. (1) is a convolution-type relation that can be computed numerically using 2D FFTs, even without any further analytical transformations. However, it can be more convenient to use the quadratic approximation of the exponent phase in Eq. (1), which can be obtained by applying a Taylor series expansion assuming $L > > |{x_2} - {x_1}|$, $L > > |{y_2} - {y_1}|$:

(2)$${U_2}({x_2},{y_2}) \approx \frac{{k\exp (ikL)}}{{2\pi iL}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{U_1}({x_1},{y_1})\exp \left\{ {\frac{{ik}}{{2L}}[{{({x_2} - {x_1})}^2} + {{({y_2} - {y_1})}^2}]} \right\}d{x_1}d{y_1}} } ,$$

where ${U_1}({x_1},{y_1}) \equiv U({x_1},{y_1},{z_1})$ and ${U_2}({x_2},{y_2}) \equiv U({x_2},{y_2},{z_1} + L)$. Equation (2) is often referred to as the Fresnel free-space propagator. The kernel of Eq. (2) and its Fourier transform from the coordinate $(x,y)$ to the angular $({\theta _x},{\theta _y})$ domain are:

(3)$$\begin{array}{c} K(x,y) = \frac{k}{{2\pi iL}}\exp \left[ {\frac{{ik}}{{2L}}({x^2} + {y^2})} \right],\\ \tilde{K}({\theta _x},{\theta _y}) = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {K(x,y)\exp [ - ik({\theta _x}x + {\theta _y}y)]dxdy = \exp \left[ { - \frac{{ikL}}{2}(\theta_x^2 + \theta_y^2)} \right]} } . \end{array}$$

By applying the convolution theorem, ${U_2}$ can be calculated numerically using the inverse Fourier transform:

(4)$$\begin{array}{c} {U_2}(x,y) \approx \frac{{{k^2}\exp (ikL)}}{{{{(2\pi )}^2}}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{{\tilde{U}}_1}({\theta _x},{\theta _y})\exp \left[ { - \frac{{ikL}}{2}(\theta_x^2 + \theta_y^2)} \right]} } \\ \times \exp [ik(x{\theta _x} + y{\theta _y})]d{\theta _x}d{\theta _y}, \end{array}$$

where ${\tilde{U}_1}({\theta _x},{\theta _y}) = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{U_1}(x,y)\exp [ - ik({\theta _x}x + {\theta _y}y)]dxdy} } $ is the Fourier transform of the electric field (i.e. the angular domain electric field) before the propagation (cf. §4.2 in [1]).

The numerical computation of Eq. (2) via Eq. (4) requires sufficiently dense sampling of the initial electric field ${U_1}$ versus the transverse coordinates x₁ and y₁ to ensure that all its oscillations within the ranges of x₁ and y₁, where the amplitude is significant, are resolved. If these electric field oscillations are due to the quadratic phase terms, then the required number of equidistant points over each of the two coordinates can be estimated as:

(5)$$N = [{{p\Delta {\theta ^2}R} \mathord{\left/ {\vphantom {{p\Delta {\theta^2}R} \lambda }} \right.} \lambda }],$$

where ${\Delta }\theta$ is the angular size of a wavefront that is assumed to be symmetric with respect to the optical axis in the horizontal or vertical plane, R is its corresponding radius of curvature, $\lambda$ is the radiation wavelength, p is a “sampling parameter”: p = 1 corresponds to nominal sampling, when phase advance between any adjacent mesh points is not larger than π, p < 1 corresponds to under-sampling, and p > 1 to over-sampling regime. Note that Eq. (5) coincides with the expression for the Fresnel number for a wavefront observed at distance R at p = 1 (disregarding the square brackets meaning integer part of a number). In practice, however, the number of points in numerical propagation calculations may need to be larger than what follows from Eq. (5) to have sufficient sampling of the electric field after the propagation. As we will illustrate in Sect. 4, the requirement for large numbers of points versus the horizontal and vertical positions represents a significant disadvantage of the numerical calculation method based on Eq. (4) in wave propagation simulations.

3. Analytical treatment of the quadratic radiation phase terms

Assume that the electric field before the propagation has quadratic terms in its phase, defined by the radii of the wavefront curvature in the horizontal and vertical planes R_x, R_y and the transverse ordinates of the center point (x₀, y₀):

(6)$${U_1}({x_1},{y_1}) = {F_1}({x_1},{y_1})\exp \left\{ {\frac{{ik}}{2}\left[ {\frac{{{{({x_1} - {x_0})}^2}}}{{{R_x}}} + \frac{{{{({y_1} - {y_0})}^2}}}{{{R_y}}}} \right]} \right\}.$$

The substitution of Eq. (6) into Eq. (2) gives:

(7)$$\begin{array}{c} {U_2}({x_2},{y_2}) \approx \frac{{k\exp (ikL)}}{{2\pi iL}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{F_1}({x_1},{y_1})} } \\ \times \exp \left\{ {\frac{{ik}}{2}\left[ {\frac{{{{({x_1} - {x_0})}^2}}}{{{R_x}}} + \frac{{{{({y_1} - {y_0})}^2}}}{{{R_y}}} + \frac{{{{({x_2} - {x_1})}^2} + {{({y_2} - {y_1})}^2}}}{L}} \right]} \right\}d{x_1}d{y_1}, \end{array}$$

which, after re-grouping the terms in the exponent, can be re-written as:

(8)$${U_2}({x_2},{y_2}) = {F_2}({x_2},{y_2})\exp \left\{ {\frac{{ik}}{2}\left[ {\frac{{{{({x_2} - {x_0})}^2}}}{{{R_x} + L}} + \frac{{{{({y_2} - {y_0})}^2}}}{{{R_y} + L}}} \right]} \right\}$$

with:

(9)$$\begin{array}{c} {F_2}({x_2},{y_2}) \approx \frac{{k\exp (ikL)}}{{2\pi iL}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{F_1}({x_1},{y_1})} } \\ \times \exp \left\{ {\frac{{ik}}{{2L}}\left[ {\frac{{{R_x} + L}}{{{R_x}}}{{\left( {\frac{{{R_x}{x_2} + L{x_0}}}{{{R_x} + L}} - {x_1}} \right)}^2} + \frac{{{R_y} + L}}{{{R_y}}}{{\left( {\frac{{{R_y}{y_2} + L{y_0}}}{{{R_y} + L}} - {y_1}} \right)}^2}} \right]} \right\}d{x_1}d{y_1}. \end{array}$$

Like Eq. (2), Eq. (9) is a convolution-type integral, though it applies to the function F₁ that can be considered as the initial electric field after the “subtraction” of its quadratic phase terms according to Eq. (6). Note that F₂ can be considered as a propagated electric field with its (evolved) quadratic terms subtracted, in line with Eq. (8).

The kernel of the convolution-type integral in Eq. (9) and its Fourier transform are:

(10)$$\begin{array}{c} G(x,y) = \frac{k}{{2\pi iL}}\exp \left[ {\frac{{ik}}{{2L}}\left( {\frac{{{R_x} + L}}{{{R_x}}}{x^2} + \frac{{{R_y} + L}}{{{R_y}}}{y^2}} \right)} \right],\\ \tilde{G}({\theta _x},{\theta _y}) = {\left[ {\frac{{{R_x}{R_y}}}{{({R_x} + L)({R_y} + L)}}} \right]^{1/2}}\exp \left[ { - \frac{{ikL}}{2}\left( {\frac{{{R_x}}}{{{R_x} + L}}\theta_x^2 + \frac{{{R_y}}}{{{R_y} + L}}\theta_y^2} \right)} \right]. \end{array}$$

After applying the convolution theorem, we obtain:

(11)$$\begin{array}{c} {F_2}(x,y) \approx \frac{{{k^2}\exp (ikL)}}{{{{(2\pi )}^2}}}{\left[ {\frac{{{R_x}{R_y}}}{{({R_x} + L)({R_y} + L)}}} \right]^{1/2}}\int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{{\tilde{F}}_1}({\theta _x},{\theta _y})} } \\ \times \exp \left[ { - \frac{{ikL}}{2}\left( {\frac{{{R_x}}}{{{R_x} + L}}\theta_x^2 + \frac{{{R_y}}}{{{R_y} + L}}\theta_y^2} \right)} \right]\exp \left[ {ik\left( {\frac{{{R_x}x + L{x_0}}}{{{R_x} + L}}{\theta_x} + \frac{{{R_y}y + L{y_0}}}{{{R_y} + L}}{\theta_y}} \right)} \right]d{\theta _x}d{\theta _y}, \end{array}$$

where ${\tilde{F}_1}({\theta _x},{\theta _y}) = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } {{F_1}(x,y)\exp [ - ik({\theta _x}x + {\theta _y}y)]dxdy} } $. To obtain F₂ from F₁, two 2D FFTs must be computed, in full correspondence with Eq. (4).

The quadratic terms are in many cases the leading terms of the complex electric field phase. These terms play a fundamental role in the propagation of a radiation beam (e.g. by defining whether the beam is diverging or converging, and at which rate). These terms are typically responsible for the phase growth as a function of horizontal and vertical position, and for the related oscillatory behavior of the complex electric field outside of the beam waist. To ensure the accuracy of digital Fourier transforms of the complex electric field, it must be sampled densely enough to resolve all its oscillations, and this typically drives the memory requirements in the standard Fresnel propagator computation method, see Eq. (5). Therefore, the possibility of “subtraction” of the quadratic phase terms before applying FFTs, without compromising the accuracy of the Fourier optics based calculations in the paraxial approximation, means a considerable reduction of memory requirements and increase of CPU-efficiency of the numerical calculations, which is particularly rewarding in the 2D case considered in this paper.

It is important to emphasize one more feature of the calculation method defined by Eqs. (6)–(11). As follows from Eq. (9), the convolution takes place with respect to the scaled transverse coordinates after the propagation, ${{({R_x}{x_2} + L{x_0})} \mathord{\left/ {\vphantom {{({R_x}{x_2} + L{x_0})} {({R_x} + L)}}} \right.} {({R_x} + L)}}$ and ${{({R_y}{y_2} + L{y_0})} \mathord{\left/ {\vphantom {{({R_y}{y_2} + L{y_0})} {({R_y} + L)}}} \right.} {({R_y} + L)}}.$ This means that, after formally applying Eq. (11), without any additional re-sampling or interpolation, we obtain the propagated electric field sampled over the transverse ranges that are scaled proportionally to the geometrically-estimated size of the propagated radiation beam: $\Delta {x_2} = {{\Delta {x_1}({R_x} + L)} \mathord{\left/ {\vphantom {{\Delta {x_1}({R_x} + L)} {{R_x}}}} \right.} {{R_x}}},$ $\Delta {y_2} = {{\Delta {y_1}({R_y} + L)} \mathord{\left/ {\vphantom {{\Delta {y_1}({R_y} + L)} {{R_y}}}} \right.} {{R_y}}},$ where $\Delta {x_1},\Delta {y_1}$ and $\Delta {x_2},\Delta {y_2}$ are the horizontal and vertical grid dimensions before and after the propagation. This feature is illustrated in Fig. 1. It enables additional economy of memory and CPU resources at wave propagation simulations.

Fig. 1. Illustration of the scaling of transverse coordinates at the wave propagation over a distance L when using the method with the analytical treatment of the quadratic phase terms.

Download Full Size | PDF

The numerical propagator with the analytical treatment of the quadratic phase terms described above does have limitations that need to be discussed, however. Equation (11) has singularities at ${R_x} + L = 0$ and ${R_y} + L = 0$, and therefore cannot be applied for simulation of the electric field propagation directly to a waist if R_x and R_y are the exact radii of wavefront curvature before the propagation. Also, the application of Eq. (11) can be inefficient if the initial electric field is considered at a radiation beam waist, where there are no significant quadratic phase terms. To circumvent these issues, several simple techniques can be used. For example, one can use in Eqs. (6) and (11) values of R_x and R_y $({R_x} \ne 0,{R_y} \ne 0)$ that somewhat deviate from the exact wavefront radii of curvature, making sure that ${R_x} + L \ne 0,{R_y} + L \ne 0.$ Such adjustment of R_x and R_y values can be done automatically in software (this method was implemented in the SRW code and applied in the simulations described in Sect. 4). Alternatively, one can use the standard method based on Eq. (4) for the propagation over a small distance in the vicinity of the waist (where the quadratic phase terms are small), and apply the method based on Eq. (11) at a larger distance from the waist. If the propagation to/from waist takes place over a very large distance L, one can use another special Fourier optics propagator based on the far-field Fraunhofer approximation requiring one FFT (cf. §4.3 and §6.2 in [1]).

4. Economy of memory and CPU resources

Numerical simulations of electric field propagation through drift spaces in the paraxial approximation using the method with the analytical treatment of the quadratic phase terms (see Eqs. (6)–(11)) typically require much less memory and CPU time than the simulations based on the standard method (see Eqs. (2)–(4)). In this section, we illustrate this using two representative wave propagation examples related to Undulator Radiation (UR) and bending magnet SR sources. In both cases, the initial radiation was calculated using the same general method based on the retarded potentials, see Eqs. (1) and (2) in [8]. The emission and wavefront propagation calculations were performed with the SRW code.

4.1 Undulator radiation example

In the first example, we simulated propagation of ∼15 keV UR (at 9th spectral harmonic, generated by a 3 GeV electron beam in a 20-mm period 3-m long planar undulator) through a hypothetic thin lens with 10 m focal length, located at 30 m from the center of the undulator, and through a drift space, see Fig. 2. The length of the drift was equal to 14.5 m in one case (for the propagation to a position 0.5 m before the geometrical image plane), and 15 m in another case (for the propagation exactly to the geometrical image plane). The simulations were performed both using the standard Fresnel free-space propagation method and the method with the analytical treatment of the quadratic phase terms. The results of these simulations are presented in Fig. 3, where: the distributions of the intensity and phase of the initial UR are shown in Fig. 3(a), the distributions at 0.5 m before the image plane in Fig. 3(b), and the ones in the image plane in Fig. 3(c). The continuous and smooth phase cuts vs. horizontal and vertical positions in the graphs on the right in Figs. 3(a) and 3(b) were reconstructed from the computed principal phase values, whereas the phase cuts in Fig. 3(c) show the principal phase value without any further processing.

Fig. 2. Optical scheme used in the UR simulations presented in Fig. 3.

Download Full Size | PDF

Fig. 3. Intensity and phase distributions of undulator radiation by a filament 3 GeV electron beam at ∼15 keV photon energy (at ∼6 eV “red shift” with respect to the resonant photon energy of 9th spectral harmonic of a 20-mm period 3 m long planar undulator): (a) at 30 m from the undulator center; (b) propagated to a position 0.5 m before the image plane of the lens with 10 m focal length; (c) propagated to the image plane (see Fig. 2). Graphs in the center show cuts, by the horizontal and vertical mid-planes, of the intensity distributions presented by image plots on the left. Dots in these graphs correspond to the calculations with the standard Fresnel free-space propagator, and lines to the calculations with the “special” free-space propagator making the analytical treatment of the quadratic phase terms. Graphs on the right show the corresponding horizontal and vertical cuts of the radiation phase; the first of these graphs (case a) contains, besides the total phase, sampled as for the “special” propagator, also the residual phase after the subtraction of the quadratic phase terms; the other graphs (cases b, c) illustrate the total phase after the propagation using the two different methods.

Download Full Size | PDF

The initial UR was calculated with a ∼6 eV “red shift” with respect to the on-axis resonant photon energy of the harmonic, which results in a “doughnut-shape” intensity distribution with a dip in the center, see the image plot on the left and graph in the center in Fig. 3(a). The radiation flux is higher in such a mode compared to the mode corresponding to the exact on-axis resonant photon energy [3,22]. Note that despite of a considerable non-homogeneity of the UR intensity distribution (and the amplitude of the complex electric field), the residual radiation phase, remaining after the subtraction of the quadratic terms, is nearly constant over the area where the amplitude is most significant (see graph on the right in Fig. 3(a)). This UR feature was discussed in detail in [22].

The initial wavefront of ∼30 × 30 µrad² angular size was calculated on a grid with 80 points in the horizontal and vertical direction, though Eq. (5) suggests 326 points in this case. For the propagation using the standard method, the initial electric field was re-sampled by interpolation to reduce the grid step size and increase the number of points by a factor of 15 in each direction. This factor was chosen to obtain nominally resolved intensity distributions after the propagation. The quadratic phase terms were subtracted before the re-sampling and added back after the re-sampling, which is the standard procedure for re-sampling electric fields in SRW.

As this was discussed in Sect. 3, in the propagation method with the analytical treatment of the quadratic phase terms, transverse ranges of the computation grid change according to the rules of geometrical optics (see Fig. 1), which helps to reduce possible losses of information at the propagation simulations. However, since coherent radiation propagates according to the laws of wave optics, rather than the laws of geometrical optics, this geometrical scaling is still not sufficient in many cases, including the case under consideration. Due to this reason, the range of the initial wavefront had to be increased by factor 1.4 in each direction by “padding zeros” (without changing the step sizes of the initial grid) in order to obtain the intensity distributions of the propagated radiation fully covered by the final grid. The simulation results obtained using the method with the analytical treatment of the quadratic phase terms are shown by curves marked as “special” in Figs. 3(b)–3(c).

The calculation using the standard method required ∼23 MB of memory (total for the two transverse field components) and took ∼700 ms on a laptop computer equipped with Intel Core i7 7^th generation processor operating at 2.8 GHz. The calculation of the propagation with the analytical treatment of the quadratic phase terms required only ∼0.2 MB of memory and took less than 10 ms on the same computer. We note that, despite the considerably smaller amount of memory and CPU time spent on the calculation with the latter method, the resulting intensity distributions appeared to be more densely sampled than at the calculation using the standard method.

4.2 Bending magnet SR example

In the second example, propagation of 10 keV SR generated by 3 GeV electron beam in 0.4 T bending magnet was simulated using the method with the analytical treatment of the quadratic phase terms. The optical scheme in these simulations was the same as in the previous example, see Fig. 4; the simulation results are presented in Fig. 5. The initial wavefront was calculated within the 30 mm horizontal and 15 mm vertical coordinate range at 30 m distance from the tangential point on the circular electron trajectory in the bending magnet, see Fig. 5(a). These ranges make angular dimensions of the initial wavefront relatively large (e.g. as compared to the previous example): 1 mrad in the horizontal and 0.5 mrad in the vertical plane, which corresponds to typical acceptances of a hard X-ray bending magnet or wiggler beamline at a synchrotron light source facility.

Fig. 4. Optical scheme used in the bending magnet SR simulations presented in Fig. 5.

Download Full Size | PDF

Fig. 5. Intensity and phase distributions of SR from a 0.4 T bending magnet generated by a filament 3 GeV electron beam at 10 keV photon energy: (a) at 30 m from tangential source point in the bending magnet; (b) propagated to a position 10 mm before the geometrical image plane of the lens with 10 m focal length; (c) propagated to the geometrical image plane (see Fig. 4). Graphs in the center show cuts, by the horizontal and vertical planes, of the intensity distributions presented by image plots on the left. Graphs on the right show the corresponding horizontal and vertical cuts of the radiation phase. The calculations were performed by making use of the free-space propagator with the analytical treatment of the quadratic phase terms.

Download Full Size | PDF

The bending magnet SR intensity distribution (see the image plot on the left and graph in the center in Fig. 5(a)) is practically uniform in the horizontal direction (due to a relatively large distance from the tangential source point) and has a Gaussian-like profile in the vertical direction, which corresponds to the emission at photon energies higher than the critical energy of the bending magnet SR. The phase of the initial wavefront is very strongly dominated by the quadratic phase terms, see the upper graph on the right in Fig. 5(a). However, as different from the previous example, the residual phase after the subtraction of the quadratic terms is not constant in this example, even though this can’t be seen in the same linear scale together with the total phase. This residual phase has cubic dependence on the horizontal position (and angle), as illustrated in the zoomed graph for small phase values, located under the main phase graph in Fig. 5(a) on the right. This feature of the phase of bending magnet SR was discussed in other works, see e.g. [22], where an analytical approximation for the residual phase was given by Eq. (8).

Figures 5(b) and (c) show the intensity distributions at 10 mm before and in the geometrical image plane, respectively (see Fig. 4). The intensity distribution before the image plane is dominated by the diffraction at the input slit, defined by the dimensions of the initial wavefront, with modulation over the horizontal position. The contribution of the residual phase to the total radiation phase in the horizontal direction can be seen in the graph on the right in Fig. 5(b), since the quadratic phase terms at that location are much smaller than in the initial wavefront (Fig. 5(a)). The intensity distribution in the image plane has a specific shape of a focused bending magnet SR, strongly affected by this residual phase [22].

The electric field of the initial wavefront was calculated on a grid with 294 × 78 points versus horizontal and vertical position, and re-sampled on a finer grid to reduce the step sizes down by factors 50 in the horizontal and 25 in the vertical direction (the re-sampling was made after subtraction of the quadratic phase terms). Even after this re-sampling, the numbers of points in the wavefront (14700 × 1950) were much smaller (by factor of ∼20 in the horizontal direction and by factor of ∼40 in the vertical direction) than the nominal numbers of points required for the standard Fresnel free-space propagator, as estimated by Eq. (5). The simulation of propagation of the two transverse electric field components required ∼460 MB of memory and took ∼14 s on the same computer that was used in the previous example. The corresponding calculation using the standard method would require at least ∼370 GB of memory and would have lasted much longer. The actual calculation with the standard free-space propagator was not performed in this example because of lack of the required hardware.

5. Conclusion

The numerical free-space Fourier optics propagator with the analytical treatment of the quadratic phase terms allows for considerable economy of memory and CPU resources as compared to the standard Fresnel free-space propagator. This economy was shown to be between a factor of 100 and 1000 both in terms of the memory and the execution time in the synchrotron light source wave-optics simulation examples related to the hard X-ray spectral range. The method was presented for fully-coherent radiation beams, however, it is extensively applied in the SRW code for more complicated and more practically-important partially-coherent simulations for storage ring sources, and for 3D time-/frequency-dependent simulations for XFELs [23,24]. Such simulations typically require more computational resources than the fully-coherent / steady-state ones, making the performance gains offered by this method extremely valuable. We also note that the very significant memory economy offered by the method with the analytical treatment of the quadratic phase terms facilitates execution of the wave propagation simulations on ultra-fast GPU based computer systems (that usually have smaller memory than the CPU based ones), resulting in further significant boost of the computation performance. This may open up new possibilities for applications of these simulations, such as real-time synchrotron beamline and XFEL instrument tuning for user experiments and fast experimental data processing.

Funding

Office of Science (DE-SC0011237, DE-SC0012704, FWP PS-017).

Acknowledgments

Authors are grateful to C. Jacobsen, S. A. Syed (Northwestern University), L. Rebuffi (ANL / APS), M. S. del Rio, R. Barrett (ESRF), J. Krzywinski (SLAC), M. E. Couprie (SOLEIL), B. Nash (RadiaSoft) and R. Coles (BNL / NSLS-II) for encouragement, comments, and help.

References

1. J. W. Goodman, Introduction to Fourier Optics (W. H. Freeman and Company, 2017), 4th ed.

2. J. Bahrdt, “Wave-front propagation: design code for synchrotron radiation beam lines,” Appl. Opt. 36, 4367–4381 (1997). [CrossRef]

3. O. Chubar and P. Elleaume, “Accurate and efficient computation of synchrotron radiation in the near field region,” Part. accelerator. Proceedings, 6th Eur. conference, EPAC’98C980622, 1177–1179 (1998).

4. Open-source repository of “Synchrotron Radiation Workshop” code: https://github.com/ochubar/SRW

5. D. Spiga and L. Raimondi, “X-ray optical systems: from metrology to point spread function,” Proc. SPIE 9209, 920901 (2014). [CrossRef]

6. X. Shi, R. Reininger, M. Sanchez del Rio, and L. Assoufid, “A hybrid method for X-ray optics simulation: combining geometric ray-tracing and wavefront propagation,” J. Synchrotron Radiat. 21, 669–678 (2014). [CrossRef]

7. R. Chernikov and K. Klementiev, “Recent progress of the xrt: ray tracing and wave propagation toolkit (conference presentation),” Proc. SPIE 10388, 1038806 (2017). [CrossRef]

8. O. Chubar, “Precise computation of electron-beam radiation in nonuniform magnetic fields as a tool for beam diagnostics,” Rev. Sci. Instrum. 66, 1872–1874 (1995). [CrossRef]

9. P. Dumas, F. Polack, B. Lagarde, O. Chubar, J. Giorgetta, and S. Lefrançois, “Synchrotron infrared microscopy at the French Synchrotron Facility SOLEIL,” Infrared Phys. Technol. 49, 152–160 (2006). [CrossRef] International Workshop on Infrared Microscopy and Spectroscopy with Accelerator-Based Sources.

10. O. Chubar, J. Susini, M. Cotte, F. Polack, B. Lagarde, K. Scheidt, P. Elleaume, and P. Dumas, “Simulation and optimization of synchrotron infrared micro-spectroscopic beamlines using wave optics computation: ESRF and SOLEIL cases,” AIP Conf. Proc. 879, 607–610 (2007). [CrossRef]

11. A. Giuliani, F. Jamme, V. Rouam, F. Wien, J.-L. Giorgetta, B. Lagarde, O. Chubar, S. Bac, I. Yao, S. Rey, C. Herbeaux, J.-L. Marlats, D. Zerbib, F. Polack, and M. Réfrégiers, “DISCO: a low-energy multipurpose beamline at synchrotron SOLEIL,” J. Synchrotron Radiat. 16, 835–841 (2009). [CrossRef]

12. N. Canestrari, V. Bisogni, A. Walter, Y. Zhu, J. Dvorak, E. Vescovo, and O. Chubar, “Wavefront propagation simulations for a UV/soft x-ray beamline: Electron Spectro-Microscopy beamline at NSLS-II,” Proc. SPIE 9209, 920901 (2014). [CrossRef]

13. O. Chubar, Y. S. Chu, K. Kaznatcheev, and H. Yan, “Application of partially coherent wavefront propagation calculations for design of coherence-preserving synchrotron radiation beamlines,” Nucl. Instrum. Methods Phys. Res., Sect. A 649(1), 118–122 (2011). [CrossRef]

14. O. Chubar, M.-E. Couprie, M. Labat, G. Lambert, F. Polack, and O. Tcherbakoff, “Time-dependent FEL wavefront propagation calculations: Fourier optics approach,” Nucl. Instrum. Methods Phys. Res., Sect. A 593, 30–34 (2008). [CrossRef] FEL Frontiers 2007.

15. O. Chubar, “High-accuracy partially- and fully-coherent wavefront propagation calculations for storage ring and FEL sources,” presentation at Workshop on Softw. for Opt. Simulations, Trieste, Italy, Oct. 3–7 (2016). https://www.elettra.trieste.it/Conferences/2016/SOS/uploads/Main/Abstract1226.pdf.

16. D. P. Kelly, “Numerical calculation of the Fresnel transform,” J. Opt. Soc. Am. A 31, 755–764 (2014). [CrossRef]

17. M. Hillenbrand, D. P. Kelly, and S. Sinzinger, “Numerical solution of nonparaxial scalar diffraction integrals for focused fields,” J. Opt. Soc. Am. A 31, 1832–1841 (2014). [CrossRef]

18. U. Vokinger, Propagation, Modification and Analysis of Partially Coherent Light Fields, UFO-Dissertation (UFO, Atelier für Gestaltung und Verlag, 2000).

19. W. H. Southwell, “Validity of the fresnel approximation in the near field,” J. Opt. Soc. Am. 71, 7–14 (1981). [CrossRef]

20. G. W. Forbes, “Validity of the Fresnel approximation in the diffraction of collimated beams,” J. Opt. Soc. Am. A 13, 1816–1826 (1996). [CrossRef]

21. M. Born, E. Wolf, A. B. Bhatia, P. C. Clemmow, D. Gabor, A. R. Stokes, A. M. Taylor, P. A. Wayman, and W. L. Wilcock, Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light (Cambridge University, 1999), 7th ed.

22. O. Chubar, P. Elleaume, and A. Snigirev,“Phase analysis and foFundincusing of synchrotron radiation,” Nucl. Instrum. Methods Phys. Res., Sect. A 435, 495–508 (1999). [CrossRef]

23. O. Chubar, L. Berman, Y. S. Chu, A. Fluerasu, S. Hulbert, M. Idir, K. Kaznatcheev, D. Shapiro, Q. Shen, and J. Baltser, “Development of partially-coherent wavefront propagation simulation methods for 3rd and 4th generation synchrotron radiation sources,” Proc. SPIE 8141, 814107 (2011). [CrossRef]

24. M. S. Rakitin, P. Moeller, R. Nagler, B. Nash, D. L. Bruhwiler, D. Smalyuk, M. Zhernenkov, and O. Chubar, “Sirepo: an open-source cloud-based software interface for X-ray source and optics simulations,” J. Synchrotron Radiat. 25, 1877–1892 (2018). [CrossRef]

Memory and CPU efficient computation of the Fresnel free-space propagator in Fourier optics simulations

Abstract

1. Introduction

2. Fresnel free-space propagator

3. Analytical treatment of the quadratic radiation phase terms

4. Economy of memory and CPU resources

4.1 Undulator radiation example

4.2 Bending magnet SR example

5. Conclusion

Funding

Acknowledgments

References

Cited By

Figures (5)

Equations (11)

Optics Express