## Abstract

We have generalized finite-difference (FD) simulations for time-dependent field propagation problems, in particular in view of ultra-short x-ray pulse propagation and dispersion. To this end, we first derive the stationary paraxial (parabolic) wave equation for the scalar field envelope in a more general manner than typically found in the literature. We then present an efficient FD implementation of propagators for different dimensionality for stationary field propagation, before we treat time-dependent problems by spectral decomposition, and suitable numerical sampling. We prove the validity of the numerical approach by comparison to analytical theory, using simple tractable propagation problems. Finally, we apply the framework to the problem of modal dispersion in X-ray waveguide. We show that X-ray waveguides can be considered as non-dispersive optical elements down to sub-femtosecond pulse width. Only when considering resonant absorption close to an X-ray absorption edge, we observe pronounced dispersion effects for experimentally achievable pulse widths. All code used for the work is made available as supplemental material.

© 2017 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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 $\left|\frac{{\partial}^{2}u}{\partial {z}^{2}}\right|\ll \left|k\frac{\partial u}{\partial z}\right|$ to obtain a parabolic wave equation for the envelope

*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

*ψ̃*=

*ℱ*

_{x⃗⊥}[

*ψ*](

*k⃗*

_{⊥}) and $\beta :=\sqrt{{k}^{2}{n}^{2}-{k}_{\perp}^{2}}$. The differential operator may be written as the commutative product of two operators

*ψ̃*solves either of the differential Eqs.

*k⃗*is oriented in −

*e⃗*direction and the right equation applies to a beam where the wave vector is oriented in

_{z}*e⃗*direction. For paraxial beams, where wave vectors nearly parallel to the

_{z}*z*axis, the support of

*ψ̃*is contained in a small region where ${k}_{\perp}^{2}\ll {k}^{2}$. We can therefore approximate the square root in

*β*through its first-order Taylor series around ${k}_{\perp}^{2}=0$ Substituting this into Eq. (3) (right) and transforming back to real space results in the following paraxial Helmholtz equation 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*), 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.

#### 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*^{2}, 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

*ψ*

_{0}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.

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.

## 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:={\mathcal{F}}_{\omega}^{-1}[\psi ](t)$ defined by it its spectrum *ψ*(*ω*), as well as the corresponding complex-valued field ${A}^{\prime}:={\mathcal{F}}_{\omega}^{-1}\left[\psi \hspace{0.17em}\hspace{0.17em}{\mathbf{1}}_{\left\{\omega \ge 0\right\}}\right](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 *ω*_{0} and corresponding wavevector *k*_{0} = *ω*_{0}/*c*, we consider the time-dependent envelope of the complex field *u*_{ω0} (*t*) = *A′*/exp(*i*(*ω*_{0}*t* − *k*_{0}*z*)). The time-dependent envelope is related to the stationary envelope *u*(*ω*) = *ψ*/*e*^{−ik0z} by (App. C)

*u*is much smoother in

*z*than the carrier frequency $\left|\frac{\partial {u}_{{\omega}_{0}}}{\partial z}\right|\ll \left|{k}_{0}{u}_{{\omega}_{0}}\right|$, 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*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\left(t-\frac{z}{c}\right)$ with a spacial length of *σ _{z}*, the coordinate transform results in beam with the stretched length

*σ′*=

_{z}*σ*/(1 −

_{z}*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}*σ*= 1/(1 −

_{z}*a*). We thus introduce the

*stretched complex envelope*as a function of

*t′*

*z*the stretched envelope

*u*

_{ω0}(

*t′*) is identical to the actual envelope

*u*

