Finite difference methods for stationary and time-dependent X-ray propagation

Lars Melchior; Tim Salditt

doi:10.1364/OE.25.032090

1. Introduction

X-ray optics has undergone substantial progress over the last decade, from advanced sources of ultra-fast and/or highly brilliant radiation, improved optical fabrication, to novel schemes of time-resolved and coherent diffractive imaging [1–3]. With the augmented experimental capabilities the demand for accurate and high performance field simulations and optical design has increased. However, to simulate and predict the propagation of x-ray fields in optical devices and objects is an extremely challenging task, as one quickly realizes when considering the requirement of sampling of macroscopic objects or optical devices on scales of typical x-ray wavelength λ. At the same time, x-ray optics is characterized by an intrinsically weak interaction with matter, and many relevant problems address forward-directed paraxial propagation involving only small angles with respect to the optical axis. Further, as the x-ray index of refraction, conventionally written as n = 1 − δ + iβ, deviates only very weakly from vacuum with δ, β ≪ 1, parabolic approximations to the stationary wave equation and envelope approximations can be used. For this reason, finite difference (FD) calculations [4, 5] are well suited for many x-ray optical propagation problems, solving approximative parabolic wave Eqs. by forward-stepping on a grid. FD simulations of the parabolic wave equation have been used for example to study focused fields of Fresnel zone plates (FZP) [6], and propagation of x-ray guided beams [7–12]. While most implementations use propagation in two dimensions (2D) as described in detail in [9], numerically schemes for more demanding propagation in three dimensions (3D) have also been presented [10]. Non-stationary time-dependent propagation (4D), on the other hand, has become increasingly relevant in view free electron laser (FEL), and high harmonics generation sources (HHG), and the need to design the corresponding experimental schemes of short pulses. While this challenge has been addressed based on a spectral decomposition of free space Fresnel propagation [13], as well as with advanced finite-difference time domain (FDTD) including the non-linear optical response [14], a simple time-dependent FD framework is still missing. It offers the advantage to compute pulse propagation efficiently in the limit of small λ up to macroscopic system size.

Optical simulation of x-ray pulses is required for example to study pulse dispersion in focusing optics such as FZP or compound refractive lenses (CRL), or in beamsplitters, apertures, monochromators and are also required to design and optimize ultra-fast coherent imaging experiments. A simple example of a device which is to date nearly impossible to study numerically is the propagation of pulses in x-ray waveguide chips, with a path length of a pulse in the chip extending several mm’s, while at the same time extremely high spatial and temporal sampling is required, for example to cover the modulations on nanometer scale in the lateral dimension, as well as the atto- to femto-second modulations of the envelope.

The scope of the present work is to extend finite-difference (FD) simulations of beam propagation to cover time-dependent problems by spectral decomposition. While we are fully aware that ultra-brilliant femtosecond FEL pulses offer unique opportunities for non-linear optics, surprisingly many problems of hard x-ray propagation can be still well treated by linear optical propagation of quasi-monochromatic light, down to ultra-short pulse duration. Conceptually simple, the main challenge is in the numerical performance and implementation. Note that efficient stationary numerical schemes are a prerequisite to perform calculations based on spectral decomposition, since numerical complexity increases by a factor given by the spectral sampling points. To this end, we present a set of FD and Fresnel different propagators for 2D, 3D and problems with circular symmetry (CS) in the plane orthogonal to the optical axis, all implemented in an open-source high-performance simulation toolbox for Python. We also present a propagator for paraxial propagation along arbitrary smooth optical paths, and present a unified implementation for FD and free space (Fesnel) propagation. We then validate the stationary propagators by comparison to analytical results, and compare the numerical performance and accuracy with multislice approaches [15]. Next, we extend the scope to time-dependent problems. Finally, we present pulse propagation in x-ray waveguides as a first example, and show that these can be considered as nearly dispersion free optics down to femtosecond pulses, while dispersion effects start to become visible in from of mode seperation for attosecond pulses.

The manuscript is organized as follows: After this introduction, section 2 briefly recapitulates the paraxial Helmholtz equation for scalar diffraction as the analytical foundation for the numerical propagation, which is presented in section 3 including several different versions of propagators. The numerical methods are validated by comparison to analytically solvable problems as well as to the multi-slice method. Section 4 then presents the numerical approach to time-dependent propagation, which is illustrated by first applications in section 5, before the manuscript closes with a discussion. For notational clarity and completeness, some fundamental aspects of scalar diffraction theory are included in the appendices, and code to generate all figures in this work is included in Code 1 (Ref. [16]).

2. The paraxial Helmholtz equation

To describe the propagation of paraxial beams in matter one usually solves the parabolic approximation of the elliptic Helmholtz equation, due to its lower order and corresponding advantages in numerics and specifying the boundary conditions. A common derivation is to insert the free-space paraxial beam ψ = u exp(−ikz) into the Helmholtz equation and to neglect second-order derivatives $| \frac{\partial^{2} u}{\partial z^{2}} | ≪ | k \frac{\partial u}{\partial z} |$ to obtain a parabolic wave equation for the envelope

\frac{\partial u}{\partial z} = \frac{1}{2 i k} Δ_{⊥} u - i k \frac{n^{2} - 1}{2} u .

This form of the paraxial wave equation is by now well established in x-ray optics ([6,8,10]). It is limited to refractive indices very close to one n ≃ 1, to keep the spatial modulation of the envelope sufficiently small.

Here we derive a more general, slightly different form of the paraxial wave equation, by applying the paraxial approximation in Fourier space, analogous to the approach used in [18] for nonlinear propagation problems. To this end, we take a two-dimensional Fourier transform with respect to the perpendicular space directions x⃗_⊥:= (x y)^T, and transform the scalar Helmholtz Eq. (11) to a second-order ordinary differential equation

\frac{\partial^{2} \tilde{ψ}}{\partial z^{2}} + (k^{2} n^{2} - k_{⊥}^{2}) \tilde{ψ} = (\frac{\partial^{2}}{\partial z^{2}} + β^{2}) \tilde{ψ} = 0,

where ψ̃ = ℱ_{x⃗_⊥}[ψ](k⃗_⊥) and

β : = \sqrt{k^{2} n^{2} - k_{⊥}^{2}}

. The differential operator may be written as the commutative product of two operators

(\frac{\partial}{\partial z} + i β) (\frac{\partial}{\partial z} - i β) \tilde{ψ} = (\frac{\partial}{\partial z} - i β) (\frac{\partial}{\partial z} + i β) \tilde{ψ} = 0 .

The differential equation is then fulfilled if ψ̃ solves either of the differential Eqs.

\frac{\partial \tilde{ψ}}{\partial z} = i β \tilde{ψ} \frac{\partial \tilde{ψ}}{\partial z} = - i β \tilde{ψ}

and the general solution is a linear combination of both solutions. In the context of paraxial beams, we see that the left equation applies to a beam where the wave vector k⃗ is oriented in −e⃗_z direction and the right equation applies to a beam where the wave vector is oriented in e⃗_z direction. For paraxial beams, where wave vectors nearly parallel to the z axis, the support of ψ̃ is contained in a small region where

k_{⊥}^{2} ≪ k^{2}

. We can therefore approximate the square root in β through its first-order Taylor series around

k_{⊥}^{2} = 0