_{ω0}(

*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}*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*_{ω0}|_{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*_{ω0}|_{z=0}:

- Determine the time-dependent envelope of the scalar field
*u*(*ω*) =*ℱ*[_{t}*u*_{ω0}|_{z=0}] (*ω*−*ω*_{0}) - Use a suitable propagator to calculate the spatially discretized envelope
*u*(${\overrightarrow{x}}_{ef}^{g}$,*ω*) for a set of_{h}*N*equidistant frequencies_{ω}*ω*centered around_{h}*ω*_{0}, 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}_{{\omega}_{0}}\left({\overrightarrow{x}}_{ef}^{g},{{t}^{\prime}}_{h}\right)={\text{DFT}}_{h}^{-1}\left[{u}_{ef}^{g}({\omega}_{h}){e}^{-iz{\omega}_{h}/cs}\right]\left({{t}^{\prime}}_{h}\right)$.

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_{x}N_{y}N_{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*_{1}, where the pulses have propagated for distances *l*_{1} and *l*_{2} = *l*_{1}/cos(*α*), respectively, corresponding to a delay of Δ*t* = *l*_{1} (cos(*α*)^{−1} − 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*_{0} = 12 keV, a lateral width of *σ _{x}* = 8.5

*μ*m, and a pulse duration of FWHM

*= 0.3 fs. Propagation is computed with discretization step sizes Δ*

_{t}*x*= 0.8 nm, Δ

*z*= 3.6

*μ*m and Δ

*ω*= 0.07 FWHM

*≈ 1.3/fs. Given the geometric parameters*

_{ω}*α*= 20 mrad and

*l*

_{1}= 4 FWHM

*/(cos(*

_{t}c*α*)

^{−1}− 1) ≈ 1.8 mm and free space propagation, the time difference between the pulses at the intersection point is Δ

*t*= 4 FWHM

*= 1.2 fs, so that the pulses cannot interfere. The length of the straight pulse in*

_{t}*z*direction is FWHM

*=*

_{z}*c*FWHM

*≈ 90 nm. The angled beam reaches the intersection point*

_{t}*z*=

*l*

_{1}at

*t*=

*τ*:=

*l*

_{2}/

*c*≈ 60 ps = 200000 FWHM

*. In order to accurately track the propagation of the pulses temporally and spatially we chose a stretching factor*

_{t}*s*= 2000, which elongates the pulse length to approximately FWHM′

*≈ 180*

_{z}*μ*m and changes the moment of intersection to

*τ′*=

*τ*−

*a l*

_{1}/

*c*≈ 4.2 fs. In Fig. 4(b) the locally averaged instantaneous intensity is shown at

*z*=

*l*

_{1}, 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′

*= 180(9)*

_{z}*μ*m, consistent with design and prediction.

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*_{0} = 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%.

#### 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.

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.

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*_{0} ≈ 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.

## 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⃗*

*Helmholtz equation*for the electric field can be derived from Maxwell’s Eqs. for stationary fields 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.

*ψ⃗*[V/m] by taking the curl of

*ψ⃗*This directly follows from the vector calculus identities

*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 we can construct solutions of the Maxwell system Eq. (8) and Eq. (9) by setting

*ψ⃗*=

*ψ⃗e⃗*for any unit vector

_{p}*e⃗*for which

_{p}*ψ*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}*ψ*so that

_{z}*ψ⃗*=

*ψ*

_{x}*e⃗*+

_{x}*ψ*+

_{y}e⃗_{y}*ψ*. The electromagnetic field vectors are then given by the sum of the fields determined individually from the potentials.

_{z}e⃗_{z}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.

*k⃗*satisfying the condition |

*k⃗*| =

*k*. If we choose

*k⃗*‖

*e⃗*and

_{z}*e⃗*=

_{p}*e⃗*, the electromagnetic field vectors are well approximated by $\widehat{\overrightarrow{E}}\approx n\psi {\overrightarrow{e}}_{x}$ and $\widehat{\overrightarrow{B}}\approx \frac{{n}^{2}}{c}\psi {\overrightarrow{e}}_{y}$.

_{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

where*μ*

_{0}= 4

*π*10

^{−7}

*N/A*

^{2}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

*ω*for

_{i}*i*∈ ℕ

*ω*) sin(

_{i}t*ω*)〉 = 0 or 〈cos(

_{j}t*ω*) sin(

_{i}t*ω*)〉 =

_{j}t*δ*/2, where

_{ij}*δ*denotes the Kronecker delta symbol. The double sum is therefore reduced to a single sum

_{ij}*ψ*which oscillates slowly in

*e⃗*direction and therefore fulfills the condition

_{p}*n*) ≫ = ℑ(

*n*). The condition Eq. (15) is valid for two dimensional beams which are constant in

*e⃗*direction and for paraxial beams when

_{p}*e⃗*is chosen perpendicular to the direction of propagation. By a straightforward calculation it follows that

_{p}*e⃗*direction as

_{p}*k*=

_{i}*ω*/

_{i}*c*. For paraxial beams as previously defined the modulus of the phase gradient is approximated by ‖∇⃗ arg(

*ψ*(

*ω*))‖ ≈ ℜ((

_{i}*n*(

*ω*))

_{i}*k*, and hence the time-averaged Poynting vector becomes

_{i}*I*, proportional to the squared scalar field (denoted as scalar potential in appendix A)

## 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∊*_{0} *A*^{2}, where *A* is the *time-dependent field amplitude*. We define the time-dependent field amplitude by its spectrum *ψ*(*ω*)

*A′*denotes the spectrum restricted to positive frequencies only

**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

*ω*

_{0}in time and

*k*

_{0}=

*ω*

_{0}/

*c*in space domain), we consider the time-dependent envelope

*u*

_{ω0}(

*t*) =

*A′*/exp(

*i*(

*ω*

_{0}

*t*−

*k*

_{0}

*z*)) of the complex field amplitude. The relationship of the time-dependent envelope with the stationary envelope

*u*(

*ω*) =

*ψ*/

*e*

^{−ik0z}, as used by the stationary propagation schemes mentioned before, is given by

*ω*

_{0}.

## 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).