β = \sqrt{k^{2} n^{2} - k_{⊥}^{2}} \approx k n - \frac{k_{⊥}^{2}}{2 k n} .

Substituting this into Eq. (3) (right) and transforming back to real space results in the following paraxial Helmholtz equation

\frac{\partial ψ}{\partial z} = \frac{1}{2 i k n} Δ_{⊥} ψ - i k n ψ .

The paraxial Helmholtz equation is valid for forward propagating paraxial beams in the linear optical regime. In view of higher numerical stability, we again formulate the paraxial Helmholtz equation for the envelope u, defined by ψ = u exp(−ikz),

\frac{\partial u}{\partial z} = \frac{1}{2 i k n} Δ_{⊥} u - i k (n - 1) u .

This equation has enhanced applicability compared to Eq. (1) as no additional constraints are set on the refractive index n. Note that in the range of hard x-rays Eq. (5) and Eq. (1) are essentially identical. Next we will address the numerical solution of Eq. (5) by finite difference calculations on a grid.

In the above derivation we tacitly assumed homogeneous matter, but the results carry over to sufficiently weakly varying index distributions n→ n(r⃗), as well as to piecewise homogeneous media, fulfilling the boundary conditions of a continuously differentiable scalar envelope.

3. Numerical propagation

The paraxial Helmholtz equation can be solved numerically using a variety of well known propagation algorithms. For small |n − 1|, a versatile and rather simple method is multislice Fresnel propagation [1,15], based on alternating steps of free space propagation and multiplication of the field with the projected optical indices. Complex distributions of matter are hence solved by subsequent calculation of thin slices, with free space propagation in between, efficiently computed by fast Fourier transformations. An advantage of the multi-slice approach is that large free space distances can be computed in a single step. Alternatively, we can use well-known Finite-Difference methods such as the Crank-Nicolson method [4,5] to directly solve Eq. (5) very efficiently and accurately on a suitably sampled grid. Using Finite-Differences and the paraxial Helmholtz equation there is no requirement for |n − 1| to be small. In the simplest formulation Dirichlet boundary conditions are used.

3.1. PyPropagate

We have implemented the multislice and finite-difference propagation algorithms in the numerical simulation framework PyPropagate (source code available at github [17], installation through pip). The algorithms are implemented for systems with two or three spacial dimensions and for systems with cylindrical symmetry around the z-axis. PyPropagate is implemented in the high-performance language C++ and uses Python 2.7 as a front-end scripting language. It is recommended to be used inside a jupyter notebook [19] environment. The object-oriented approach, that decouples simulation parameters from the propagation algorithms, allows users to set up simulations in very few lines of code and solve them using any of the implemented propagators. An overview of the implemented propagators and their time-complexity is shown in Tab. 1.

Table 1. Overview of the propagators implemented by PyPropagate. The algorithms are abbreviated as Crank-Nicolson method (CN), Alternating Direction Implicit method (ADI), Fast Fourier Transform based multislice method (FFT) and Discrete Hankel transform with resampling [20] based multislice method (DHT). The boundary conditions at the simulation box edges differ depending on the used propagation scheme. As an addition to the time complexity, an example runtime is shown, where the simulation box is divided into a grid with N_x = N_y = 1024 or N_r = 1024 and N_z = 1000 nodes. The example runtime has been determined on a 2.7 GHz dual-core consumer notebook from 2011.

View Table | View all tables in this article

3.2. Comparison of finite difference and Fresnel propagation

While both multislice and Finite-Difference methods essentially solve the paraxial Helmholtz equation Eq. (5), the methods’ accuracy behave very differently as a function the propagation step size Δz. In free space the accuracy of finite difference methods scales with Δz², while the accuracy of Fresnel methods is completely independent of Δz. We visualize the implications by comparing numerical results with the analytical solution of a paraxial Gaussian beam, which is an exact solution of the paraxial Helmholtz equation obtained by analytically solving the Fresnel convolution integral for a Gaussian-shaped initial field distribution. In three dimensions the Gaussian beam is

ψ_{Gauss} = ψ_{0} \frac{σ_{r}^{2}}{w_{r} {(z)}^{2}} exp (- i k n z - \frac{r^{2}}{2 w_{r} {(z)}^{2}}),

where

r = \sqrt{x^{2} + y^{2}}

,

w_{r} (z) = \sqrt{\frac{z}{i k n} + σ_{r}^{2}}

and ψ₀ is the beam amplitude in the focal spot at r, z = 0. The comparison is presented in Fig. 1 and Fig. 2. As expected, using large step sizes the Fresnel methods are still accurate, while finite difference methods produce inaccurate intensities and energy flows.

Fig. 1 Comparison of different propagation algorithms for free space propagation and propagation in cylindrical waveguides. Optical intensity is normalized by |ψ₀|² and color coded, while the time-averaged direction of energy flow is indicated by white streamlines in the lower half of the images. Left: For free space propagation, we set the initial field distribution of a two or three dimensional monochromatic Gaussian beam with 12 keV photon energy, a waist size of σ_r = 0.25 μm, at 1 mm from the focus. The step sizes are chosen as Δx = Δy = Δr = 10 nm and Δz = 400 μm. At these step sizes the results of Fresnel propagation agree well with the analytical solution, while finite difference algorithms produce significant numerical errors. Right: For the waveguide propagation we choose a plane wave initial condition. The simulated waveguide consists of a vacuum core and a Germanium cladding at room temperature. The diameter is chosen as 50 nm and the initial beam is Gaussian distributed with a waist size of σ_r = 100 nm. The analytical solution is obtained by eigenmode projection onto guiding modes that do not dissipate into the cladding. Therefore modes that rapidly dissipate into the cladding present in numerical simulations are not taken into account. The step sizes of the numerical propagators are Δx = Δy = Δr = 0.7 nm and Δz = 0.8 μm. The waveguide edge is indicated by a dashed white line in the upper image half. At these step sizes multislice methods produce significant numerical artifacts visible as horizontal lines that prevent accurate determination of energy flow. Finite difference methods produce significantly better results that are nearly identical with simulations of much smaller step sizes.

Download Full Size | PDF

Fig. 2 Convergence behaviour of different propagation algorithms for free space propagation (a)–(b) and propagation in matter (c)–(d). The largest magnitude of the difference of the propagated field envelope u_p from the analytical solution u_a normalized by the largest magnitude of the analytical solution is shown as a function of the number of propagation steps N_z. The step sizes are chosen as Δx = Δy = Δr ≈ 1 nm (a), 10 nm (b), 0.1 nm (c), 1 nm (d) and Δz = s_z/N_z with the propagation distance s_z = 20 mm (a)–(b), 1 mm (c)–(d). The Gaussian beam and waveguide parameters are identical to the parameters given in Fig. 1. For Gaussian beams the analytical solutions are exact, while the analytical solution for waveguides is based on eigenmode projection. To avoid bias caused by the non-guiding modes, for (c) and (d) the maximum is taken after propagating 0.5 mm inside the waveguide. While accuracy of Fresnel methods is independent of the propagation step size for free space propagation, they are outperformed by finite difference methods when propagating in matter. For cylindrically symmetric problems, the cylindrical symmetric propagators perform almost equally well as the 3D propagators at a significantly reduced computational cost.

Download Full Size | PDF

When the matter distribution is inhomogeneous, finite difference propagation produces much better results at larger step sizes than Fresnel propagation. This is shown in Fig. 1 and Fig. 2, where we compare numerical results of propagation in slab and cylindrical waveguides with analytical solutions gained by eigenmode projection [10].

Finally, we generalize paraxial propagation schemes to systems where the paraxial propagation direction contiguously changes with the propagation distance. In these systems the paraxial approximation is valid on nearby planes that are perpendicular to direction of propagation. After each propagation step the direction of propagation is slightly tilted to match the next plane. This tilting can be numerically implemented by resampling the envelope in the tilted plane. For a plane with a local coordinate system system (x′, y, z′) tilted by a small angle α with regard to (x, y, z), the resampled envelope u′ is obtained by ψ = u e^−ikz = u′ e^−ikz′ = u′ e^{−ik(x sin(α)+z cos(α))} from which it follows that u′ ≈ u e^ikx ^sin(α). Additionally, for such small angles the envelope is approximately identical in both coordinate systems u′(x′) ≈ u′(x). The new envelope in the tilted coordinate system is therefore approximately obtained by multiplying exp (ikx sin(α)) with the original envelope.

To verify the piecewise paraxial propagation scheme, results by the finite difference 2D propagator in a curved slab waveguide are compared in cartesian and curved coordinates. The waveguide consists of a curved vacuum guiding layer with 200 nm diameter, 40 mm radius of curvature and a tantalum cladding at room temperature. The incident field is a plane wave with 7.9 keV photon energy. Figure 3 shows the comparison of fields computed by finite difference propagation in cartesian and curved coordinates.

Fig. 3 Above: normalized intensity distribution of a 7.9 keV plane wave in a tantalum curved slab waveguide with 200 nm diameter and a radius of curvature of 40 mm. The waveguide edges are illustrated by dashed white lines. The step sizes are Δx = 350 pm and Δz = 50 nm. Center: intensity distribution from above simulation, straightened by a coordinate system that curves with the waveguide. Below: the intensity distribution of the curved coordinate system obtained directly by piecewise paraxial propagation. Due to the straightened geometry the size of the simulation box can be chosen much smaller with identical computational effort, resulting in a step size of Δx′ = 40 pm and thus improved accuracy.

Download Full Size | PDF

4. Time-dependent propagation

Next, we discuss propagation of non-stationary, i.e. time-dependent fields in matter, extending the implementation described above to the time domain by spectral decomposition.

4.1. Framework of the numerical approach

For time-dependent problems, we use the time-dependent real-valued field amplitude A $A : = ℱ_{ω}^{- 1} [ψ] (t)$ defined by it its spectrum ψ(ω), as well as the corresponding complex-valued field $A^{'} : = ℱ_{ω}^{- 1} [ψ 1_{{ω \geq 0}}] (t)$ obtained by restricting the Fourier integral to positive frequencies, using the indicator function 1_{ω≥0}, as explained in App. C. The (measurable) real-valued field A′(t) is then computed from A′(t) by A = 2 ℜ(A′). Further, for narrow-banded high-energy fields oscillating around a central frequency ω₀ and corresponding wavevector k₀ = ω₀/c, we consider the time-dependent envelope of the complex field u_ω₀ (t) = A′/exp(i(ω₀t − k₀z)). The time-dependent envelope is related to the stationary envelope u(ω) = ψ/e^−ik₀z by (App. C)

u_{ω_{0}} (t) = ℱ_{ω}^{- 1} [u (ω + ω_{0}) e^{- i k z} 1_{{ω \geq - ω_{0}}}] (t) .

The instantaneous intensity can then be expressed by the complex envelope as

\begin{array}{l} I_{inst} & \approx c ∊_{0} A^{2} = 2 c ∊_{0} ℜ {(u_{ω_{0}} (t) e^{i (ω_{0} t - k_{0} z)})}^{2} \\ = 2 c ∊_{0} {(ℜ (u_{ω_{0}} (t)) cos (ω_{0} t - k_{0} z) - ℑ (u_{ω_{0}} (t)) sin (ω_{0} t - k_{0} z))}^{2} . \end{array}

For x-ray photon energies, the envelope u is much smoother in z than the carrier frequency

| \frac{\partial u_{ω_{0}}}{\partial z} | ≪ | k_{0} u_{ω_{0}} |

, so that the instantaneous intensity Ī_inst can be locally averaged in z over many oscillations. Equivalently the argument can be averaged over many temporal periods, again since the envelope is essentially constant over one period, yielding the locally averaged time-dependent intensity

I_{inst, avg} \approx c ∊_{0} {| u_{ω_{0}} |}^{2} .

We see that the squared magnitude of the time-dependent envelope is proportional to the locally averaged instantaneous intensity.

When numerically solving the wave equation for very short pulses and also for plotting the results, we found it useful to ‘stretch’ the solution in z direction, as the pulse length would otherwise be too small to sample accurately. This can be achieved by a coordinate transform t → t′ = t − az/c, where a ∈ [0, 1] is the stretching parameter. When the beam is described by an envelope $f (t - \frac{z}{c})$ with a spacial length of σ_z, the coordinate transform results in beam with the stretched length σ′_z = σ_z/(1 − a). For a = 0, the beam’s length remains unchanged while for a = 1 the beam’s length is stretched to infinity. It is convenient to introduce the stretching factor s:= σ′_z/σ_z = 1/(1 − a). We thus introduce the stretched complex envelope as a function of t′

\begin{array}{l} u_{ω_{0}} (t^{'}) & = ℱ_{ω}^{- 1} [u (ω + ω_{0}) e^{- i k z} 1_{{ω \geq - ω_{0}}}] (t^{'} + a z / c) \\ = ℱ_{ω}^{- 1} [u (ω + ω_{0}) \cdot e^{- i k z / s} 1_{{ω \geq - ω_{0}}}] (t^{'}) . \end{array}

Note that for fixed z the stretched envelope u_ω₀ (t′) is identical to the actual envelope u_ω₀(t) with a temporal offset. We use the envelope stretching technique to accurately sample the a pulse in z-direction, following the propagation of a single pulse through the discretized simulation box. The group velocity of an ultra-short pulse, for example, is determined from

v_{g} = \frac{Δ z}{Δ t} = \frac{Δ z}{Δ t^{'} + a \frac{Δ z}{c}} = \frac{c}{a + c / v_{g}^{'}},

based on the equivalent group verlocity v′_g = Δz/Δt′ in the stretched coordinate system.

In numerical simulation, the spectrum is discretized with suitable sampling, and correspondingly the signal becomes time-periodic. Hence also the initial signal u_ω₀|_z=0, which is entered in the propagator, is periodic, but with a time interval large enough to prevent temporal interference of the pulse with itself. For the numerical implementation, we then use the following algorithm to determine the stretched time-dependent envelope from the initial signal u_ω₀|_z=0:

Determine the time-dependent envelope of the scalar field u(ω) = ℱ_t[u_ω₀|_z=0] (ω − ω₀)
Use a suitable propagator to calculate the spatially discretized envelope u( ${\vec{x}}_{e f}^{g}$ , ω_h) for a set of N_ω equidistant frequencies ω_h centered around ω₀, where e, f, g, h are the according discretization indices of x, y, z, ω.
Choose a suitable stretching factor s to determine the stretched time-dependent envelope $u_{ω_{0}} ({\vec{x}}_{e f}^{g}, {t^{'}}_{h}) = {DFT}_{h}^{- 1} [u_{e f}^{g} (ω_{h}) e^{- i z ω_{h} / c s}] ({t^{'}}_{h})$ .

The fast Fourier transform is used to efficiently calculate the discrete Fourier transform. The time-dependent finite difference propagator then has a total numerical complexity of 𝒪(N_ωN_xN_yN_z), with typically N_ω ≪ N_z.

4.2. Validation of numerical propagation

The implementation of the periodic envelope propagation algorithm into the PyPropagate framework was validated by its consistency with analytically solvable problems and by tests against a simple test-case, see Fig. 4. Two Gaussian pulses start from the same plane at z = 0 and t = 0 with wavevectors inclinded by an angular offset of α with respect to each other. The paths intersect at z = l₁, where the pulses have propagated for distances l₁ and l₂ = l₁/cos(α), respectively, corresponding to a delay of Δt = l₁ (cos(α)⁻¹ − 1)/c. For pulse durations larger than Δt the pulses interfere, while they will miss each other for shorter durations. The (2D) Gaussian pulses are chosen with central photon energy E₀ = 12 keV, a lateral width of σ_x = 8.5 μm, and a pulse duration of FWHM_t = 0.3 fs. Propagation is computed with discretization step sizes Δx = 0.8 nm, Δz = 3.6 μm and Δω = 0.07 FWHM_ω ≈ 1.3/fs. Given the geometric parameters α = 20 mrad and l₁ = 4 FWHM_tc/(cos(α)⁻¹ − 1) ≈ 1.8 mm and free space propagation, the time difference between the pulses at the intersection point is Δt = 4 FWHM_t = 1.2 fs, so that the pulses cannot interfere. The length of the straight pulse in z direction is FWHM_z = c FWHM_t ≈ 90 nm. The angled beam reaches the intersection point z = l₁ at t = τ:= l₂/c ≈ 60 ps = 200000 FWHM_t. In order to accurately track the propagation of the pulses temporally and spatially we chose a stretching factor s = 2000, which elongates the pulse length to approximately FWHM′_z ≈ 180 μm and changes the moment of intersection to τ′ = τ − a l₁/c ≈ 4.2 fs. In Fig. 4(b) the locally averaged instantaneous intensity is shown at z = l₁, confirming the expected temporal separation of the two beams. The fitted delay between the pulses is Δt = 1.21 (2) fs and the stretched pulse width is FWHM′_z = 180(9) μm, consistent with design and prediction.

Fig. 4 (a) Schematic illustration of the angled pulse test-case. Two Gaussian pulses start at z = 0, one propagating downward at an angle α, and reaching the intersection point at a corresponding time delay. (b) Normalized locally averaged instantaneous intensity of the pulsed Gaussian beams at z = l₁ as a function of x and t′, confirming the expected delay. The dashed line indicates the expected intersection time t′ = τ′. (c,d) Time-averaged normalized intensity around the intersection point, for (c) two monochromatic Gaussian beams, and (d) the two periodic Gaussian pulses with pulse duration of 0.3 fs. Due to the different path lengths, the pulses in (d) do not interfere, as the angled pulse reaches the center 1.2 fs after the straight pulse.

Download Full Size | PDF

Next, the algorithm is validated by determining the dispersion of an ultrashort pulse propagating in water, and comparison with the analytical result [21]. The initial beam is a (1D) Gaussian pulse with central energy E₀ = 12 keV and a duration of FWHM_t = 10 as. As the analytical solution is only valid when the second derivative of n is negligible, we approximate n through its Taylor series up to first order around the central energy, as shown in Fig. 5. Using 2D Fresnel propagation the time dependent field is computed with step sizes Δz = 100 μm and Δω = 0.06 FWHM_ω, and the pulse duration is subsequently determined for all values of z by least-square fitting of Gaussians. The comparison to the analytical solution [21] in Fig. 5 shows a maximum deviation of about 1%.

Fig. 5 Left: The refractive index of water n is approximated by n′, the Taylor series of first order around E₀, to match the analytical approximation. Right: Ratio of the pulse duration Σ_t and the initial pulse duration σ_t as a function of propagation distance for a Gaussian pulse in water at E₀ = 12 keV base photon energy and FWHM_t = 10 as. Analytical values calculated as in [21].

Download Full Size | PDF

4.3. Pulse propagation in waveguides

Just like waveguides in other spectral ranges, x-ray waveguides support a finite number of guided modes, each with its own propagation constant and effective absorption index. Therefore, the modes propagate with different effective dispersion and group velocity, depending on the derivatives of the effective refractive indices. Since the derivatives are extremely small, however, we expect very small effects which shall be quantified by the numerical methods presented here, as an illustration and first application of time-dependent finite difference propagation. To this end, the time dependent field was computed for Gaussian pulses propagating in slab and cylindrical waveguides with 100 nm diameter, composed of silicon cladding and a vacuum (Va) core. In order to conveniently separate different modes by their group velocity despite the small expected dispersion effects, we chose extremely short pulse duration of FWHM_t = 5 as. The spectrum of the pulse and the refractive index of silicon is shown in Fig. 6.

Fig. 6 Left axis: real and imaginary part of the refractive index of silicon. Right axis: Spectrum of the Gaussian pulse with base energy E₀ = 12 kev and FWHM_t = 5 as pulse duration.

Download Full Size | PDF

Figure 7(a)–Fig. 7(c) shows the instantaneous intensity of the beam at different propagation distances, for the case of a slab (planar) waveguide, computed with the 2D propagator at step sizes Δx = 0.2 nm, Δz = 0.6 μm and Δω = 0.07 FWHM_ω. The three guided modes are observed to separate spatially and temporally during propagating. As can be seen, however, the pulse separation is only a few nm after several mm of propagation distances. After switching to stretched coordinates as seen in 7 (d,e), single wave packets can be accurately tracked as they propagate through the waveguide. By fitting Gaussian distributions to the averaged intensities at z = 5 mm, the exact group velocities and dispersion of the modes was determined, and tabulated in Tab. 2.

Fig. 7 (a,b,c) Normalized instantaneous intensity of the 5 attosecond beam for different propagation distances z ≈ 0.5 mm (a), z ≈ 2 mm (b), z ≈ 4 mm (c) in the silicon slab waveguide of 100 nm diameter. The waveguide’s edges are indicated by black lines. (d,e) Using the stretched time coordinate t′, we can track individual pulses as they propagate through the whole system, shown in locally averaged intensity for different time points t′ (d) and again averaged over the waveguide cross-section for all time-points t′ (e) to accurately determine the group velocity. The time-points of (d) are indicated in (e) by striped black lines. The separation of the wave packets into groups traveling with different velocities is clearly visible.

Download Full Size | PDF

Table 2. Group velocity and pulse duration for the first three symmetrical guided modes in the 100 nm diameter silicon waveguide for a Gaussian pulse with FWHM_t = 5 as. The modes denoted as (slab) and (cyl) correspond to slab and cylindrical geometry, respectively.

View Table | View all tables in this article

Analogous simulations were carried out for cylindrical waveguides with 100 nm diameter (data not shown), yielding slightly larger dispersion effects, see the lowest three lines in Tab. 2, owing to the higher percentage of intensity penetrating into the cladding.

Much stronger dispersion effects can be observed when the pulse spectrum is located at an absorption edge of the cladding material as shown in Fig. 8 for the nickel K absorption edge at E ≈ 8.33 keV. The steep increase in the magnitude of the imaginary part of the refractive index, can result in a substantial suppression of the high energy tail of the pulse spectrum. To illustration, a Gaussian pulse with a base energy of E₀ ≈ 8.33 keV and FWHM_t = 300 as was propagated in different nickel structures. The pulse’s duration was found to initially increase before saturating at about 1.82 times of the initial duration, reflecting complete absorption of the high energy tail (beam softening). The dispersion follows similar saturation values and lineshapes with propagation, but on entirely different length scales. In fact, Fig. 8 shows that compared to bulk nickel, the propagation length in the slab waveguide can be thousandfold larger before a reaching comparable pulse broadening.

Fig. 8 (a) Left axis: real and imaginary part of the refractive index of nickel around the K absorption edge. The data is taken from xraylib [22]. Right axis: normalized initial spectrum of the 300 as pulse with a base energy of E₀ ≈ 8.33 keV. (b) Normalized spectrum and (c) locally centered and averaged intensity (right) of an 300 as pulse with a base energy of E₀ ≈ 8.33 keV propagating in nickel at the propagation distances z = 0 and z = 42 mm. (d–f) Normalized pulse width Σ_t(z) and normalized maximal intensity A(z) of 300 as pulses with a base energy of E₀ = 8.3328 keV in homogeneous nickel (d), a 100 nm Va slab waveguide (e) and cylindrical waveguide (f) in a nickel cladding. The pulse widths and maximum intensities are determined by fitting Voigt profiles to the locally averaged intensity profiles. The sudden rise in pulse intensity at the waveguide entrance is due to propagation of the initial plane wave from outside into the waveguide.

Download Full Size | PDF

5. Summary and conclusion

In summary, we have extended finite difference calculation for x-ray propagation from stationary to time-dependent problems. While this is up to now restricted to quasi-monochromatic and linear propagation within scalar theory, the assumptions made hold for a wide range of x-ray propagation problems. Extensions to non-linear effects of high brilliance FEL pulses can be addressed in future work. The conceptual simplicity of the problem treated here, is however contrasted with a tremendous numerical challenge, in view of the demanding sampling constraints for the hard x-ray spectral range. The use of suitable envelope Eqs. is therefore as critical as the FD scheme and numeric implementation. We have first significantly extended the performance of FD propagators and have adapted them to a wider range of dimensionality and symmetry. After a rigorous derivation of scalar propagation theory we designed a single consistent and high-performance numerical framework for FD and Fresnel propagation and adapted them to systems with cylindrical symmetry or curved paraxial beam trajectories. This has then allowed us to treat time-dependent problems based on a suitable sampling of the pulse spectrum. Validation of the numerical implementation by comparison to analytical theory has been a primary concern of this work, since we convinced ourselves how numerical errors can easily propagate and result in flawed field distributions, in particular when considering subtle effect on logarithmic scales [23]. All of the code used is made publicly available as supplemental material.

The numerical implementation and validation then allowed us to explore first applications of ultra-fast pulse propagation in x-ray waveguides. The results showed, that as expected [24], x-ray waveguide support nearly dispersion-free pulse propagation down to ultra-short pulse width in the range of 0.1 fs. This makes them interesting optical devices for the emerging field of ultra-fast x-ray optics at FEL and HHG sources. Only for pulse width of 5 attoseconds we observe the dispersion of different pulses, which now separate spatially as the propagation distance along the optical axis is increased. Compared to the off-resonance case (far away from absorption edges), the dispersion effect is much more pronounced when the pulse spectrum covers an absorption edge of the cladding material. In this case modal dispersion is observable also for a pulse width of 0.3 fs, which is close to be reached experimentally. These results illustrate that dispersion effects for x-ray pulses propagating in x-ray waveguides are extremely small, when waveguides are fabricated by lithographic methods resulting in open air/vacuum channels with a metal or semiconductor cladding [25].

Beyond the presented waveguide application, the numerical approach presented can also be used to treat dispersion and pulse shaping by dispersive optical elements such as gratings or zone plates. To this end, numerical simulation of FEL optics [13] and in particular propagation of ultra-short FEL pulses offers an enhanced tool to design time-resolved imaging as well as x-ray quantum optical experiments [26], by providing an precise understanding how the pulses propagate in matter. We have recently found, for example that special resonant beam couplers can be designed so that two equally strong reflected beams are observed, one displaced along the surface with respect to the other, a kind of giant Goos-Hänchen effect [12]. While that work was limited to the stationary case (both simulation and experiment), one can exploit such phenomena for control of ultra-short and focused x-ray pulses. For example, such multi-guide resonant beam couplers could serve as time-delay FEL or HHG beam splitters with attosecond delay, orders of magnitude smaller than current macroscopic pulse delay stages [3,27]. This will require precise numerical simulation of pulse propagation as provided by the current work.

A. Scalar diffraction theory and paraxial beams

In solving optical diffraction problems one typically takes recourse to a single complex scalar field called the light disturbance [28,29], which obeys the scalar Helmholtz Eq. (11) and whose squared magnitude corresponds to the optical intensity. This is a rather simplified approach compared to determining the full solution of the vectorial Maxwell system. The approach is usually justified by the assumption that the components of the electric field are uncoupled and can therefore be regarded independently, e.g. [30,31]. While the scalar approximation is certainly very well justified for forward-directed x-ray diffraction, it is however, strictly speaking, in direct violation of the Maxwell Eq. (9), which requires the divergence of the electric field to vanish. A rigorous derivation by H. S. Green and E. Wolf shows that real-valued divergence-free fields can indeed be fully described by a complex scalar potential and that under certain circumstances the optical intensity corresponds to the squared magnitude of such a potential [32,33]. However the applicability is limited to free-space (or absorption-free) propagation, as only then the frequency-domain field vectors can be real-valued everywhere. In general, this is not warranted due to the complex nature of the refractive index n in the Helmholtz equation Alternatively, one can use the fact that an arbitrary vector field may be represented by two scalar Debye potentials [34], which is often used for Mie scattering and in multipole expansion problems. However, the construction of the electromagnetic fields from the potentials is quite exhaustive and the direct connection to the optical intensity remains unclear. In view of the x-ray diffraction problems regarded in this work, we can safely assume that scalar theory is an excellent approximation. However, we still need a rigorous derivation how to compute electric and magnetic fields (which obey Maxwell’s Eqs.) from the scalar propagation quantity ψ, which is denoted as the ‘scalar field’ or the scalar potential. To this end, we consider the time-domain Fourier transform of the electric vector field E⃗ and the magnetic vector field B⃗

\hat{\vec{E}} = ℱ_{t} [\vec{E} (ω)] = \frac{1}{\sqrt{2 π}} \int_{ℝ} \vec{E} e^{- i ω t} d t .

When the propagation medium is homogeneous, isotropic, non-magnetic, and does not contain free charges and the field intensity is sufficiently low to disregard non-linear response of polarisation, the well-known Helmholtz equation for the electric field can be derived from Maxwell’s Eqs. for stationary fields

Δ \hat{\vec{E}} + k^{2} n^{2} \hat{\vec{E}} = 0,

where k = ω/c is the wave number, c is the speed of light and n is the complex refractive index of the propagation medium. Under these conditions, the dynamics of the electromagnetic field are fully characterized by the Helmholtz equation and the Maxwell Eqs.

\vec{\nabla} \cdot \hat{\vec{E}} = 0 and \hat{\vec{B}} = \frac{1}{i ω} \vec{\nabla} \times \hat{\vec{E}} .

A divergence-free solution

\hat{\vec{E}}

of the Helmholtz Eq. (8) can be constructed from any (not necessarily divergence-free) solution ψ⃗ [V/m] by taking the curl of ψ⃗

\hat{\vec{E}} = \frac{i}{k} \vec{\nabla} \times \vec{ψ} .

This directly follows from the vector calculus identities

\begin{array}{l} \vec{\nabla} \cdot \hat{\vec{E}} = \frac{i}{k} \vec{\nabla} \cdot (\vec{\nabla} \times \vec{ψ}) = 0, \\ Δ \hat{\vec{E}} = \frac{i}{k} Δ (\vec{\nabla} \times \vec{ψ}) = \frac{i}{k} \vec{\nabla} \times (Δ \vec{ψ}) = - n^{2} k^{2} \hat{\vec{E}} . \end{array}

Note that the prefactor i/k is chosen to simplify subsequent expressions. For any vector e⃗, the projection (ψ⃗·e⃗) e⃗ is also a solution of the Helmholtz equation In particular, given a solution ψ of the scalar Helmholtz equation

Δ ψ + k^{2} n^{2} ψ = 0,

we can construct solutions of the Maxwell system Eq. (8) and Eq. (9) by setting ψ⃗ = ψ⃗e⃗_p for any unit vector e⃗_p for which

\hat{\vec{E}} = \frac{i}{k} (\vec{\nabla} ψ) \times {\vec{e}}_{p} and \hat{\vec{B}} = \frac{1}{c} (\frac{1}{k^{2}} \vec{\nabla} ({\vec{e}}_{p} \cdot \vec{\nabla} ψ) + n^{2} ψ {\vec{e}}_{p}) .

We will denote ψ the scalar potential of the electromagnetic field, but in scalar diffraction theory it is mostly simply denoted as the ‘scalar field’. General solutions can be constructed from linear combinations of three independent scalar potentials ψ_x, ψ_y, ψ_z so that ψ⃗ = ψ_x e⃗_x + ψ_ye⃗_y + ψ_ze⃗_z. The electromagnetic field vectors are then given by the sum of the fields determined individually from the potentials.

For the special case of paraxial beams, solutions of the scalar Helmholtz equation are locally approximated by the product of a slowly varying envelope u and a fast oscillating term exp (−ink⃗·x⃗), i.e.

ψ = u e^{- i n \vec{k} \cdot \vec{x}},

and the wave vector k⃗ satisfying the condition |k⃗| = k. If we choose k⃗ ‖ e⃗_z and e⃗_p = e⃗_y, the electromagnetic field vectors are well approximated by

\hat{\vec{E}} \approx n ψ {\vec{e}}_{x}

and

\hat{\vec{B}} \approx \frac{n^{2}}{c} ψ {\vec{e}}_{y}

.

B. Time averaged energy flux

To obtain the correct prefactors and for notational clarity, we briefly derive the time-averaged energy flux for the numerically computed scalar field, starting from the full electromagnetic field and the Poynting vector

\vec{S} = μ_{0}^{- 1} \vec{E} \times \vec{B},

where μ₀ = 4π 10⁻⁷N/A² is the free space permeability. The energy flux is indirectly measurable by detectors as the integral of the energy flux over the detector surface and the given exposure time. Thus, one considers the energy flux averaged (or integrated) over the exposure time, which is simply referred to as the optical intensity. When the signal is time-periodic, we can write the field as as a discrete sum of frequencies ω_i for i ∈ ℕ

\begin{array}{l} \vec{E} & = ℱ_{ω}^{- 1} [\hat{\vec{E}} (t)] \\ = \sqrt{\frac{2}{π}} \sum_{i \in ℕ} ℜ (\hat{\vec{E}} (ω_{i})) cos (ω_{i}) + ℑ (\hat{\vec{E}} (ω_{i})) sin (ω_{i}), \end{array}

and the Poynting vector equation 14 can be written as a double sum containing cross-products of fields and products of trigonometric functions. In the time-average

〈 f 〉 = {lim}_{T \to \infty} T^{- 1} \int_{0}^{T} f d t

these trigonometric products reduce to either 〈cos(ω_it) sin(ω_jt)〉 = 0 or 〈cos(ω_it) sin(ω_jt)〉 = δ_ij/2, where δ_ij denotes the Kronecker delta symbol. The double sum is therefore reduced to a single sum

〈 \vec{S} 〉 = \frac{1}{μ_{0} π} \sum_{i \in ℕ} ℜ (\hat{\vec{E}} (ω_{i})) \times ℜ (\hat{\vec{B}} (ω_{i})) + ℑ (\hat{\vec{E}} (ω_{i})) \times ℑ (\hat{\vec{B}} (ω_{i})) .

When the fields are constructed from a single scalar potential ψ which oscillates slowly in e⃗_p direction and therefore fulfills the condition

‖ \vec{\nabla} ({\vec{e}}_{p} \cdot \vec{\nabla} ψ) ‖ ≪ | n^{2} k^{2} ψ |,

we can approximate

\begin{array}{l} ℜ (\hat{\vec{E}}) \times ℜ (\hat{\vec{B}}) & = ℜ (\frac{i}{k} \vec{\nabla} ψ \times {\vec{e}}_{p}) \times ℜ (\frac{1}{c} (\frac{1}{k^{2}} \vec{\nabla} ({\vec{e}}_{p} \cdot \vec{\nabla} ψ) + n^{2} ψ {\vec{e}}_{p})) \\ = - \frac{1}{c k} ℑ (\vec{\nabla} ψ + {\vec{e}}_{p}) \times ℜ (\frac{1}{k^{2}} \vec{\nabla} ({\vec{e}}_{p} \cdot \vec{\nabla} ψ) + n^{2} ψ {\vec{e}}_{p}) \\ \approx - \frac{1}{ω} ℑ (\vec{\nabla} ψ \times {\vec{e}}_{p}) \times ℜ (n^{2} ψ {\vec{e}}_{p}) \\ \approx - \frac{ℜ (n^{2})}{ω} ℜ (ψ) ℑ ((\vec{\nabla} ψ \times {\vec{e}}_{p}) \times {\vec{e}}_{p}) \\ = \frac{ℜ (n^{2})}{ω} ℜ (ψ) ℑ (\vec{\nabla} ψ) \end{array}

and, analogously

ℑ (\hat{\vec{E}}) \times ℑ (\hat{\vec{B}}) \approx - \frac{ℜ {(n)}^{2}}{ω} ℑ (ψ) ℜ (\vec{\nabla} ψ),

where we used the fact that for x-rays we usually have ℜ(n) ≫ = ℑ(n). The condition Eq. (15) is valid for two dimensional beams which are constant in e⃗_p direction and for paraxial beams when e⃗_p is chosen perpendicular to the direction of propagation. By a straightforward calculation it follows that

ℜ (ψ) ℑ (\vec{\nabla} ψ) - ℑ (ψ) ℜ (\vec{\nabla} ψ) = {| ψ |}^{2} \vec{\nabla} arg (ψ) .

Hence, we can approximate the time-averaged Poynting vector for periodic signals with negligible derivative in e⃗_p direction as

〈 \vec{S} 〉 \approx \frac{c ∊_{0}}{π} \sum_{i \in ℕ} ℜ {(n (ω_{i}))}^{2} {| ψ (ω_{i}) |}^{2} \frac{\vec{\nabla} arg (ψ (ω_{i}))}{k_{i}},

where k_i = ω_i/c. For paraxial beams as previously defined the modulus of the phase gradient is approximated by ‖∇⃗ arg(ψ(ω_i))‖ ≈ ℜ((n(ω_i)) k_i, and hence the time-averaged Poynting vector becomes

〈 \vec{S} 〉 \approx \frac{c ∊_{0}}{π} \sum_{i \in ℕ} ℜ {(n (ω_{i}))}^{3} {| ψ (ω_{i}) |}^{2} \frac{\vec{\nabla} arg (ψ (ω_{i}))}{‖ \vec{\nabla} arg (ψ (ω_{i})) ‖} .

For a monochromatic paraxial beam, the time-averaged Poynting vector has hence a magnitude, or optical intensity I, proportional to the squared scalar field (denoted as scalar potential in appendix A)

I = ‖ 〈 \vec{S} 〉 ‖ \approx \frac{c ∊_{0}}{π} ℜ {(n)}^{3} {| ψ |}^{2}

and is oriented in the direction of the phase gradient.

C. Time-dependent fields: definitions and notation

For narrow-banded paraxial beams in the linear optical regime the the time-dependent energy density in free space, or instantaneous intensity I_inst, is approximated by I_inst = c∊₀ A², where A is the time-dependent field amplitude. We define the time-dependent field amplitude by its spectrum ψ(ω)

A : = ℱ_{ω}^{- 1} [ψ] (t),

while A′ denotes the spectrum restricted to positive frequencies only

\begin{array}{l} A^{'} & : = ℱ_{ω}^{- 1} [ψ 1_{{ω \geq 0}}] (t) \\ = \frac{1}{\sqrt{2 π}} \int_{0}^{\infty} ψ e^{i ω t} d ω, \end{array}

where 1_{ω≥0} is the indicator function of the set {ω ≥ 0} and is defined as 1_{ω≥a} = 1 if ω ≥ a and 1_{ω≥a} = 0 otherwise, for any a ∈ ℝ. Given the Hermitian symmetry of ψ(ω), A(t) is real-valued, while A′(t) is a complex-valued field amplitude, calculated from the positive spectrum. Note that by use of the complex field amplitude A′(t), calculations of time dependent propagation are considerably simplified. The actual real-valued field A(t) is then computed from A′(t) by

A = ℱ_{ω}^{- 1} [ψ (t)] = 2 ℜ (A^{'}) .

As narrow-banded high-energy fields oscillate at very high frequencies (we define the central frequencies as ω₀ in time and k₀ = ω₀/c in space domain), we consider the time-dependent envelope u_ω₀(t) = A′/exp(i(ω₀t − k₀z)) of the complex field amplitude. The relationship of the time-dependent envelope with the stationary envelope u(ω) = ψ/e^−ik₀z, as used by the stationary propagation schemes mentioned before, is given by

u_{ω_{0}} (t) = ℱ_{ω}^{- 1} [u (ω + ω_{0}) e^{- i k z} 1_{{ω \geq - ω_{0}}}] (t),

which shows that the time-dependent envelope may be determined by the Fourier transform of the positive-frequency stationary envelopes, centered around the base frequency ω₀.

Funding

Deutsche Forschungsgemeinschaft (DFG/SFB 755); Bundesministerium für Bildung und Forschung (BMBF) (05K16MGB).

Acknowledgments

We thank Markus Osterhoff and Jan Goeman for discussion and advice in computing.

References and links

1. D. M. Paganin, Coherent X-ray Optics (Oxford University, 2006). [CrossRef]

2. H.N. Chapman and K.A. Nugent, “Coherent lensless x-ray imaging,” Nat. Photon 4(12), 833–839 (2010). [CrossRef]

3. A. Barty, “Time-resolved imaging using x-ray free electron lasers,” J. Phys. B: At. Mol. Opt. Phys. 43(19), 194014 (2010). [CrossRef]

4. J. Crank and P. Nicolson, “A practical method for numerical evaluation of solutions of partial differential Eqs. of the heat-conduction type,” Math. Proc. Cam. Phil. Soc. 43(1), 50–67 (1947). [CrossRef]

5. W. H. Press, S.A. Teukolsky, W.T. Vetterling, and B. P. Flannery, Numerical Recipes, 3rd ed. (Cambridge University, 2007).

6. Y. V. Kopylov, A. V. Popov, and A. V. Vinogradov, “Application of the parabolic wave equation to X-ray diffraction optics,” Opt. Commun. 118, 619–636 (1995). [CrossRef]

7. J. H. H. Bongaerts, C. David, M. Drakopoulos, M. J. Zwanenburg, G. H. Wegdam, T. Lackner, H. Keymeulen, and J. F. van der Veen, “Propagation of a partially coherent focused X-ray beam within a planar X-ray waveguide,” J. Synchrotron Rad. 9, 383–393 (2002). [CrossRef]

8. C. Bergemann, H. Keymeulen, and J. F. van der Veen, “Focusing X-Ray Beams to Nanometer Dimensions,” Phys. Rev. Lett. 91, 204801 (2003). [CrossRef] [PubMed]

9. C. Fuhse and T. Salditt, “Finite-difference field calculations for one-dimensionally confined X-ray waveguides,” Phys. B: Condensed Matter 357(1–2), 57–60 (2005). [CrossRef]

10. C. Fuhse and Tim Salditt, “Finite-difference field calculations for two-dimensionally confined x-ray waveguides,” Appl. Opt. 45(19), 4603–4608 (2006). [CrossRef] [PubMed]

11. S. P. Krüger, K. Giewekemeyer, S. Kalbfleisch, M. Bartels, H. Neubauer, and T. Salditt, “Sub-15 nm beam confinement by twocrossed x-ray waveguides,” Opt. Express 18(13), 13492–13501 (2010). [CrossRef] [PubMed]

12. Q. Zhong, M. Osterhoff, M. W. Wen, Z. S. Wang, and T. Salditt, “X-ray waveguide arrays: tailored near fields by multi-beam interference,” X-Ray Spectrom. 46(2), 107–115 (2017). [CrossRef]

13. L. Samoylova, A. Buzmakov, O. Chubar O, and H. Sinn , “WavePropaGator: interactive framework for X-ray free-electron laser optics design and simulations,” J. Appl. Cryst. 49(4), 1347–1355 (2016). [CrossRef]

14. I. Barke, H. Hartmann, D. Rupp, L. Flückiger, M. Sauppe, M. Adolph, S. Schorb, C. Bostedt, R. Treusch, C. Peltz, S. Bartling, T. Fennel, K. Meiwes-Broer, and T. Möller, “The 3D-architecture of individual free silver nanoparticles captured by X-ray scattering,” Nat. Commun. 6, 6187 (2015). [CrossRef] [PubMed]

15. L. Li, M. Wojcik, and C. Jacobsen, “Multislice does it all – calculating the performance of nanofocusing X-ray optics,” Opt. Express 25(3), 13107–13124 (2017).

16. L. Melchior and T. Salditt, “Finite difference methods for stationary and time-dependent x-ray propagation - Simulations and Figures,” Figshare (2017) [retrieved 16 November 2017], https://doi.org/10.6084/m9.figshare.5449939.

17. L. Melchior, “PyPropagate on Github,” Github (2017) [retrieved 27 September 2017], https://github.com/TheLartians/PyPropagate.

18. A. Husakou, “Nichtlineare Phaenomene spektral ultrabreiter Strahlung in Photonischen Kristallfasern und Hohlen Wellenleitern,” Ph.D. dissertation, Freie Universität Berlin, Germany (2002).

19. F. Pérez and B.E. Granger, “IPython: a System for Interactive Scientific Computing,” Computing in Science and Engineering 9(3), 21–29 (2007). [CrossRef]

20. A.W. Norfolk and E.J. Grace, “Reconstruction of optical fields with the Quasi-discrete Hankel transform,” Opt. Express 18(10), 10551–10556 (2010). [CrossRef] [PubMed]

21. B.E.A Saleh and M.C. Teich, “Fundamentals of photonics,” Wiley (1991). [CrossRef]

22. T. Schoonjans, A. Brunetti, B. Golosio, M. Sanchez del Rio, A. Sole, C. Ferrero, and L. Vincze, “xraylib 3.1.0,” xraylib 3.1.0. Zenodo https://doi.org/10.5281/zenodo.12378 (2014).

23. T. Salditt, S. Hoffmann, M. Vassholz, J. Haber, M. Osterhoff, and J. Hilhorst, “X-Ray Optics on a Chip: Guiding X Rays in Curved Channels, Phys. Rev. Lett. 115, 203902 (2015). [CrossRef]

24. A. Jarre, T. Salditt, T. Panzner, U. Pietsch, and F. Pfeiffer, “White beam x-ray waveguide optics,” Appl. Phys. Lett. 85(2),” 161–163 (2004). [CrossRef]

25. S. Hoffmann-Urlaub and T. Salditt, “Miniaturized beamsplitters realized by X-ray waveguides,” Acta Crystallographica Section A , 72(5), 515–522 (2016). [CrossRef]

26. R. Röhlsberger, K. Schlage, B. Sahoo, S. Couet, and R. Rüffer, “Collective Lamb Shift in Single-Photon Superradiance,” Science 328(5983), 1248–1251 (2010). [CrossRef] [PubMed]

27. C. K. Schmising, B. Pfau, M. Schneider, C. M. Günther, M. Giovannella, J. Perron, B. Vodungbo, L. Müller, F. Capotondi, E. Pedersoli, N. Mahne, J. Lüning, and S. Eisebitt, “Imaging Ultrafast Demagnetization Dynamics after a Spatially Localized Optical Excitation,” Phys. Rev. Lett. 112(21), 217203 (2014). [CrossRef]

28. A. S. Marathay and G. B. Parrent, “Use of Scalar Theory in Optics,” J. Opt. Soc. Am. 60(2), 243–245 (1970). [CrossRef]

29. M. Ornigotti and A. Aiello, “The Hertz vector revisited: A simple physical picture,” J. Opt. 16(10), 105705 (2014). [CrossRef]

30. J.W. Goodman, Introduction to Fourier Optics (McGraw-Hill, 1996).

31. J.D. Schmidt, “Numerical Simulation of Optical Wave Propagation with Examples in MATLAB,” SPIE monographs (2010).

32. H.S. Green and E. Wolf, “A Scalar Representation of Electromagnetic Fields,” Proc. Phys. Soc. A 66(12), 1129 (1953). [CrossRef]

33. E. Wolf, “A Scalar Representation of Electromagnetic Fields: II,” Proc. Phys. Soc. A , 74(3), 269 (1959). [CrossRef]

34. C. G. Gray and B.G. Nickel, “Debye potential representation of vector fields,” Am. J. Phys. 47(8), 736 (1979).

Propagator	Coordinates	Edge condition	Algorithm	Time complexity	Runtime

Finite Differences 2D	(x, z)	Dirichlet	CN	𝒪(N_x)	0.3 s
Finite Differences 3D	(x, y, z)	Dirichlet	ADI	𝒪(N_x N_y)	120 s
Finite Differences CS	(r, z)	Dirichlet	CN	𝒪(N_r)	0.3 s
Fresnel 2D	(x, z)	Periodic	FFT	𝒪(N_x log(N_x))	0.6 s
Fresnel 3D	(x, y, z)	Periodic	FFT	𝒪(N_xN_y log(N_xN_y))	260 s
Fresnel CS	(r, z)	Dirichlet (= 0)	DHT	$𝒪 (N_{r}^{2})$	30 s

Mode	v_G/c	Σ_t(z = 5 mm)/Σ_t(z = 0)

slab 1	1 – 0.18 × 10⁻⁶	1
slab 2	1 – 1 × 10⁻⁶	1.02
slab 3	1 – 2.6 × 10⁻⁶	1.12
cyl 1	1 – 0.33 × 10⁻⁶	1
cyl 2	1 – 1.36 × 10⁻⁶	1.04
cyl 3	1 – 3.07 × 10⁻⁶	1.18

Propagator	Coordinates	Edge condition	Algorithm	Time complexity	Runtime

Finite Differences 2D	(x, z)	Dirichlet	CN	𝒪(N_x)	0.3 s
Finite Differences 3D	(x, y, z)	Dirichlet	ADI	𝒪(N_x N_y)	120 s
Finite Differences CS	(r, z)	Dirichlet	CN	𝒪(N_r)	0.3 s
Fresnel 2D	(x, z)	Periodic	FFT	𝒪(N_x log(N_x))	0.6 s
Fresnel 3D	(x, y, z)	Periodic	FFT	𝒪(N_xN_y log(N_xN_y))	260 s
Fresnel CS	(r, z)	Dirichlet (= 0)	DHT	$𝒪 (N_{r}^{2})$	30 s

Mode	v_G/c	Σ_t(z = 5 mm)/Σ_t(z = 0)

slab 1	1 – 0.18 × 10⁻⁶	1
slab 2	1 – 1 × 10⁻⁶	1.02
slab 3	1 – 2.6 × 10⁻⁶	1.12
cyl 1	1 – 0.33 × 10⁻⁶	1
cyl 2	1 – 1.36 × 10⁻⁶	1.04
cyl 3	1 – 3.07 × 10⁻⁶	1.18

Finite difference methods for stationary and time-dependent X-ray propagation

Abstract

1. Introduction

2. The paraxial Helmholtz equation

3. Numerical propagation

3.1. PyPropagate

3.2. Comparison of finite difference and Fresnel propagation

4. Time-dependent propagation

4.1. Framework of the numerical approach

4.2. Validation of numerical propagation

4.3. Pulse propagation in waveguides

5. Summary and conclusion

A. Scalar diffraction theory and paraxial beams

B. Time averaged energy flux

C. Time-dependent fields: definitions and notation

Funding

Acknowledgments

References and links

Supplementary Material (1)

Cited By

Figures (8)

Tables (2)

Equations (35)

Optics Express