Nano-optical theory of planar x-ray waveguides

Leon M. Lohse; Leon M. Lohse; Petar Andrejić

doi:10.1364/OE.504206

1. Introduction

Half a century after the first demonstration of x-ray waveguiding [1], x-ray waveguides have been advanced [2–4] to the point that they can now be routinely used as optical elements for imaging [5,6]. They can bend, split and recombine x-rays on chips [7,8], and they are employed as cavity resonators for quantum-optical experiments [9–12]. And soon, due to the increasing availability of x-ray free electron lasers, new kinds of phenomena that require extremely high intensities facilitated by x-ray waveguides, such as transient nuclear inversion [13], might come into reach.

X-ray waveguides and cavities have enabled a multitude of experiments in different fields. For phase-contrast imaging and tomography, x-ray waveguides provide highly coherent [14–18] nanometer-sized secondary point-sources [5,6,19,20]. These are used to acquire 3d images with resolutions well below 100 nm for research ranging from medicine, biology to chemical catalysis. Recently, Vassholz and Salditt have demonstrated the effective generation of x-rays inside a multilayer waveguide [21] which could lead to more brilliant x-ray sources for the laboratory. For fundamental quantum-mechanics, x-ray waveguides provide one-dimensional confinement and enhanced coupling to atoms embedded into the waveguide, thereby forming a cavity resonator. Exploiting this, Röhlsberger and coworkers embedded Mössbauer nuclei into a waveguide and reported the first observation of the collective Lamb shift [9]. This paved the way for many more demonstrations of quantum-optical phenomena in the hard x-ray regime with Mössbauer nuclei (see the review [12]) and atomic resonances, including on artificial two-level systems and demonstrating spectral control of inner-shell resonances [10,11]. In front coupling, Lee, Ahrens and Liao [22] have proposed using the longitudinal Rabi oscillations between two modes coupled to Mössbauer nuclei as a sensitive measure of the gravitational red-shift, while we have proposed using a similar setup as a platform for implementing geometry-dependent super-radiance in Mössbauer nuclei [23].

Ongoing developments promise to reach new regimes with the help of x-ray waveguides. X-ray free electron lasers are much brighter than synchrotron sources so that multi-photon phenomena with resonant nuclei [12] or non-linear phenomena with resonant inner-shell electrons [13] become feasible. Moreover, beam-splitters and curved waveguides [8,24] combined with resonant nuclei could lead to more complex quantum-optical experiments on millimeter-sized chips.

A sound theoretical foundation facilitates interpretation of results and the design of new experiments. Diffraction of hard x-rays, electromagnetic waves with wavelengths in the order of 10⁻¹⁰ m, largely follow the same principles as light in the visible and infrared wavelengths and can be well described in terms of refractive indices $n$. Theory for dielectric waveguides, mainly focussing on optical wavelengths, has been extensively studied already in the 1960s (see for example [25]). Nowadays, also due to their important applications in telecommunication, dielectric waveguides are well-understood by the optics and electrical engineering communities.

However, x-ray waveguides have important differences to their longer-wavelength counterparts, because the refractive indices, typically written as

(1)$$n = 1 - \delta + i \beta, \quad \delta, \beta {\lesssim} 10^{-5},$$

in the x-ray regime, are generally very close to unity. As a first consequence, the index contrast is tiny so that x-ray waveguides are generally only weakly guiding resulting in the extent of guided modes being many times larger than the wavelength. Second, absorption due to photoionization is appreciably large, $\textrm{Im}\{n\} = \beta > 0$, and the mean free path of photons typically ranges from a mere micrometers in heavy materials up to a few millimeters in lighter materials.

For longer wavelengths, absorption is mostly of little concern because other losses dominate. As a result, the theory of optical waveguides has largely been formulated for real refractive indices and losses are considered separately – notable exeptions are [26–29]. The restriction to real indices has tremendous advantages because the theory can now be formulated in terms of self-adjoint Sturm-Liouville theory, closely related to potential scattering theory in quantum mechanics. For x-ray waveguides, the existing theory of waveguides with real indices has also been the method of choice and absorption has been included ad hoc [14,16,30–33]. In addition, semi-analytical calculations [16] and numerical finite-differences propagation have been performed [34–37] to accurately simulate the performance of waveguides. Notably, semi-analytical calculations based on the full complex refractive indices have also been reported [38]. The design and optimization of parameters for x-ray waveguides has also been studied from the perspective of applied mathematics [39]. In that work, important theorems about the existence of resonances and their dependence on the refractive index profile have been rigorously proven. Several other papers covering specific theoretical aspects of x-ray waveguides have been published, including dispersion in 3-layer waveguides [31] and wave-field formation in front-coupling geometry [40].

Waveguiding is closely linked to a material structure supporting guided modes, so one of the main tasks has been to find those guided modes for given parameters. For simple 3-layer waveguides and negligible material absorption, the modes can readily be written down analytically [25]. This becomes unfeasible for more complicated layer structures and is no longer possible when absorption is included. For optical wavelengths, various analytical approaches and numerical techniques have been developed that can reliably find and characterize the modes both for 1d (planar) and 2d confinement (see for example [41]) for (in principle) arbitrary layer structures. In the x-ray regime, however, the aforementioned weak guiding and finite absorption poses its own challenge and the established techniques can not trivially been applied. Here, the guided modes have been characterized with finite-difference propagation simulations [34,42] and with Numerov’s method [43].

Planar waveguides that are probed in grazing incidence, historically referred to as resonant beam couplers, are much simpler to describe because, assuming that the incoming fields are plane waves, fields and refractive indices are translationally invariant and the problem becomes effectively one-dimensional. The field profiles and reflected intensities for arbitrary complex layer structures can then be readily computed with transfer matrix calculations [44]. Based on this method, a number of software packages have been developed to cater the popularity of x-ray reflectivity measurements (see for example [45]).

When classical or quantum scatterers that are not modeled macroscopically by the complex indices are embedded in a waveguide, the system can no longer be described purely in terms of the fields resulting from the index profile. However, these scattering problems can be formulated [46–48] in terms of the classical electromagnetic Green’s function (see for example [49]), which is fully determined by the refractive indices. This is also being used in the emerging field of waveguide quantum electrodynamics (see for example [50]). Several authors have derived algorithms to compute the Green’s function for arbitrary complex index profiles [49,51,52].

Notably, the importance of the characteristic Green’s function of waveguides has already been discussed in the electrical engineering literature (see for example Vassallo [26], Chew [28], van Stralen [29], and the references therein). In fact, the mode structure of the fields in a waveguide naturally emerges from poles and branch cuts of the in-plane spatial Fourier transform of its time-harmonic Green’s functions. The Green’s function is thus closely connected to the source-free solutions of the electromagnetic field equations.

While waveguide theory has been extensively studied in the past, the focus has been on longer wavelengths and the results applicable to – or relevant for – x-ray waveguides are sparsely scattered throughout the literature. On the other hand, the unavoidable material absorption prevalent in x-ray waveguides means that the usual Hermitian mode expansions and orthogonality that are commonly used in the non-absorbing case do not apply, and the non-Hermitian generalisations of these mode expansions have received little attention so far.

In this paper, we give a revised self-contained presentation of the nano-optical theory of planar x-ray waveguides. To that end, we present the theory from the ground up, giving derivations where appropriate and also translate some key results from the electrical engineering literature with application to the hard x-ray spectrum. The basic theory is universally applicable to short and longer wavelengths with the most notable features that it consistently assumes complex permittivities and arbitrary layer structures. We derive asymptotic expansions of the electromagnetic Green’s functions in terms of resonant modes using a novel generalised bi-orthogonality relation, which is particularly useful in the presence of strong material absorption. We then apply this general theory and discuss three common experimental geometries for x-ray waveguides (see Fig. 1), namely front-coupling in forward-incidence, evanescent coupling in grazing indicence, and emission from buried dipoles. Using the mode expansions yields compact expressions for some key quantities, such as absorption length, wavelength, and coupling coefficients. In addition to the analytical discussion, we also present numerical code to find resonant waveguide modes and the corresponding coefficients as well as to compute the full Green’s functions in stratified media. The code is made publicly available under the GPLv3 license [53].

Fig. 1. (a) Evanescent coupling in grazing incidence. (b) Front-coupling in forward-incidence. (c) Emission from a dipole within the waveguide.

Download Full Size | PDF

After this introduction, the paper is structured as follows. Section 2 introduces the vector wave equations and Green’s functions in inhomogeneous media. In section 3, we derive solutions to the vector wave equations in the absence of sources and discuss their classification into resonant and non-resonant modes. In section 4 we derive and discuss the fields caused by line sources. Point sources are considered in section 5. In section 6, we discuss a matrix formalism and our numerical approach that can be used to compute the fields and Green’s functions in layered media. In section 7, these general results are then applied. For three important geometries corresponding to prototypical experiments with x-ray waveguides, we derive analytical expressions and present numerical calculations. We conclude the paper in section 8.

2. Equations for the electromagnetic fields

Starting from Maxwell’s equations, we remind ourselves of the derivation of the time harmonic vector wave equations and their solutions in terms of electromagnetic Green’s functions within macroscopic electrodynamics. We then discuss transverse solutions in translationally invariant stratified media.

2.1 Vector wave equations

Consider linear, isotropic, and possibly lossy materials characterized by their position-dependent complex relative permittivity $\epsilon (\vec r)$ and permeability $\mu (\vec r)$. For electromagnetic fields with harmonic time dependence $e^{-i \omega t}$, the macroscopic Maxwell equations read

(2)$$\vec \nabla \times \vec{E} ={-}\vec J_\mathrm{m} + i \omega \vec B ,$$

(3)$$\vec \nabla \times \vec{H} = \vec J_\mathrm{e} - i \omega \vec D .$$

Here, $\vec J_\mathrm {e}$ is the electric and $\vec J_\mathrm {m}$ the magnetic current density. The latter is commonly introduced to make the equations symmetric from a mathematical viewpoint (see for example [28] or [54]). Even though there are no microscopic magnetic charges or true magnetic currents, magnetic dipoles due to microscopic electric current loops or atomic or nuclear transitions can be described by this magnetic current density.

In linear and isotropic media, we have the constitutive relations $\vec {D} = \epsilon _0 \epsilon \vec {E}$ and $\vec {B} = \mu _0 \mu \vec {H}$. Substituting the constitutive relations into the Maxwell equations and taking the curl gives inhomogeneous equations for the electric and magnetic fields, the vector wave equations

(4a)$$\left\{ \mu \vec \nabla \times \frac{1}{\mu} \vec \nabla \times{-} k^2 n^2 \right\} \vec{E} =i \omega \mu_0 \mu \vec{J}_\mathrm{e} - \mu \vec \nabla \times \frac{\vec{J}_\mathrm{m}}{\mu} ,$$

(4b)$$\left\{ \epsilon \vec \nabla \times \frac{1}{\epsilon} \vec \nabla \times{-} k^2 n^2 \right\} \vec{H} = i \omega \epsilon_0 \epsilon \vec{J}_\mathrm{m} + \epsilon \vec \nabla \times \frac{\vec{J}_\mathrm{e}}{\epsilon} ,$$

where $k = \omega /c = \omega \sqrt {\mu _0 \epsilon _0}$ and $n^2 = \epsilon \mu$ is the complex refractive index. Noting the symmetry of (4a) and (4b), a manifestation of the duality principle of electric and magnetic fields, we write

(5)$$\left\{ s_a \vec \nabla \times \frac{1}{s_a} \vec \nabla \times{-} k^2 n^2 \right\} \vec{U}_a = s_a \vec{F}_a,$$

with fields $\vec U_\mathrm {e} = \vec E$, $\vec U_\mathrm {m} = H$, source terms

(6a)$$\vec F_\mathrm{e} = i \omega \mu_0 \vec J_\mathrm{e} - \vec \nabla \times \mu^{{-}1} \vec J_\mathrm{m},$$

(6b)$$\vec F_\mathrm{m} = i \omega \epsilon_0 \vec J_\mathrm{m} + \vec \nabla \times \epsilon^{{-}1} \vec J_\mathrm{e},$$

and $s_\mathrm {e} = \mu$ and $s_\mathrm {m} = \epsilon$. In the following we omit the field index $a$ whenever the meaning is clear.

2.2 Green’s functions

The field produced by a source $\vec F$ can be solved in terms of the dyadic Green’s functions,

(7)$$\vec U_a(\vec r) = \int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_a(\vec r, \vec r') \cdot \vec F_a(\vec r')$$

with the Green’s function being a solution of the vector wave equations for a point source,

(8)$$\left\{ s_a(\vec r) \vec \nabla \times \frac{1}{s_a(\vec r)} \vec \nabla \times{-} k^2 n(\vec r)^2 \right\} \overleftrightarrow{\mathbf{G}}_a(\vec{r}, \vec{r}') = s_a(\vec r) \delta(\vec r - \vec r') \overleftrightarrow{\mathbf{1}} ,$$

where $a = \mathrm {e}, \mathrm {m}$. The Green’s functions are unique when boundary conditions are specified. In an unbounded and absorbing medium, the solutions must vanish at infinity. The validity of (7) can be easily verified by inserting (7) into (5) and using (8). Note that (7) is only valid for $\vec r$ not lying in the support of the source term $\vec F$, because the Green’s function has a singularity at $\vec r = \vec r'$ that requires special treatment [28,54]. To simplify the discussion, we assume here and in the following that the observation point lies outside the source region.

The dyadic Green’s functions for the electric and magnetic field are not independent, but are related by

(9a)$$\overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') ={-}\frac{1}{\epsilon(\vec r) k^2} \delta^3(\vec r - \vec r') - \frac{1}{\epsilon(\vec r) k} \vec \nabla \times \overleftrightarrow{\mathbf{G}}_\mathrm{m} \times \overleftarrow{ \nabla}' \frac{1}{\epsilon(\vec r') k}$$

(9b)$$\overleftrightarrow{\mathbf{G}}_\mathrm{m}(\vec r, \vec r') ={-}\frac{1}{\mu(\vec r) k^2} \delta^3(\vec r - \vec r') - \frac{1}{\mu(\vec r) k} \vec \nabla \times \overleftrightarrow{\mathbf{G}}_\mathrm{e} \times \overleftarrow{ \nabla}' \frac{1}{\mu(\vec r') k} $$

(compare [49], eq. 2.204). While many authors focus on the electric Green’s function for that reason, both have their use in the present discussion.

The fields can be expressed in terms of the electric and magnetic currents. Here, we consider the electric field caused by an electric source current. The general case is derived in appendix D. Substituting $\vec F$ in (7), assuming a purely electric source current, yields

(10)$$\vec E(\vec r) = \frac{i \omega \epsilon_0}{c^2}\int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') \cdot \vec J_\mathrm{e}(\vec r').$$

This allows to directly compute the electric and magnetic field caused by oscillating dipoles, such as $\vec J_\mathrm {e}(\vec r) = - i \omega \vec p \delta ^3(\vec r - \vec r_0)$.

2.3 Transverse solutions in stratified media

The equations simplify significantly in stratified media, and even further for fields that are constant along $y$.

Consider a stratified medium in a coordinate system where the relative permeability $\epsilon$ and relative permeability $\mu$ depend on $z$ only. Further, assume that the fields are constant along $y$. In this case, the term in curly braces in (5) and (8) becomes

(11)$$\left\{ \begin{pmatrix} -s \partial_z s^{{-}1} \partial_z & 0 & s \partial_z s^{{-}1} \partial_x \\ 0 & - \partial_x^2 - s \partial_z s^{{-}1} \partial_z & 0 \\ \partial_x \partial_z & 0 & -\partial_x^2 \\ \end{pmatrix} -k^2 n^2 \right\}$$

so that the $y$-component decouples from the $x$ and $y$ components. In source-free regions, since the fields are coupled by (2) and (3), only 2 out of the 6 components of $\vec E$ and $\vec H$ are independent. A solution for $H_y$ and $E_y$ fully specifies the solution for the remaining components. Solutions with $H_y = 0$ and $E_y$ are often called transverse electric (TE) and transverse magnetic (TM), respectively. The derivatives for the two polarizations are summarized in Table 1. Since the two polarizations decouple, we can analyze them independently.

Table 1. Transverse field components and their derivatives.

View Table

2.4 Power flow

Power flow is given by the real value of the Poynting vector, which is defined as

(12)$$\vec{S} = \frac12 \vec E \times \overline{\vec H} = \frac{-i}{2\omega \mu_0 \bar \mu} \vec E \times (\nabla \times \overline{\vec{E}}).$$

For fields with $\partial _y \vec E = 0$, it separates into TE and TM contributions like

(13)$$\vec{S}^\mathrm{TE} = \frac{-i}{2 \omega \mu_0 \bar \mu} \left\{ E_y \partial_x \bar E_y \hat{x} + E_y \partial_z \bar E_y \hat{z} \right\}$$

(14)$$\vec{S}^\mathrm{TM} = \frac{-i}{2 \omega \mu_0 \bar \mu} \left\{ E_z (\partial_x \bar E_z -\partial_z \bar E_x) \hat{x} - E_x (\partial_x \bar E_z - \partial_z \bar E_x) \hat{z} \right\}.$$

They can be expressed in terms of $U_y$ using Tab. 1.

3. Free fields in stratified media

We discuss the structure of solutions for the source-free case of the vector wave Eqs. (5). Suppose that permittivity and permeability independent of $x$ and $y$ and take constant values $\epsilon _-$, $\mu _-$ and $\epsilon _+$, $\mu _+$ when $z < z_-$ and $z > z_+$, respectively, for some positions $z_-$, $z_+$ in bottom and top half-space.

3.1 Separating the dimensions

The TE- and TM-polarized components of (5) decouple as in (11), so that the $y$-component simplifies to

(15)$$\left\{ s \partial_z s^{{-}1}(z) \partial_z + \partial_x^2 + k^2 n^2(z) \right\} U_y = 0.$$

The components $E_y$ and $H_y$ as well as $E_x$ and $H_x$ are continuous across interfaces, which implies that $U_y$ and $s^{-1} \partial _z U_y$ with $U_y = E_y$ (TE) and $U_y = H_y$ (TM) are continuous.

We use separation of variables on (15) and write $U_y(x,z) = u_{x}(x) u_{z}(z)$, obtaining a pair of eigenvalue problems for each polarization,

(16a)$$- \partial_x^2 u_{x}(x) = k^2 \nu^2 u_{x}(x),$$

(16b)$$\left\{ s \partial_z s^{{-}1} \partial_z + k^2 n^2(z) \right\} u_{z}(z) = k^2 \nu^2 u_{z}(z),$$

such that $u_x$, $(u_x)'$, $u_z$, as well as $s^{-1}(u_z)'$ are continuous, complex-valued functions. The unit-less quantity $\nu \in \mathbb {C}$ is known as the effective (refractive) index and the wave-number $q = k \nu$ as the propagation constant in the context of waveguides. The general solution for $u_{x}$ is easily found to be a superposition of forward-propagating $\exp (i k \nu x)$ and backward-propagating $\exp (-i k \nu x)$ components. The transversal equation will occupy us in the following section.

3.2 Transversal eigenvalue problem

While the longitudinal eigenvalue problem (16a) has a continuous spectrum, the transversal eigenvalue problem (16b) is of Sturm-Liouville form and generally has a discrete and a continuous spectrum of eigenvalues [54]. The discrete spectrum of proper solutions correspond to the guided (or surface-wave) modes, whereas the continuous spectrum of improper solutions correspond to the radiative modes in the context of classical waveguide theory.

It is worthwhile to take a closer look at the operator on the left-hand side of (16b),

(17)$$L = s \partial_z s^{{-}1} \partial_z + k^2 n^2.$$

Importantly, $L$ is generally not self-adjoint because the refractive index $n$ is complex-valued. While self-adjoint Sturm-Liouville operators were extensively studied in the literature, among other things due to their manifold applications in quantum mechanics, the non-self-adjoint theory is much less developed (see for example [55]). In particular, analagous to the situation occuring in the eigensystems of non-symmetric matrices, many tools of self-adjoint operator theory such as the $L^2$-orthogonality of eigenfunctions cannot be applied. Nevertheless, we can find similar relations in many cases.

Analagous to non-symmetric matrices having distinct left and right eigenvectors, the eigenfunctions and adjoint eigenfunctions of $L$ are distinct. In the case that the following integrals converge, integration by parts gives

(18)$$\left\langle \overline{s^{{-}1}u},L v\right\rangle = \left\langle \overline{s^{{-}1}v}, L u\right\rangle + \left[s^{{-}1}(uv' - u'v)\right]_{-\infty}^{\infty},$$

where $\left \langle f, g\right \rangle$ denotes the usual $L^2$ inner product with complex conjugation in the first argument. If $v$ and its derivative vanish sufficiently rapidly compared with $u$ and its derivative, and/or vice versa, the boundary term vanishes. Thus, for for a given eigenfunction $u$ with eigenvalue $\lambda$, the corresponding adjoint eigenfunction is $\overline {s^{-1} u}$, in the sense that for any sufficiently rapidly decaying test function $\psi$ we have

(19)$$\langle \overline{s^{{-}1} u}, L\psi\rangle = \langle \overline{s^{{-}1} \psi}, Lu\rangle = \lambda\langle \overline{s^{{-}1} \psi}, u\rangle = \lambda\langle \overline{s^{{-}1} u}, \psi\rangle.$$

In particular, taking $u,v$ in (18) to be two normalizable eigenfunctions $u_i, u_j$ with eigenvalues $\lambda _i, \lambda _j$ we can show the following relation,

(20)$$(\lambda_i - \lambda_j)\langle \overline{s^{{-}1} u_n}, u_m\rangle =0.$$

Thus, (normalizable) eigenfunctions are orthogonal to adjoint eigenfunctions corresponding to distinct eigenvalues, and we will choose to normalize these eigenfunctions such that

(21)$$\langle \overline{s^{{-}1} u_n}, u_m\rangle = \delta_{nm} .$$

This biorthogonality replaces the usual orthogonality of the eigenfunctions. It should be noted, that this fact is well-established, also in optics (see for example [56]), yet not relevant for non-absorbing waveguides. An alternative way to formalize (21) is the introduction of 2-form such that the eigenfunctions are orthogonal with respect to this form [57].

In general, not all (generalized) eigenfunctions are $L^2$-normalizable, and the integrals may fail to converge. Nevertheless, we will demonstrate that a suitable generalization of (21) exists that allows to normalize certain non-normalizable eigenfunctions.

3.3 Resonant modes

Any formal solution to $(L - k^2 \nu ^2) u = 0$ can be written for $z < z_-$ and $z > z_+$, for some points $z_\pm$ in the bottom and top layer respectively, as

(22)$$u(z) = a_\pm e^{i k p_\pm (z-z_\pm)} + b_\pm e^{{-}i k p_\pm (z-z_\pm)},$$

where $p_\pm ^2 = n_\pm ^2 - \nu ^2$. The coefficients $a_\pm, b_\pm \in \mathbb {C}$ depend on the boundary conditions. An important class of solutions are those that can be written as a single exponential function in both of the regions, so that $a_- = 0$ and $b_+ = 0$. Equivalently, for some points $z_-$ and $z_+$ in the bottom and top layer, respectively,

(23a)$$u'(z) ={-}i k p_- u(z), \text{ for all } z< z_-,$$

(23b)$$u'(z) = i k p_+ u(z), \text{ for all } z \ge z_+.$$

We call the solutions satisfying both (23a) and (23b) resonant modes. They are equivalent to the "resonances" that have been studied by Schenk [39].

Note that, since $p_\pm$ are complex numbers with arbitrary sign, the resonant modes are not necessarily bounded at infinity, and not all resonant modes are normalizable.

3.4 Modified inner product

We can modify the inner product to be applicable to functions belonging to the solution class of resonant modes and to avoid the divergences at infinity. Inspired by a similar definition for quasi-normal modes [57], we define

(24)$$\left\langle u,v\right\rangle_\partial := \int_{z_-}^{z_+} \overline{u(z)} v(z) \mathrm{d} z - \left[ \frac{\bar u v}{ \bar u' / \bar u + v' / v } \right]_{z_-}^{z_+}.$$

More generally, (24) merely requires the functions $u$ and $v$ to satisfy (23a) and (23b) and not necessarily requires them to be solutions of (16b). For convenience, we will refer to (24) as the modified inner product, although we emphasize that it is not, strictly, an inner product. Equation (24) naturally emerges in the semi-spectral representation of the Green’s function in section 4.5 and appendix C. Given (23a) and (23b), this definition does not depend on the choice of $z_\pm$ as long as they are points in the bottom and top layer. Moreover, if both $u$ and $v$ decay for $\left | z \right | \to \infty$, definition (24) reduces to the standard $L^2$ inner product, as the boundary terms correspond to the integrals from $z_\pm$ to $\pm \infty$.

In appendix A we show for any functions $u,v$ which can be written as an exponential function in the top and bottom layers, the following integration by parts identity holds,

(25)$$\left\langle u',v\right\rangle_\partial + \left\langle u,v'\right\rangle_\partial =0.$$

Thus, the adjoint of $L$ is well defined with respect to the modified inner product, and in particular every resonant mode has a corresponding adjoint eigenfunction with respect to the modified inner product,

(26)$$\left\langle \overline{s^{{-}1} u}_i, L v\right\rangle_\partial = \lambda_i \left\langle \overline{s^{{-}1} u}_i, v\right\rangle_\partial.$$

This in turn shows that even for divergent resonant modes, we have the following bi-symmetry relation

(27)$$\left\langle \overline{s^{{-}1} u}_m,L u_n\right\rangle_\partial = \left\langle \overline{s^{{-}1} u}_n, L u_m\right\rangle_\partial.$$

The bi-symmetry relation guarantees that resonant modes, even diverging ones, are bi-orthogonal with respect to (24), and thus we can normalize all resonant modes to obey

(28)$$\left\langle \overline{ s^{{-}1} u_n}, {u_m}\right\rangle_\partial = \delta_{m,n}.$$

3.5 Wronskian resonance condition

Next, we present a practical method to classify whether a given complex number $\nu \in \mathbb {C}$ admits a resonant-mode solution with that eigenvalue $\nu$. To that end, we use a relatively general method (see [58], Ch IV.6) to relate the resonance condition (23a) and (23b) to a condition about the solutions of a regular Sturm-Liouville initial value problem on the compact interval $(z_-, z_+)$. It leads to a one-to-one correspondence of the complex roots of a holormorphic function, and the set of resonant modes.

There are at most a discrete set of resonant modes, and possibly none at all [39, Thm. 2.12]. Certainly, it is clear that (23a) and (23b) can only simultaneously be fulfilled for special values of $\nu$. However, since the initial value problem on the interval $(z_-, z_+)$ has a solution for all possible initial values and this solution is unique [55], we can always find unique solutions up to a global scale that satisfy either (23a) or (23b). We construct solutions $u_+(z; \nu )$ and $u_-(z; \nu )$ by setting

(29)$$\begin{cases} u_\pm(z_\pm; \nu) = 1 \\ \partial_z u_\pm(z_\pm; \nu) ={\pm} ik p_\pm(\nu). \end{cases}$$

These solutions are resonant modes if and only if they are linearly dependent, that is, if their Wronskian vanishes,

(30)$$\operatorname{Wron}\left\{u_-(z), u_+(z)\right\} = u_-(z) \frac{\partial u_+(z)}{\partial_z} - \frac{\partial u_-(z)}{\partial_z} u_+(z) = 0.$$

Here, the identically vanishing Wronskian implies linear dependence because the solutions are different from zero throughout the interval $(z_-, z_+)$ (see [59]). Instead of the bare Wronskian, it is useful to work with the weighted Wronskian defined as

(31)$$W(\nu) = \frac{1}{s} \left( u_-(\nu) \frac{\partial u_+(\nu)}{\partial_z} - \frac{\partial u_-(\nu)}{\partial_z} u_+(\nu) \right),$$

which is independent of $z$, as can be easily shown by inserting (16b).

We thus have a one-to-one correspondence between the set of resonant modes and the complex roots of the weighted Wronskian $W(\nu )$.

The weighted Wronskian takes particularly simple forms at $z_-$ or $z_+$,

(32)$$\begin{aligned} W &= s^{{-}1}(z_+) i k p_+ u_-(z_+) - s^{{-}1}(z_+) \partial_z u_-(z_+)\\ &= s^{{-}1}(z_-) i k p_- u_+(z_-) + s^{{-}1}(z_-) \partial_z u_+(z_-). \end{aligned}$$

We will see in section 6 how it can be readily computed in the transfer matrix formalism.

A resonant mode $u_\pm$ for fixed $\nu _m$, with $W(\nu _m) = 0$, is unique up to a global prefactor. We normalize the solutions by setting

(33)$$u_m(z) := \frac{u_\pm(z)}{\left\langle \overline{ s^{{-}1} u_\pm}, u_\pm\right\rangle_\partial}.$$

Note that (33) may give a different overall sign for $u_+$ and $u_-$, so a choice has to be made for consistency.

3.6 Riemann sheets and branch cuts

The choice of sign in the mappings $\nu \mapsto p_\pm (\nu )$ that solve $n_\pm ^2 - \nu ^2 = p_\pm ^2$ for the resonance condition (23a) and (23b) dictates the asymptotic behavior of the corresponding modes and thereby what kind of modes are found when solving $W(\nu ) = 0$. Each of the two quadratic expressions has a 2-valued solution set for every complex number $\nu$, the Riemann surface of the complex square root. In the same way, $W(\nu )$ can be seen as a 4-valued Riemann surface, discriminated by the two signs for $p_+(\nu )$ and $p_-(\nu )$. As a consequence, the complex function $W(\nu )$ has 4 algebraic branch points where the mode index coincides with the material index of the bottom or top layer, $\nu \in \{ \pm n(z_-), \pm n(z_+)\}$. In general, resonant modes can be located on all 4 sheets of the Riemann surface. A single-valued function $W$ is obtained by choosing a single-valued complex square root, which introduces 4 branch cuts, as shown in Fig. 2.

Fig. 2. Sketch of branch points (stars) and branch cuts (lines) for non-degenerate outer layer indices. (a) Natural branch cuts for bounded solutions. (b) Vertical branch cuts. The branch points, which are the refractive indices of the outer layers, would be indistinguishable from $1$ at this scale for hard x-rays.

Download Full Size | PDF

Clearly, a resonant mode is normalizable if and only if $\textrm{Im}\{ p_- \} > 0$ and $\textrm{Im}\{ p_+\} > 0$. Setting $p_\pm = p_\mathrm {dec}(\nu ; n_\pm )$, where

(34)$$p_\mathrm{dec}(\nu; n) = i \sqrt{\nu^2 - n^2}$$

with the principal complex square root, ensures that $\textrm{Im}\{p_\pm \} \ge 0$.

Each single-valued mapping $\nu \mapsto W(\nu )$ is a holomorphic function in the entire complex plane $\mathbb {C}$ except for the two branch cuts. This is guaranteed because solutions of the regular Sturm-Liouville initial value problem are complex differentiable everywhere as a function of their data and initial values [55]. The branch cuts enter through the initial values in (29).

In addition to the decaying choice of branch cut made in (34), there are other useful choices.

(35)$$p_\mathrm{out}(\nu) := \sqrt{n^2 - \nu^2}$$

with the standard choice of complex square root guarantees $\textrm{Re}\{p\} \ge 0$, so that the corresponding modes are outwards-propagating. It is easy to see that for absorbing materials, $ \operatorname{Im}\{v\}>0 $, both definitions coincide on the real axis so that $p_\mathrm {out}(\nu ) = p_\mathrm {dec}(\nu )$ for $ \operatorname{Im}\{v\}=0 $.

Neither definitions, (34) and (35), are ideal for numerical purposes, because they define branch cuts which are curved in the complex $\nu$ plane. The inverse Fourier integral along the branch cut involves a plane-wave factor of the form $e^{i k \nu x}$, and a curved branch cut will result in a rapdily oscillating integral, which is challenging to numerically evaluate.

Setting

(36)$$p_\mathrm{vert}(\nu) := i \sqrt{-i (n-\nu)} \sqrt{-i(n+\nu)}$$

results in branch cuts that are parallel to the imaginary axis. Along such a branch cut, the real part of $\nu$ is constant, and thus the plane wave factor $e^{i k \nu x}$ in the inverse Fourier integral will be purely exponentially decaying, and the integral will rapidly converge.

Finally, the Riemann surface can be unfolded to eliminate the branch points using conformal mappings [60,61]. Setting

(37)$$p_\pm(\xi) = \xi \pm \frac{n_+^2 - n_-^2}{4 \xi},$$

eliminates the branch cuts [60] for $W(\xi ) = W(\nu (\xi ))$. The other quantities then become

(38)$$p(\xi) = i \sqrt{- p_+^2(\xi) -(n^2 - n_+^2)}, \quad \nu(\xi) = \sqrt{n_+^2 - p_+^2(\xi)}.$$

Using the mapped coordinates, $W(\xi ) = 0$ will be solved by all resonant modes.

3.7 Guided and leaky modes

In the literature, the resonant modes are classified into guided and leaky modes (see for example [62] for a recent discussion). Guided modes are the normalizable resonant modes. In contrast to guided modes, leaky modes are diverging for either $z \to -\infty$ (bottom-leaky), $z \to \infty$ (top-leaky), or both (double-leaky).

Depending on the choice of branch cuts, the complex roots of $W(\nu )$ correspond to different types of resonant modes. For the "natural" choice defined in (34) (see Fig. 2(a)), the transversal solutions are bounded and the roots hence only correspond to guided modes. For other branch cuts, such as (36) (see Fig. 2(b)), on the other hand, the solution set in general consists of both guided and leaky modes.

3.8 Radiative modes

The non-nomalizable but bounded improper solutions of $(L - k^2 \nu ^2) u = 0$ form the set of radiative (or radiation) modes. As discussed by Vassallo [26], they can be divided into two groups that are either radiative in the top or bottom layer. The groups become degenerate when top and bottom materials are identical, which can be avoided by changing one of the indices by an infinitesimal imaginary amount.

A solution is radiative in the top or bottom layer, if $\textrm{Im}\{p_+\} = 0$ or $\textrm{Im}\{p_-\} = 0$, respectively. Setting $p_+ = p_\mathrm {out}(\nu )$ ensures that $\textrm{Im}\{ p_+\} = 0$ when $\nu \in \left \{ \nu : \textrm{Im}\{ \nu ^2 \} = \textrm{Im} \{n_+^2\}, \left | \textrm{Re} \{ \nu \} \right | \le \textrm{Re} \{ n_+ \} \right \}$, and $p_-$ accordingly. The mode indices $\nu$ of radiative modes have $ \operatorname{Im}\{v\}=\operatorname{Re}\{v\} \operatorname{Im}\left\{n_{ \pm}\right\} / \operatorname{Re}\left\{n_{ \pm}\right\} $. For any such $\nu$, the radiative modes are simply given by $u_\pm$ as defined in (29), and using (34) for the non-radiative layer. If $n_+ \neq n_-$, no solution can be radiative in the top and bottom layer simultaneously.

3.9 Mode dispersion

So far we have looked at the strictly monochromatic case only. X-ray light sources such as synchrotrons and x-ray free electron lasers, however, produce pulsed light, which has a finite bandwidth, typically in the order of 1 eV. When discussing pulses, it is important to know how the mode indices depend on the angular frequency $\omega$. The group velocity is approximated to first order by

(39)$$v_\mathrm{g} = \frac{c}{\operatorname{Re} \left \{\nu_m + \omega \frac{\partial \nu_m}{\partial \omega} \right\}} \approx \frac{c}{\operatorname{Re} \nu_m}.$$

In appendix B, we show that

(40)$$\nu_m \frac{\partial \nu_m}{\partial \omega} = \frac{1}{2}\left\langle \frac{\partial \ln(s)}{\partial \omega} (\nu_m^2 - n^2) \right\rangle_{u_m} + \frac{c^2}{2\omega^2} \left\langle \frac{\partial \ln(s)}{\partial \omega}\right\rangle_{\partial_z u_m} + \left\langle \frac{n^2 - \nu_m^2}{\omega}\right\rangle_{u_m} + \left\langle n \frac{\partial n}{\partial \omega}\right\rangle_{u_m},$$

for any resonant mode with index $\nu _m(\omega )$. Here, the ‘expectation value’ is defined as $\left \langle A \right \rangle _u = \langle \overline {s^{-1} u}, A u \rangle _\partial$. For TE-modes, we have $s \equiv \mu \equiv 1$ so that the first 2 terms vanish. Using $\nu _m \approx n \approx 1$, this can be well approximated by

(41)$$\frac{\partial \nu_m(\omega)}{\partial \omega} \approx \int \left\{ \frac{2 (n - \nu)}{\omega} + \frac{\partial n(z,\omega)}{\partial\omega} \right\} \frac{u_m(z,\omega)^2}{s(z, \omega)} \mathrm{d} z.$$

As a consequence, the mode dispersion is on the scale of the material dispersion and hence small in the absence of absorption edges or other resonances. Corrections to the group velocity can be computed using (41).

It should be noted that the decomposition into guided and leaky modes is only valid locally in the frequency spectrum. For broad-band pulses, or in the vicinity of strong material resonances, the mode classification is no longer unique [29].

4. Line sources

The fields in the presence of sources can be conveniently solved in terms of the dyadic Green’s functions, defined in (8). An important special case are line sources, which extend homogeneously along one dimension, here $y$. This case was previously discussed in [26].

4.1 Two-dimensional Green’s function

When the source terms are constant over $y$, the integral over $y$ in (7) can be performed and the Green’s function for line sources can be defined as

(42)$$\overleftrightarrow{\mathbf{G}}^{\mathrm{ls}}(x-x', z, z') = \int \mathrm{d} y' \overleftrightarrow{\mathbf{G}}(x-x', y-y', z,z').$$

While the 3-dimensional Green’s functions have the dimension of inverse length, these 2-dimensional Green’s functions are dimensionless. By integrating both sides of (8) over $y$ and using the fact that the functions vanish at infinity, it follows that

(43)$$\left\{ \mathrm{D}(z) - k^2 n(\vec r)^2 \right\} \overleftrightarrow{\mathbf{G}}^\mathrm{ls}_a(\vec{r}, \vec{r}') = s_a(\vec r) \delta(x-x') \delta(z-z') \overleftrightarrow{\mathbf{1}} ,$$

with

(44)$$\mathrm{D}(z) = \begin{pmatrix} -s \partial_z s^{{-}1} \partial_z & 0 & s \partial_z s^{{-}1} \partial_x \\ 0 & -\partial_x^2 - s \partial_z s^{{-}1} \partial_z & 0 \\ \partial_x \partial_z & 0 & -\partial_x^2 \\ \end{pmatrix}.$$

Consequently, the $yy$-components decouple as in (11) and we can analyze the two scalar $yy$-components independently. Defining $g_a = \overleftrightarrow {\mathbf {G}}^\mathrm {ls}_{a,yy}$, we obtain

(45)$$\left\{ \partial_x^2 + s \partial_z s_a^{{-}1} \partial_z + k^2 n^2(z) \right\} g_a(x - x', z, z') ={-} s_a(z) \delta(z-z') \delta(x-x') .$$

This shows that the $yy$-components of the Green’s function for line sources are simply the two-dimensional Green’s functions of the medium. In the following, we will discuss solutions of (45). We will drop the index $a$ for convenience and keep in mind that the expressions are valid for both electric and magnetic Green’s functions.

4.2 Homogeneous Green’s function

We start by deriving the two-dimensional Green’s function for homogeneous media $g^{(0)}(x-x', z-z')$. Writing the Green’s function as a Fourier transform and exploiting its rotational symmetry, we obtain

(46)$$\begin{aligned} g^{(0)}(\vec \rho) &= \int \frac{\mathrm{d}^2 q}{(2 \pi)^2} e^{i \vec q \cdot \vec \rho} g^{(0)}(q)\\ &=\int_0^\infty \frac{\mathrm{d} q}{2 \pi} J_0(q \rho) q g^{(0)}(q), \end{aligned}$$

where $J_0$ is the Bessel function of first kind. Inserting (46) into (45) and using that $s$ and $n$ are constants, we immediately obtain the spectral representation

(47)$$g^{(0)}(q) = \frac{s}{q^2 - k^2 n^2}.$$

Inserting (47) into (46) yields

(48)$$g^{(0)}(\vec \rho) = \frac{i s}{4} H_0^{(1)}(n k \rho),$$

where $H_0^{(1)}$ is the Hankel function of first kind.

4.3 Semi-spectral representation

The general stratified medium can be solved similarly using Fourier transforms. While the solution is well known, we briefly remind ourselves of its construction. Taking the one-dimensional Fourier transform of (45) with respect to $x-x'$ yields

(49)$$\left\{ s \partial_z s^{{-}1} \partial_z + k^2 n^2(z) - q^2 \right\} g(q, z, z') ={-} s(z) \delta(z-z'),$$

where $g(q,z, z')$ is a continuous function of $z$. Identifying $q^2$ with $k^2 \nu ^2$, Eq. (49) becomes identical to the source-free eigenvalue problem (16b), except for the $\delta$-contribution at $z=z'$. As a consequence, we can construct a bounded solution piece-wise for $z < z'$ and $z > z'$,

(50)$$g(q, z, z') = a_- u_-(z; q) \Theta(z' - z) + a_+ u_+(z; q) \Theta(z - z'),$$

where we have inserted the solutions defined in (29). The symbol $\Theta (z)$ denotes the unit step function. It remains to find the two coefficients $a_\pm$. In the following the $q$-dependence is not displayed. Integrating (49) from $z'-0$ to $z'+0$ and using that $g(q,z)$ is continuous yields

(51)$$\partial_z g(q, z'+0, z') -\partial_z g(q, z'-0, z') ={-} s(z') ,$$

which, using (50), becomes

(52)$$a_+ \partial_z u_+(z') - a_- \partial_z u_-(z') ={-} s(z').$$

Moreover, we have from continuity that $a_-u_-(z') = a_+u_+(z')$, so that

(53)$$\frac{- a_- a_+ }{s(z')} \left(u_-(z') \partial_z u_+(z') - u_+(z') \partial_z u_-(z')\right) = a_+ u_+(z') = a_- u_-(z').$$

We obtain $a_\pm = -u_\mp (z')/W$, where we have inserted $W$ as defined in (31), which is independent of $z$. Inserting the coefficients into (50) finally gives

(54)$$g(q, z, z') ={-}\frac{1}{W} \left\{u_+(z') u_-(z) \Theta(z' - z) + u_-(z') u_+(z) \Theta(z - z') \right\},$$

which we call the semi-spectral representation.

As an example, consider a homogeneous environment with refractive index $n$. Equation (54) then yields

(55)$$g^{(0)}(q,z,z') = \frac{is}{2 k p} e^{ikp \left| z -z' \right|},$$

the Fourier-transform of (48).

4.4 Derivation of real-space solutions

The real-space form of the Green’s function can be written as

(56)$$g(x-x', z, z') = \int_{-\infty}^\infty \frac{\mathrm{d} q}{2\pi} e^{i q (x-x')} g(q, z, z').$$

Numerical integration is difficult due to the rapidly oscillating integrand. We can however extend the integral over a closed contour and exploit the structure of $g(q, z, z')$. The Fourier-space Green’s function is given as a rational function (54), where the numerator $u_+(q) u_-(q)$ is holomorphic everywhere and the denominator $W(q)$ is holomorphic everywhere except at its 4 algebraic branch points. Since $W(q)$ is the denominator, its complex roots produce poles in $g(q)$.

For $x > x'$ ($x < x'$), we close the contour in the upper half (lower half) of the complex $q$ plane (see Fig. 2) such that that arc portion of the integral vanishes due to the exponential decay of $\exp (i q (x-x'))$. We must, however, be careful to choose the contours in a way to avoid the poles and branch cuts of the integrand.

Note, however, that the branch cuts can be chosen (compare Section 3.6) and its choice influences which poles are exposed. The branch cut corresponding to decaying solutions is curved and only exposes the poles corresponding to guided modes, a finite number $N_\mathrm {guided}$. We deform the contour so that it runs next to the branch cut from infinity to the branch point and back to infinity on the other side (see Fig. 2). We obtain

(57)$$\begin{aligned} g(x-x', z, z') &= \sum_{m=1}^{N_\mathrm{guided}} i \operatorname{Res} \{ g; q = k \nu_m \} e^{i k \nu_m \left| x- x' \right|}\\ &+ r_\mathrm{top}(x-x', z, z')+ r_\mathrm{bot}(x-x', z, z'). \end{aligned}$$

The non-resonant terms $r_\mathrm {top}$ and $r_\mathrm {bot}$ correspond to the parts of the contour wrapping around the branch cuts caused by the top and bottom layer, respectively. Since the branch cuts coincide with the radiation modes, these terms can be understood as contribution by radiation modes. Performing the integrals for these non-resonant terms is non-trivial, because the branch cuts run almost parallel to the real axis, causing rapid oscillations. However, in many cases the guided modes alone provides good asymptotic approximations of the Green’s function.

But not all systems even support guided modes. To obtain asymptotic expressions for those systems, we deform the branch cuts to be parallel to the imaginary axis. This does not change the value of the integral, because it leaves the real axis unaffected. It changes, however, the non-resonant contributions to the integral and exposes the poles corresponding to leaky modes, if they exist. We obtain

(58)$$\begin{aligned} g(x-x', z, z') &= \sum_m i \operatorname{Res} \{ g; q = k \nu_m \} e^{i k \nu_m \left| x- x' \right|}\\ &+ \tilde r_\mathrm{top}(x-x', z, z')+ \tilde r_\mathrm{bot}(x-x', z, z'), \end{aligned}$$

where now $m$ runs over the (possibly infinite) number of exposed guided and leaky modes. The non-resonant terms $\tilde r_\mathrm {top}$ and $\tilde r_\mathrm {bot}$ now correspond to the contour integrals following the deformed vertical branch cuts, which makes them easier to compute numerically. Due to the change of branch cut, these non-resonant modes are diverging for $\left | z \right | \to \infty$. We call them diverging radiation modes. Given that the total Green’s function decays at infinity and the deformation of the branch cuts does not change the value of the integral (56), the leaky modes and diverging radiation modes must interfere destructively outside the central layers in such a way that the total Green’s function does not diverge [62].

4.5 Computation of the residues

The residues in (58) can be expressed in a more useful way. From (54) and using that numerator and denominator are holomorphic functions of $q$, we can read off

(59)$${\textrm{Res}} \{ g; q = k \nu_m \} ={-} \frac{u_+(z') u_-(z)}{\partial_q W(q=k \nu_m)},$$

where we used that $u_+ \propto u_-$ for zeros of $W$. As derived in appendix C, we obtain

(60)$$\partial_q W(\nu_m) ={-}2 \nu_m k \left\langle \overline{s^{{-}1} u_-}, u_+\right\rangle_\partial,$$

where $\left \langle \cdot, \cdot \right \rangle _\partial$ is the modified inner product as defined in (24). It then follows that

(61)$${\textrm{Res}} \{ g; q = k \nu_m \} = \frac{u_+(z') u_-(z)}{ 2 \nu_m k \left\langle \overline{s^{{-}1} u_-}, u_+\right\rangle_\partial}.$$

Using $u_+ \propto u_-$, we can insert the normalized modes to obtain

(62)$${\textrm{Res}} \{ g; q = k \nu_m \} = \frac{u_m(z') u_m(z)}{ 2 \nu_m k}.$$

Inserting (62) into (58) or (57), we have

(63)$$g(x-x', z, z') = i \sum_m \frac{u_m(z') u_m(z)}{ 2 \nu_m k} e^{i k \nu_m \left| x- x' \right|} + \text{non-resonant},$$

where the sum goes either over the guided modes or all the resonant modes and the non-resonant part changes accordingly. Knowing the resonant terms, the semi-spectral representation of (63) is easily seen to be

(64)$$g(q, z, z') = \sum_m \frac{u_m(z') u_m(z)}{ q^2 - (\nu_m k)^2} + \text{non-resonant},$$

where the non-resonant terms are the Fourier-transforms of the non-resonant terms in (63). Note that due to its branch points, $g(q, z, z')$ is not meromorphic as a function of $q$ so that Mittag-Leffler’s theorem is not applicable and the non-resonant terms are in general not holomorphic.

4.6 Dyadic Green’s functions for line sources

From (43) it follows that the dyadic Green’s function for line sources has 5 non-vanishing components,

(65)$$\overleftrightarrow{\mathbf{G}}^\mathrm{ls}_a = \begin{pmatrix} G^\mathrm{ls}_{a,xx} & 0 & G^\mathrm{ls}_{a,xz} \\ 0 & G^\mathrm{ls}_{a,yy} & 0 \\ G^\mathrm{ls}_{a,zx} & 0 & G^\mathrm{ls}_{a,zz} \\ \end{pmatrix}.$$

The index $a$ indicates the electric ($a = \mathrm {e}$) and magnetic ($a = \mathrm {m}$) Green’s functions.

While we have derived the $yy$-component in the previous section, $G^\mathrm {ls}_{a,yy} \equiv g_a$, we can express the remaining components of (65) in terms of $g_a$. Integrating both sides of (9a) over $y$ and inserting $\overleftrightarrow {\mathbf {G}}^\mathrm {ls}_\mathrm {m} = g_\mathrm {m} \hat {y} \hat {y}$, we obtain, neglecting the deltas,

(66a)$$G^\mathrm{ls}_{\mathrm{e},xx}(x-x', z, z') = \frac{\partial_z \partial_{z'} g_{\mathrm{m}}(x-x', z, z')}{k^2 \epsilon(z) \epsilon(z')}$$

(66b)$$G^\mathrm{ls}_{\mathrm{e},xz}(x-x', z, z') ={-} \frac{\partial_z \partial_{x'} g_{\mathrm{m}}(x-x', z, z')}{k^2 \epsilon(z) \epsilon(z')}$$

(66c)$$G^\mathrm{ls}_{\mathrm{e},zx}(x-x', z, z') ={-} \frac{\partial_x \partial_{z'} g_{\mathrm{m}}(x-x', z, z')}{k^2 \epsilon(z) \epsilon(z')}$$

(66d)$$G^\mathrm{ls}_{\mathrm{e},zz}(x-x', z, z') = \frac{\partial_x \partial_{x'} g_{\mathrm{m}}(x-x', z, z')}{k^2 \epsilon(z) \epsilon(z')},$$

and accordingly for $\overleftrightarrow {\mathbf {G}}^\mathrm {ls}_\mathrm {m}$.

Taking the Fourier-transform with respect to $x-x'$ yields, separating the TE and TM components,

(67)$$\overleftrightarrow{\mathbf{G}}_\mathrm{e}^\mathrm{TE}(q) = g_{\mathrm{e}}(q) \hat{y} \hat{y}$$

(68)$$\overleftrightarrow{\mathbf{G}}_\mathrm{e}^\mathrm{TM}(q) = \frac{1}{k^2 \epsilon(z) \epsilon(z')} \left \{ \partial_z \partial_{z'} \hat{x} \hat{x} + iq \partial_{z} \hat{x} \hat{z} - iq \partial_{z'} \hat{z} \hat{x} + q^2 \hat{z} \hat{z} \right\} g_{\mathrm{m}}(q).$$

We can read off the residues of the remaining components, complementing (62), at the resonances $q = k \nu _m^\mu$,

(69a)$${\textrm{Res}}\{G^\mathrm{ls}_{\mathrm{e},xx}; q = k \nu_{\mathrm{m},m}\} = \frac{1}{k^2 \epsilon(z) \epsilon(z')} \frac{\partial_{z'} u_{\mathrm{m},m}(z') \partial_z u_{\mathrm{m},m}(z)}{ 2 \nu_{\mathrm{m},m} k}$$

(69b)$${\textrm{Res}}\{G^\mathrm{ls}_{\mathrm{e},xz}; q = k \nu_{\mathrm{m},m}\} = \frac{i}{k^2 \epsilon(z) \epsilon(z')} \frac{ u_{\mathrm{m},m}(z') \partial_z u_{\mathrm{m},m}(z)}{ 2 }$$

(69c)$${\textrm{Res}}\{G^\mathrm{ls}_{\mathrm{e},zx}; q = k \nu_{\mathrm{m},m}\} ={-} \frac{i}{k^2 \epsilon(z) \epsilon(z')} \frac{\partial_{z'} u_{\mathrm{m},m}(z') u_{\mathrm{m},m}(z)}{ 2 }$$

(69d)$${\textrm{Res}}\{G^\mathrm{ls}_{\mathrm{e},zz}; q = k \nu_{\mathrm{m},m}\} = \frac{k \nu^H_m}{k^2 \epsilon(z) \epsilon(z')} \frac{ u_{\mathrm{m},m}(z') u_{\mathrm{m},m}(z)}{ 2},$$

and accordingly for $\overleftrightarrow {\mathbf {G}}^\mathrm {ls}_\mathrm {m}$.

To compute the semi-spectral representation of the full dyadic Green’s function, one can first compute $G_{\mathrm {m}, yy}(q)$ and $G_{\mathrm {e},yy}(q)$ and insert them into (67) and (68). To compute asymptotic expansions for the real-space Green’s functions, the following strategy can be used. First, find the complex roots of $W_a(\nu )$ corresponding to the guided modes and as many of the leaky modes as needed for both $E$ and $H$. Second, compute the corresponding residues for the 5 non-vanishing components. Finally, substitute the residues into (58).

5. Point sources

The Green’s functions for line sources is useful to describe 2-dimensional problems. We now discuss the 3-dimensional Green’s function, the solution of (8), which describes the fields due to a point source. It can be computed from the 2-dimensional one by exploiting the cylindrical symmetry [52]. Similar to the 2-dimensional Green’s function, we start with the semi-spectral representation and proceed with asymptotic expansions of the real-space function.

5.1 Symmetry considerations

The 3-dimensional Green’s function takes the form $\overleftrightarrow {\mathbf {G}}(\vec r, \vec r') = \overleftrightarrow {\mathbf {G}}(\vec r_\parallel - \vec r'_\parallel, z, z')$ because of the in-plane symmetry. Writing the in-plane difference vector $\vec r_\parallel - \vec r'_\parallel$ in polar coordinates $(\rho, \phi )$, we obtain

(70)$$\overleftrightarrow{\mathbf{G}}(\rho, \phi, z, z') = \mathrm{R}_z(\phi) \overleftrightarrow{\mathbf{G}}(\rho, \phi = 0, z, z') \mathrm{R}_z(\phi)^t,$$

with rotation matrix

(71)$$\mathrm{R}_z(\phi) = \begin{pmatrix} \cos \phi & -\sin \phi & 0 \\ \sin \phi & \cos \phi & 0 \\ 0 & 0 & 1 \end{pmatrix}.$$

We write the Green’s function as a Fourier transform with respect to $\vec r_\parallel - \vec r'_\parallel$, assuming it exists,

(72)$$\begin{aligned} \overleftrightarrow{\mathbf{G}}(\rho, 0, z, z') &= \int \frac{\mathrm{d}^2q}{(2\pi)^2} e^{i \vec q \cdot (\vec r_\parallel{-} \vec r'_\parallel)} \overleftrightarrow{\mathbf{G}}(q, \phi_q, z, z')\\ &= \int_0^\infty \frac{\mathrm{d} q}{2\pi} q \int_0^{2\pi} \frac{\mathrm{d} \phi_q}{2\pi} e^{i q \rho \cos \phi_q}\overleftrightarrow{\mathbf{G}}(q, \phi_q, z, z'). \end{aligned}$$

Substituting (72) into (8) yields

(73)$$\left \{ \mathrm{R}_z(\phi_q)\mathrm{D}(q) \mathrm{R}_z(\phi_q)^t - k^2 n^2(z) \right\} \overleftrightarrow{\mathbf{G}}(q, \phi_q, z, z') = s(z) \delta(z-z'),$$

with

(74)$$\mathrm{D}(q) = \begin{pmatrix} -s \partial_z s^{{-}1} \partial_z & 0 & s \partial_z s^{{-}1} iq \\ 0 & q^2 - s \partial_z s^{{-}1} \partial_z & 0 \\ iq \partial_z & 0 & q^2 \\ \end{pmatrix}$$

being the Fourier transform of (11). Applying the rotation matrix to both sides of the equation, we obtain

(75)$$\left \{ \mathrm{D}(q) - k^2 n^2(z) \right\} \overleftrightarrow{\mathbf{G}}(q, \phi_q = 0, z, z') = s(z) \delta(z-z').$$

So the rotated Fourier-space equation is just the equation for the line source (43), and consequently $\overleftrightarrow {\mathbf {G}}(q, \phi _q = 0, z, z') \equiv \overleftrightarrow {\mathbf {G}}^\mathrm {ls}(q, z, z')$.

5.2 Derivation of the 3d real-space Green’s function

We can thus insert the known expressions for the 2-dimensional Green’s function into (72) and perform the integral. The Fourier-space Green’s functions have 5 non-vanishing components (65). Computing the matrix-product and performing the angular integrals in (72) yields [52]

(76)$$\begin{aligned}G_{xx}(\rho, z, z') &= \int_0^\infty \frac{\mathrm{d} q}{2\pi} \bigg\{ \left[q J_0(q\rho) - \frac{J_1(q\rho)}{\rho}\right] G^\mathrm{ls}_{xx}(q, z, z')\\ &+ \frac{J_1(q \rho)}{\rho} G^\mathrm{ls}_{yy}(q, z, z') \bigg\} \end{aligned}$$

(77)$$G_{xz}(\rho, z, z') = i \int_0^\infty \frac{\mathrm{d} q}{2\pi} q J_1(q\rho) G^\mathrm{ls}_{xz}(q, z, z')$$

(78)$$G_{zz}(\rho, z, z') = \int_0^\infty \frac{\mathrm{d} q}{2\pi} q J_0(q\rho) G^\mathrm{ls}_{zz}(q, z, z').$$

The remaining components can be obtained by replacing $x$ with $y$ in the indices.

Computing the integrals involving the Bessel functions is difficult due to their rapid oscillations and the divergence of the Bessel functions for complex $q$ complicates contour integration. This can be avoided by splitting the Bessel function into Hankel functions [52]. Using

(79)$$J_n(z) = \frac{1}{2} \left\{ H_n^{(1)}(z) + H_n^{(2)}(z) \right\},$$

the reflection formula, $H_n^{(2)}(z e^{-i \pi i}) = - e^{i n \pi i} H_n^{(1)}(z)$, and the symmetry and anti-symmetry with respect to $q$ of the diagonal and non-diagonal elements of the Green’s tensor, respectively, we obtain

(80)$$\begin{aligned}G_{xx}(\rho, z, z') &= \frac12 \int_{-\infty}^\infty \frac{\mathrm{d} q}{2\pi} \bigg\{ \left[q H^{(1)}_0(q\rho) - \frac{H^{(1)}_1(q\rho)}{\rho}\right] G^\mathrm{ls}_{xx}(q, z, z')\\ &+ \frac{H^{(1)}_1(q \rho)}{\rho} G^\mathrm{ls}_{yy}(q, z, z') \bigg\} \end{aligned}$$

(81)$$G_{xz}(\rho, z, z') = \frac{i}{2} \int_{-\infty}^\infty \frac{\mathrm{d} q}{2\pi} q H^{(1)}_1(q\rho) G^\mathrm{ls}_{xz}(q, z, z')$$

(82)$$G_{zz}(\rho, z, z') = \frac12 \int_{-\infty}^\infty \frac{\mathrm{d} q}{2\pi} q H^{(1)}_0(q\rho)G^\mathrm{ls}_{zz}(q, z, z').$$

The path of integration adopted by Sommerfeld runs slightly above the negative, and slightly below the positive real axis.

Even though the expressions look relatively complicated, the trace takes a much simpler form,

(83)$$\operatorname{Tr}\overleftrightarrow{\mathbf{G}} = \sum_{j = x,y,z} \int_{-\infty}^\infty \frac{\mathrm{d} q}{2\pi} q H^{(1)}_0(q\rho)G^\mathrm{ls}_{jj}(q, z, z').$$

Moreover, two of the three terms of the $xx$ and $yy$-components decay with $1/\rho$ and can be neglected in many cases.

Similarly to the 2-dimensional case, we use contour integration to expand the integrals in terms of their residues. Here we give the derivation for the $zz$-component. Closing the contour in (82), the integral has two contributions: due to the poles and due to the branch cuts. We obtain

(84)$$\begin{aligned}G_{zz}(\rho, z, z') &= \frac{i}{2} \sum_m k \nu_m H_0^{(1)}(k \nu_m \rho) {\textrm{Res}}\{G^\mathrm{ls}_{zz}; q = k \nu_m\}\\ &+ \frac12 \tilde r_\mathrm{bot}(\rho) + \frac12 \tilde r_\mathrm{top}(\rho) \end{aligned}$$

where $\tilde r_\mathrm {bot}$ are the contributions due to the branch cuts. Inserting the residues for the electric Green’s function, we are left with

(85)$$\begin{aligned} G_{\mathrm{e},zz}(\rho, z, z') &= \frac{i}{4\epsilon(z) \epsilon(z')} \sum_m \nu_{\mathrm{m},m}^2 H_0^{(1)}(k \nu_{\mathrm{m},m} \rho) u_{\mathrm{m},m}(z') u_{\mathrm{m},m}(z)\\ &+ \frac12 \tilde r_\mathrm{bot}(\rho) + \frac12 \tilde r_\mathrm{top}(\rho). \end{aligned}$$

The other integrals can be expanded similarly.

6. Transfer matrix formalism for layered media

For piece-wise continuous media, the solutions of (16b) can be conveniently formulated in a transfer matrix formalism closely related to the Abelés formalism for reflectivity [44]. Here we follow the convention of Chilwell and Hodgkinson [63].

6.1 Single layer

First, we bring the second-order equation $(L - k^2 \nu ^2) u = 0$ into first-order form, setting

(86)$$V = \begin{pmatrix} u \\ k^{{-}1} s^{{-}1} u' \end{pmatrix}$$

so that

(87)$$V' = \mathrm{Q} V = \begin{pmatrix} 0 & k s \\ - k s^{{-}1}(n^2 - \nu^2) & 0 \end{pmatrix} V.$$

Note that $V(z)$ is manifestly continuous. Diagonalizing $\mathrm {Q}$ gives

(88)$$\mathrm{Q}(z) = \mathrm{A}(z) \begin{pmatrix} i k p(z) & 0 \\ 0 & -i k p(z) \\ \end{pmatrix} \mathrm{A}^{{-}1}(z), \quad \mathrm{A} = \begin{pmatrix} 1 & 1\\ \frac{i p}{s} & \frac{-i p}{s} \end{pmatrix},$$

where $p^2 = (n^2 - \nu ^2)$. The sign of $p$ is arbitrary, as long as it is consistent in all 3 matrices. For constant $n$ and $s$, we read off that

(89)$$V(z+d) = \mathrm{A}(z+d) \begin{pmatrix} e^{i k p d} & 0 \\ 0 & e^{{-}i k p d} \\ \end{pmatrix} \mathrm{A}^{{-}1}(z) V(z).$$

This immediately generalizes to piece-wise constant media as follows.

6.2 Multiple layers

Let $N \ge 2$ be the number of layers and let $\underline {n} \in \mathbb {C}^N$ and $\underline {s} \in \mathbb {C}^N$ denote values of $n_l$ and $s_l$ in each layer. Moreover, let $\underline {z} = (z_0, z_1, z_2, \dots, z_{N}) \in \mathbb {R}^{N-1}$ denote the interface positions as shown in Fig. 3. Then

(90)$$V(z_{l+1}) = \mathrm{M}_l(z_{l+1}-z_l) V(z_l),$$

where, with $\mathrm {A}_l = \mathrm {A}(z_l)$,

(91)$$\mathrm{M}_l(d) = \mathrm{A}_l \begin{pmatrix} e^{i p_l d} & 0 \\ 0 & e^{{-}i p_l d} \end{pmatrix} \mathrm{A}^{{-}1}_l = \begin{pmatrix} \cos (p_l k d) & \frac{s_l}{p_l}\sin (p_l k d) \\ -\frac{p_l}{s_l} \sin (p_l k d) & \cos(p_l k d) \end{pmatrix}.$$

The transfer matrix from $z_n$ to $z_m > z_n$ can be seen to be the matrix product

(92)$$\mathrm{M}(z_{m}, z_n) = \mathrm{M}_{m-1}(z_{m}, z_{m-1}) \cdots \mathrm{M}_n(z_{n+1}, z_n),$$

as well as $\mathrm {M}(z_n, z_m) = \mathrm {M}^{-1}(z_m, z_n)$. Note that when allowing arbitrary position arguments, $\mathrm {M}(z, z')$, for fixed $z'$, is a fundamental matrix of (87), since its two columns are linearly independent solutions.

Fig. 3. We consider a layer system with piece-wise constant refractive indices $n_l$ and relative permeabilities and permittivities $s_l$ indicated by the colors. Outside the interval $(z_-, z_+)$, the quantities are constant.

Download Full Size | PDF

The two decaying solutions (29) can be computed as

(93)$$\begin{pmatrix} u_\pm(z)\\ k^{{-}1} s^{{-}1}(z) u'_\pm(z) \end{pmatrix} = \mathrm{M}(z, z_\pm) V_\pm, \quad V_\pm{=} \begin{pmatrix} 1 \\ \pm i p_\pm s^{{-}1}(z_\pm)) \end{pmatrix}.$$

The Wronskian function can be computed via

(94)$$\begin{aligned}W &= k \begin{pmatrix} i s^{{-}1}(z_+) p_+, & -1 \end{pmatrix} \mathrm{M}(z_+, z_-) V_-\\ &= k \begin{pmatrix} i s^{{-}1}(z_-) p_-, & 1 \end{pmatrix} \mathrm{M}(z_-, z_+) V_+. \end{aligned}$$

6.3 Numerical solutions

Having formulated matrix expressions for the solutions $u_\pm (z; \nu )$ and the Wronskian $W(\nu )$, it is straightforward to numerically compute the semi-spectral representation of the Green’s functions $\overleftrightarrow {\mathbf {G}}(q, z, z')$. For the asymptotic expansions of the real-space Green’s functions and to find the resonant modes, we need the complex roots of the Wronskian.

Since $W(\nu )$ is holomorphic up to 4 algebraic and possibly degenerate branch points, Cauchy’s argument principle can be employed to count and localize the complex roots. This approach, commonly known as Cauchy’s integration method or the argument principle method, has been used in numerical waveguide mode finders for many years (see [41] and the references therein).

Since $W(\nu )$ continuously depends on the parameters of the layer system such as layer thickness, refractive indices, and wavenumber, the location of the roots of $W(\nu )$ does, too. As a consequence, the roots can be conveniently tracked. This has been formally proven for this particular problem [39]. A root $\nu _m \in \mathbb {C}$ that is known for a certain configuration can be easily found by local root-searching methods such as Newton’s method for a slightly changed configuration.

We have implemented the transfer matrices in C++ with bindings to the Python language to efficiently compute the semi-spectral Green’s functions as well as the Wronskian function. We use the open source software package cxroots [64] to find all roots within a given complex contour. These are then used to compute the transversal mode profiles as well as the asymptotic expansions of the Green’s functions. The code is available at Ref. [53].

7. Applications

7.1 Grazing incidence

We discuss the scattering of a plane wave illuminating a layer system in grazing incidence. The geometry is translationally symmetric and as such reduces to a one-dimensional problem. By convention, we chose the coordinate system so that $z$ is pointing "down" and the "top" layer extends to $z \to -\infty$ (see Fig. 1(a)). We derive the fields generated by a horizontal sheet of source currents embedded into the top layer, which would produce a plane wave in a homogeneous environment. The fields can be expressed in terms of the semi-spectral representation of the two-dimensional Green’s function evaluated on the real axis. Resonant modes manifest as damped resonances.

Consider an electrical current distribution generating a $y$-polarized plane wave incoming under the angle $\theta \in (0, \pi /2)$ with respect to the $x$-axis,

(95)$$\vec J_\mathrm{e}(x, z) = \delta(z - z_\mathrm{src}) e^{i n_- k \cos \theta x} J_{\mathrm{e},0}(\theta) \hat{y},$$

where $n_-$ is the refractive index in the top layer. This current generates the electric field

(96)$$\begin{aligned}\vec E(x,z) &= i \omega \mu_0 \int \overleftrightarrow{\mathbf{G}}^\mathrm{ls}_\mathrm{e}(x-x', z, z') \vec J_\mathrm{e}(x', z') \mathrm{d} x' \mathrm{d} z'\\ &= i \omega \mu_0 J_{\mathrm{e},0}(\theta) \int g(x-x', z, z_\mathrm{src}) e^{i n_- k \cos \theta x'} \mathrm{d} x' \hat{y}. \end{aligned}$$

The integral takes the form of a Fourier transform and can be rewritten as

(97)$$\vec E(x,z) = \frac{i 2 k p_- E_0}{\mu_-} e^{i n_- k \cos \theta x} g(q= n_- k \cos \theta, z, z_\mathrm{src}) \hat{y},$$

where $E_0 = -\mu _0 \mu _- J_{\mathrm {e},0} c /(2 p_-)$, $p_- = n_- \left | \sin \theta \right |$, and $\mu _-$ is the permeability of the top layer.

In a homogeneous medium we obtain, using (54),

(98)$$\vec E_\mathrm{inc}(x,z) = E_0 e^{i n_- k(x \cos \theta + (z-z_\mathrm{src}) \sin \theta)} \hat{y},$$

which is a plane wave propagating in $\hat {k} = \cos \theta \hat {x} + \sin \theta \hat {z}$. An arbitrary layer system can be readily computed by using (54) to compute $g$ in (97). Substituting (54) and using (32), we bring (97) into a more explicit form,

(99)$$\vec E(x,z) = i 2 k p_- E_0 e^{i n_- k \cos \theta x} \left.\frac{u_-(z_\mathrm{src}; q) u_+(z; q)}{ik p_- u_+(z_-;q) + u_+'(z_-; q)}\right|_{q= n_- k \cos \theta} \hat{y},$$

where we have assumed that $z > z_\mathrm {src}$. For points $z$ in the top layer, we can further simplify the expression. Rewriting $u_+(z)$ using (89), we obtain

(100)$$\vec E(x,z) = \vec E_\mathrm{inc}(x,z) \left \{ 1 + r e^{- 2 i n_- k \sin\theta (z - z_-)} \right\}.$$

with

(101)$$r = \left.\frac{i k p_- u_+(z_-; q) - u_+'(z_-; q)}{i k p_- u_+(z_-; q) + u_+'(z_-; q)}\right|_{q= n_- k \cos \theta}.$$

The second term in (100) represents the reflected part of the field. For completeness, let us remark that the amplitude reflection coefficient $r$ can be readily computed using the transfer matrix formalism,

(102)$$r = \frac{C_2}{C_1}, \quad \begin{pmatrix} C_1 \\ C_2 \end{pmatrix} = \mathrm{A}^{{-}1}(z_-) \mathrm{M}(z_-, z_+) V_+,$$

using the definitions from section 6. This is equivalent to the popular matrix method proposed by Abelés [44].

Figure 4 shows the mode structure and electric field amplitude in an air/Pt(2.6 nm)/C(16 nm)/Fe(0.6 nm)/C(16 nm)/Pt layer system at a photon energy of 14.4 keV. The photon energy $\mathcal {E}$ is converted to the vacuum wavenumber via $k = \mathcal {E} / (c \hbar )$. The layer structure corresponds to a cavity used in an experiment reported in [9]. Here, we are interested in non-resonant scattering applications, and ignore the resonant nuclear contribution of the Mössbauer isotope $^{57}$Fe to the refractive index. The refractive indices were taken from [65]. Since the refractive index of carbon is lower than that of air, the system does not permit guided modes but only negative and double leaky modes, as can be seen in panel a). In panel b), the squared modulus of the 2d Green’s function is plotted as a function of $\nu = q/k$ on the real axis, which corresponds to the incidence angle via $\nu = \cos \theta$. The Green’s function $\left | g(k \nu ) \right |^2$ is proportional to the local intensity $\left | E \right |^2$ based on (97). It shows pronounced resonances around the real values of the resonant mode indices, $\nu \sim \textrm{Re}\{ \nu _m\}$. The vertical profiles in b) for real $\nu \sim \textrm{Re}\{ \nu _m\}$, closely resemble the vertical profiles at the complex resonances $\nu _m$, shown in panel c), up to the divergence for $z < 0$. Hence, one can get a good understanding of field within the waveguide from just the resonant modes.

Fig. 4. Mode structure and field amplitude of the layer system shown in (c) at 14.4 keV photon energy. (a) Mode indices in the complex plane. (b) Squared modulus of the 2d Green’s function, which corresponds to the electric field as a function of height and incidence angle $\nu = \cos \theta$. The dashed lines indicate the interfaces between platinum and carbon. (c) Layer system and transversal profiles of the negative leaky modes (squared modulus). The profiles in (c) correspond to the mode indices in (a) of matching colors.

Download Full Size | PDF

7.2 Waveguide propagation

Let us now discuss the propagation of x-rays in a waveguide layer system. To that end, we derive the fields generated by a vertical sheet of source currents that would generate a plane wave in a homogeneous medium. This situation closely corresponds to a plane wave that is coupled into the front face of a waveguide, as long as back-reflection can be neglected. Having a plane of origin perpendicular to the layer interfaces breaks the translational symmetry and makes this problem inherently two-dimensional. The fields in the waveguide can be computed using the two-dimensional Green’s function, and we will see that the resonant modes emerge as an asymptotic approximation for longer propagation distances.

Consider a vertical sheet of electric currents at $x = x_\mathrm {src}$ with a linear phase shift,

(103)$$\vec J_\mathrm{e}(x, z) = \delta(x - x_\mathrm{src}) e^{i n_0 \sin \theta z} J_{\mathrm{e},0}(\theta) \hat{y}.$$

This current generates the electric field

(104)$$\begin{aligned}\vec E(x,z) &= i \omega \mu_0 \int \overleftrightarrow{\mathbf{G}}^\mathrm{ls}_\mathrm{e}(x-x', z, z') \vec J_\mathrm{e}(x', z') \mathrm{d} z' \mathrm{d} z'\\ &= i \omega \mu_0 J_{\mathrm{e},0}(\theta) \int g(x-x_\mathrm{src}, z, z_\mathrm{src}) e^{i n_0 \sin \theta z'} \mathrm{d} x' \hat{y}. \end{aligned}$$

For brevity, we set $x_\mathrm {src} = 0$.

In a homogeneous medium with refractive index $n_\mathrm {hom}$, we get

(105)$$\vec E_\mathrm{inc}(x,z) = E_0 e^{i n_\mathrm{hom} k \hat{k} \cdot \vec r}$$

with $\hat {k} = \cos \theta \hat {x} + \sin \theta \hat {z}$ and $E_0 = - \mu _0 J_{\mathrm {e},0} \omega / (2 k n_\mathrm {hom} \cos \theta )$. For stratified media, we use (63) to obtain

(106)$$\vec E(x,z) = E_0 \sum_m C_m e^{i k \nu_m x} u_m(z) \hat{y} + \text{non-resonant},$$

where

(107)$$C_m = \int \mathrm{d} z' u_m(z') e^{i n_0 z' \sin \theta} \frac{n_\mathrm{hom} \cos \theta }{\nu_m}$$

is a complex coefficient that quantifies how the incoming field couples into the modes. The sum goes over all guided modes and the non-resonant part is the integral over the radiative modes. Note that for this infinitely extended source sheet, leaky modes are of little use due to their divergence. The coupling coefficients (107) are determined by the integral over the mode profile $u_m$.

All relevant information is contained in the mode indices $\nu _m$ and the mode profiles $u_m(z)$. A single mode decays like

(108)$$\left| \vec E(x,z)_m \right|^2 = \left| E_0 \right|^2 \left| C_m \right|^2 \left| u_m(z) \right|^2 e^{- 2 k \operatorname{Im} \nu_m x},$$

so that the $1/e$-length is given by $\lambda /\left(4 \pi \operatorname{Im}\left\{v_m\right\}\right)$. Two or more contributing modes interfere with period $\lambda / \operatorname{Re}\left\{v_m-v_{m^{\prime}}\right\}$.

Figure 5(a) shows the mode indices for a Mo/B$_4$C(20 nm)/Mo waveguide and a photon energy of 13.8 keV, (b) shows the corresponding transversal mode profiles of the 2 guided modes and the first 3 double leaky modes. Here, the leaky mode indices have significantly larger imaginary parts compared to the guided and radiative modes, where the guided mode indices have the smallest imaginary parts because of the weak absorption of B$_4$C compared to molybdenum.

Fig. 5. Mode structure and mode transmission of a Mo/B$_4$C/Mo waveguide at 13.8 keV photon energy. (a) Mode and material indices, including guided modes (filled circles), double leaky modes (empty circles) and radiative modes (dotted line). (b) Transversal mode profiles (squared modulus) corresponding to the mode indices in a. (c) Mode attenuation lengths ($1/e$) for the waveguide as a function of the B$_4$C guiding layer thickness. The drawn through and dashed lines correspond to guided and double leaky modes, respectively.

Download Full Size | PDF

Figure 5(c) shows the attenuation lengths of the resonant modes as a function of the B$_4$C core layer thickness. It demonstrates that the number of resonant modes stays invariant, while leaky modes turn into guided modes for increasing thickness.

Figure 6 shows the electric field modulus inside the same waveguide, irradiated by a $y$-polarized plane wave under normal incidence. The top row shows only the guided mode contributions computed with (106) while the bottom row shows the full solution, obtained in a finite-difference simulation of the paraxial Helmholtz equation (see [37] for details). Due to the normal incidence and symmetry, only the first guided mode contributes. The full field shows some interference in the first few 100 µm, which is caused by the radiative modes interfering with the guided modes (compare [40]). This is more evident in the Fourier-transformed fields. The full numerical solution (bottom right) shows background contributions that peak around the refractive index of the cladding material, molybdenum, in addition to the resonance of the guided mode (top right). Evidently, the non-resonant contributions decay on significantly shorter length scales, so that the field after $\sim {300}\;{\upmu}\textrm{m}$ is accurately described by only the guided modes. This highlights the validity of the guided-mode approximation as an asymptotic expansion of the field for large propagation distances, which is possible because of the large material absorption.

Fig. 6. Electric field in the waveguide for an incoming plane wave at 0° incidence angle computed from only the guided modes (top) and the full field (bottom), as a function of spatial coordinates (left) and Fourier-transformed with respect to $x$ (right).

Download Full Size | PDF

So far we have considered strictly time-harmonic fields. X-ray free electron lasers produce pulses of a few femtoseconds, so the questions arises if waveguides introduce dispersion. The bandwidth of hard x-rays produced by synchrotron or x-ray free electron laser sources is generally in the order of 1 eV. The change of refractive indices in this energy bandwidth is generally negligible, as long as the photon energy is far from absorption edges and resonances. We can thus approximate the refractive index to be independent of frequency. The wavenumber $k = \omega /c$, however, enters $W(\nu )$ and in principle causes a frequency-dependence of the mode indices $\nu _m \equiv \nu _m(\omega )$ and $u_m(z) \equiv u_m(z, \omega )$. Consequently, there is in principle intra-mode dispersion due to the frequency-dependence of $\nu _m(\omega )$ as well as inter-mode dispersion due to the difference in $\nu _m$. Since $k$ only changes on the order of the relative bandwidth, which is typically 1×10⁻⁴, the change in $u_m$ and $\nu _m$ (compare (41)), and thus also intra-mode dispersion, are typically negligible. Neglecting the time-dependence of $u_m$ and $\nu _m$ and splitting the real and imaginary part of $\nu _m$, (106) becomes

(109)$$\vec E(x,z,t) \approx \hat{y} \sum_m E_0(t - \operatorname{Re} \nu_m x /c) C_m e^{{-}2 k_0 \operatorname{Im} \nu_m x} u_m(z),$$

where $k_0 = \omega _0 / c$ is the central wave number. Pulses of length $\Delta t$ thus separate on length scales of $\Delta t c / \textrm{Re}\{\nu _m - \nu _{m'}\}$. For $\Delta t = {1}\;\textrm{fs}$ and $\textrm{Re}\{\nu _m - \nu _{m'}\} {\sim} 10^{-5}$ this dispersion length is 3 cm. Conversely, for longer pulses or shorter waveguides, dispersion can be neglected (compare [31]).

7.3 Dipole emission into waveguide modes

Finally, we discuss the fields produced by a single dipole buried in a waveguide layer structure. This corresponds to, for example, fluorescence or characteristic x-ray emission into a waveguide, as it was experimentally observed recently [21]. The fields can be described using the three-dimensional Green’s function, which asymptotically reduces to its guided modes.

For simplicity, consider a vertically oriented dipole, $\vec p = \left | p \right | \hat {z}$, so that the problem is rotationally symmetric. Using (10), the electric field produced by the dipole can be written in cylindrical coordinates as

(110)$$\vec E(\rho, \phi, z) = \omega^2 \mu_0 \overleftrightarrow{\mathbf{G}}_\mathrm{e}(\rho, \phi, z, z') \cdot \vec p.$$

The field is then TM-like in the sense that the electric field is oriented in the $\hat { \rho }$-$\hat {z}$ plane with $E_\rho = \omega ^2 \mu \mu _0 G_{xz}(\rho, 0, z, z') |p|$ and $E_z = \omega ^2 \mu \mu _0 G_{zz}(\rho, 0, z,z') |p|$. Using (85) and neglecting the non-resonant contributions, we have

(111)$$E_\rho(\rho, z) = \frac{\mu_0 \omega^2 \left| p \right|}{4\epsilon(z) \epsilon(z') k} \sum_m \nu_{\mathrm{m},m} H_1^{(1)}(k \nu_{\mathrm{m},m} \rho) u_{\mathrm{m},m}(z') \partial_z u_{\mathrm{m},m}(z)$$

(112)$$E_z(\rho, z) = i \frac{\mu_0 \omega^2 \left| p \right|}{4\epsilon(z) \epsilon(z')} \sum_m \nu_{\mathrm{m},m}^2 H_0^{(1)}(k \nu_{\mathrm{m},m} \rho) u_{\mathrm{m},m}(z') u_{\mathrm{m},m}(z).$$

Note that $\rho$ only appears in the Hankel functions. They have the asymptotic expressions

(113)$$H^{(1)}_n(\xi) \sim \sqrt{\frac{2}{\pi \xi}} e^{{-}i(1+2n) \pi /4} e^{i \xi}$$

for $\left | \xi \right | \to \infty$. Although both field components decay similarly, the global magnitude of the radial component is generally a lot smaller. For $k \rho \gg 1$, the absolute square of an individual mode becomes

(114)$$\left| E_{z,m}(\rho, z) \right|^2 \approx \frac{\mu_0 \omega^4 \left| p \right|^2 }{8 \pi} \frac{\left| u_m(z) \right|^2}{\left| \epsilon(z) \right|^2} \frac{\left| u_m(z') \right|^2}{\left| \epsilon(z') \right|^2} \frac{\exp\{{-}2k \rho \operatorname{Im} \nu_m \}}{k \rho},$$

where we have approximated $\nu _m \approx 1$.

We can read off several characteristics. The absolute square decays hyperbolically due to the planar geometry and exponentially due to material absorption, which is quantified by the complex effective index. The coupling into the mode scales with the mode function $\left | u_m(z') \right |^2$. If more than one mode is populated, the fields interfere on the length scale $\rho \sim \lambda / \operatorname{Re}\left\{v_m-v_{m^{\prime}}\right\}$.

Let us now compute the total power flow into a given mode $m$. The radial component of the Poynting vector is given by

(115)$$\hat{ \rho} \cdot \vec S(\rho, z) = \frac{-i}{2 \omega \mu_0 \bar \mu} E_z (\partial_\rho \bar E_z -\partial_z \bar E_\rho).$$

Inserting (112) and neglecting the second term, we obtain

(116)$$\frac{\omega^3 \left| p \right|^2 \mu_0 \mu(z)}{32 } \left| \nu_m \right|^4 k \bar \nu_m \left| H_0^{(1)}(k \nu_m \rho) \right|^2 \frac{\left| u_m(z') \right|^2}{\left| \epsilon(z') \right|^2} \frac{\left| u_m(z) \right|^2}{\left| \epsilon(z) \right|^2}.$$

Integrating over the cylinder surface and inserting (113) gives

(117)$$P_m(\rho) \approx \frac{\omega^3 \left| p \right|^2 \mu_0}{8} \left| \nu_m \right|^3 \bar\nu_m \exp\{{-}2 k \rho \operatorname{Im} \nu_m\} \frac{\left| u_m(z') \right|^2}{\left| \epsilon(z') \right|^2} \int_{-\infty}^\infty \frac{\left| u_m(z) \right|^2 \mu(z)}{\left| \epsilon(z) \right|^2} \mathrm{d} z$$

The remaining integral evaluates to a number close to unity. We obtain the estimate for the total radiated power

(118)$$P_m(\rho) \sim \frac{\omega \left| p \right|^2 k^3}{4 \pi \epsilon_0} \exp\{{-}2 k \rho \operatorname{Im} \nu_m\} \frac{\lambda \left| u_m(z') \right|^2}{4}.$$

Compared to the total radiated power of a dipole in vacuum [66],

(119)$$P_\mathrm{hom} = \frac{\omega \left| p \right|^2 k^3}{12 \pi \epsilon_0},$$

we notice the scaling with the coupling coefficient $\left | u_m(z') \right |^2 \lambda$. Since $u_m$ is bi-normalized, it takes values in the order of the reciprocal mode width $l$. We obtain the simple result that the radiated power scales effectively with the geometric factor $\sim \lambda / l$. For hard x-rays, $l$ is typically at least 2 orders of magnitude larger than the wavelength, so that the coupling efficiency is in the order of 0.01.

As a numerical example, consider a Ni/C(24.5 nm)/Fe(1 nm)/C(24.5 nm)/Ni layer system, in which one of the iron atoms in the central layer undergoes a radiative decay, emitting a 6.4 keV K$\alpha$ fluorescence photon. This structure corresponds to a sample used in experiments by Vassholz and Salditt [21]. We model the steady-state case by a classical electric dipole. Figure 7(a) shows the transversal guided mode profiles in the layer structure. The system admits 4 guided modes but only two of them couple to the dipole. Panel (b) shows the transversal field computed with (112) up to a global scaling and compensated by the geometrical $1/\rho$-decay. It shows pronounced interference between the two contributing modes as well as an overall attenuation due to photoelectric absorption. This interference has been experimentally observed by Vassholz and Salditt (see Fig. 4 in [21]). In their original interpretation, however, the oscillations would be caused by a change of the emission process itself due to the position of the emitter within the waveguide. Our theory clearly indicates that the oscillations can be caused by the interference of the two populated guided modes alone, without assuming a modification of the emission process such as the Purcell effect.

Fig. 7. Fields produced by an electric dipole buried in a waveguide at the iron K$\alpha$ fluorescence energy 6.4 keV. (a) Layer structure and transversal profiles (squared modulus) of the guided modes. A dipole in the central layer couples to the 1st and 3rd (solid lines) but not to the 2nd and 4th (dashed lines) guided modes. (b) Electric field amplitude produced by an iron atom in the central layer, considering only the two contributing guided modes.

Download Full Size | PDF

Next, we analyze how the waveguide thickness affects some emission characteristics. Figure 8 shows the $1/e$ attenuation length and coupling coefficient of the guided modes as a function of overall waveguide thickness. Increasing the thickness converts more leaky modes into guided modes and increases the attenuation length of each individual mode. Less desirably, increasing the thickness also decreases the coupling strength of each individual mode slightly on a linear scale. The first guided mode has a pronounced optimum between 5 nm to 10 nm, below which the coupling coefficient rapidly drops to 0. Since the coupling coefficients are approximately inversely proportional to the square root of the mode spread, the optimum coupling corresponds to the minimum mode spread. Smaller waveguides cannot effectively contain the modes due to the weak index contrast.

Fig. 8. Attenuation length (a) and mode coupling coefficient (b) as a function of waveguide thickness for an electric dipole buried in a waveguide at the iron K$\alpha$ fluorescence energy 6.4 keV. The layer structure is Ni/C($x$)/Fe(1 nm)/C($x$)/Ni and the core thickness $2 x+1 \mathrm{~nm}$ is varied.

Download Full Size | PDF

8. Conclusions

In summary, we have presented a self-contained nano-optical theory of planar x-ray waveguides combined with numerical code, that allows for the modelling of fields in arbitrary layer structures while fully including material absorption. As one of the main results, we have use a novel generalised bi-orthogonality relation to obtain asymptotic expansions into resonant modes that are particularly useful for x-ray wavelengths due to the presence of significant material absorption. Our theory allows to model, in a unified framework, both the "classical" geometries of x-ray waveguides, front-coupling and evanescent coupling geometries with their main figures of merit, as well as atom-light interactions with buried dipoles. In addition, we have discussed how interference patterns emerge in the fields caused by buried emitters, offering new interpretations of recently published experiments.

In more technical terms, we have shown how the solutions of the frequency-space vector wave Eqs. (5) of lossy stratified media, which are describing the electromagnetic fields in x-ray waveguides, can be decomposed into resonant and non-resonant modes, solutions of the transversal eigenvalue problem, and we have derived asymptotic expansions of the corresponding dyadic Green’s functions in terms of the resonant modes. The modes of the free solution naturally reappear in the semi-spectral representation of the Green’s function, where the resonant modes manifest as pole singularities.

The presence of material losses makes the underlying equations inherently non-hermitian and causes the resonant modes to violate the usual orthogonality property, which is replaced by a bi-orthogonality property (28). Due to the dielectric boundary layers, the mode cut-off is soft, and in general an infinite number of modes are supported. However, material losses in many cases rapidly attenuate all but the lowest order resonant modes so that the asymptotic mode expansions give accurate results on relevant length scales.

The resonant mode expansion is similar to but distinct from so-called quasi normal mode expansions (see for example [57,67]). While quasi-normal modes usually refer to complex-frequency resonances of the time-dependent solution, the resonant modes discussed in this paper correspond resonances of the monochromatic frequency-space solution with respect to the complex longitudinal wavenumber $q$. In particular, we note that the 3-dimensional vector wave equations for stratified media, as they were discussed here, do not have resonant frequencies, because the system is translationally invariant.

We have applied the theory to important experimental geometries, including irradiation by a plane wave in grazing incidence, propagation of a plane wave inside a waveguide, and emission into waveguide modes from a dipole embedded inside the waveguide. The unified theory highlights the similarities in the underlying mode structure of the 3 geometries but also elucidates important differences due to the different dimensionalities, symmetries, and boundary conditions.

The results presented in this paper have further applications. First, the theory together with our numerical solvers allows us to efficiently optimize layer designs for specific experimental needs. Important figures of merit such as attenuation length and coupling efficiencies are immediately obtainable from the the mode indices and profiles. This facilitates the computer-aided design of x-ray waveguides for imaging, where previously one mainly relied on finite-difference propagation simulations, which are relatively complex and provide little structural insight. Another example is the design and optimization of structured multilayer x-ray sources as proposed in Ref. [21].

Second, the dyadic Green’s functions and their asymptotic mode expansions can be used for quantum theories of, for example, coupled quantum systems in x-ray waveguides, since the classical Green’s function enters the quantization of the electromagnetic field in the framework of macroscopic quantum electrodynamics [46–49]. Our analytic theory, in particular the resonant mode expansion, will help to clarify the structure of the electromagnetic fields in x-ray waveguides. In addition, our numerical tools will complement current efforts [68] to design waveguide structures for quantum-optical experiments.

A. Integration by parts for the modified inner product

We will demonstrate that for the class of function that are of the form

(120)$$u(z) = u_0 e^{i k_\pm z}$$

in the top (+) and bottom (-) layer, that our modified inner product obeys the following integration by parts law,

(121)$$\left\langle u',v\right\rangle_\partial + \left\langle u, v'\right\rangle_\partial = 0.$$

Using the definition (24), we have

(122)$$\left\langle u',v\right\rangle_\partial + \left\langle u, v'\right\rangle_\partial = \int_{z_-}^{z_+} \bar{u}'(z) v(z) +\bar{u}(z)v'(z) \mathrm{d}{z} -\left[ \frac{\bar{u}' v +\bar{u}v'}{\bar{u}'/u + v'/v} \right]_{z_-}^{z_+},$$

where we have used the property that evaluated at $z_\pm$ we have

(123)$$u^{\prime\prime}/u' = u'/u, \quad v^{\prime\prime}/v' = v'/v.$$

We may use ordinary integration by parts to evaluate the integral, giving

(124)$$\left\langle u',v\right\rangle_\partial + \left\langle u, v'\right\rangle_\partial = \left[ \bar{u}v - \frac{\bar{u}' v +\bar{u}v'}{\bar{u}'/\bar{u} + v'/v} \right]_{z_-}^{z_+} = 0,$$

since the right hand term can be simplified to

(125)$$\frac{\bar{u}'v + \bar{u}v'}{\bar{u}'/\bar{u} + v'/v} = \bar{u}v\frac{\bar{u}'/\bar{u} + v'/v}{\bar{u}'/\bar{u} + v'/v} = \bar{u}v.$$

B. Derivation of mode dispersion

We derive identity (41) for the resonant modes.

Beginning with the eigenvalue equation

(126)$$L u = \lambda u,$$

we take a frequency derivative to obtain

(127)$$\frac{\partial L}{\partial \omega} u + L \frac{\partial u}{\partial \omega} = \frac{\partial \lambda}{\partial \omega} u + \lambda \frac{\partial u }{\partial \omega}.$$

As $u$ can be written as an exponential function in the top and bottom layers, so can $\frac {\partial L}{\partial \omega } u$ and $\frac {\partial u}{\partial \omega }$, and our modified inner product remains valid to use. Thus, we project onto the adjoint eigenfunction $\overline {s^{-1} u}$ with the modified inner product, to obtain

(128)$$\left\langle \overline{s^{{-}1} u}, \frac{\partial L}{\partial \omega} u\right\rangle_\partial + \left\langle \overline{s^{{-}1} u}, L \frac{\partial u}{\partial \omega}\right\rangle_\partial = \frac{\partial \lambda}{\partial \omega} + \lambda \left\langle \overline{s^{{-}1}u},\frac{\partial u }{\partial \omega}\right\rangle_\partial.$$

We can then note that as we have the differentiation by parts identity (25), $\overline {s^{-1}u}$ acts an adjoint eigenfunction with respect to our modified inner product, and thus we have

(129)$$\left\langle \overline{s^{{-}1} u}, L \frac{\partial u}{\partial \omega}\right\rangle_\partial = \lambda \left\langle \overline{s^{{-}1} u}, \frac{\partial u}{\partial \omega}\right\rangle_\partial.$$

We obtain the result

(130)$$\frac{\partial \lambda}{\partial \omega} = \left\langle \frac{\partial L}{\partial \omega}\right\rangle_u,$$

where for an operator $A$ we denote the ‘expectation value’ with respect to $u$ as

(131)$$\left\langle A\right\rangle_u \equiv \left\langle \overline{s^{{-}1} u}, Au\right\rangle_\partial.$$

Evaluating the frequency derivative of $L$, we have

(132)$$\frac{\partial L}{\partial \omega} = \frac{\partial s}{\partial \omega}\partial_z s^{{-}1}\partial_z -s\partial_z s^{{-}2}\frac{\partial s}{\partial \omega}\partial_z + \frac{2\omega^2}{c^2}\left( \frac{n^2}{\omega} + \frac{n}{\omega}\frac{\partial n}{\partial \omega} \right).$$

The rightmost terms of (132) are straightforward to compute with respect to the expectation value, but care must be taken when evaluating the first two terms. Taking an expecation value with respect to $u$ gives for the first term

(133)$$\left\langle \frac{\partial s}{\partial \omega}\partial_z s^{{-}1}\partial_z\right\rangle_u = \left\langle s^{{-}1}\frac{\partial s}{\partial \omega} (L - \frac{\omega^2}{c^2}n^2)\right\rangle_u = \left\langle s^{{-}1}\frac{\partial s}{\partial \omega} (\lambda - \frac{\omega^2}{c^2}n^2)\right\rangle_u.$$

For the second term we have

(134)$$\left\langle s\partial_z s^{{-}2}\frac{\partial s}{\partial \omega}\partial_z \right\rangle_u = \left\langle \overline{s^{{-}1} u}, s\partial_z s^{{-}2}\frac{\partial s}{\partial \omega}\partial_z u \right\rangle_\partial = \left\langle \overline{u}, \partial_z s^{{-}2}\frac{\partial s}{\partial \omega}\partial_z u \right\rangle_\partial,$$

where we have used the antilinearity of the inner product. We can then use integration by parts to obtain

(135)$$\left\langle s\partial_z s^{{-}2}\frac{\partial s}{\partial \omega}\partial_z \right\rangle_u ={-}\left\langle \overline{s^{{-}1}\partial_z u}, s^{{-}1}\frac{\partial s}{\partial \omega}\partial_z u \right\rangle_\partial ={-}\left\langle s^{{-}1}\frac{\partial s}{\partial \omega}\right\rangle_{\partial_z u}.$$

Putting it all together, we have

(136)$$\frac{\partial \lambda}{\partial \omega} = \left\langle s^{{-}1}\frac{\partial s}{\partial \omega}(\lambda -\frac{\omega^2}{c^2}n^2)\right\rangle_u + \left\langle s^{{-}1}\frac{\partial s}{\partial \omega}\right\rangle_{\partial_z u} + \frac{2\omega^2}{c^2} \left( \left\langle \frac{n^2}{\omega}\right\rangle_u + \left\langle n\frac{\partial n }{\partial \omega}\right\rangle_u \right).$$

Finally, we have

(137)$$\lambda = \frac{\omega^2}{c^2}\nu^2,$$

(138)$$\frac{\partial \lambda}{\partial \omega} = \frac{2\omega^2}{c^2}\left( \frac{\nu^2}{\omega} +\nu \frac{\partial \nu}{\partial \omega} \right),$$

giving us for the mode dispersion

(139)$$\nu \frac{\partial \nu}{\partial \omega} = \frac{1}{2}\left\langle \frac{\partial \ln(s)}{\partial \omega} (\nu^2 - n^2) \right\rangle_u + \frac{c^2}{2\omega^2} \left\langle \frac{\partial \ln(s)}{\partial \omega}\right\rangle_{\partial_z u} + \left\langle \frac{n^2 - \nu^2}{\omega}\right\rangle_u + \left\langle n \frac{\partial n}{\partial \omega}\right\rangle_u.$$

C. Derivation of the Wronskian derivative

We derive the identity (60). The derivation is inspired by a similar derivation for quasi-normal modes of the vector wave equations [57]. We have, using (15), that

(140)$$k^2 (\nu^2 - \nu_j^2) \int_{z_-}^{z_+} s^{{-}1} u_-(\nu_j) u_+(\nu) \mathrm{d} z = \left[ s^{{-}1} \left( u_-(\nu_j) u_+'(\nu) - u_-'(\nu_j) u_+(\nu)\right)\right]_{z=z_-}^{z_+}.$$

Using that $u_-'(\nu, z_-) = -i k p_-(\nu ) u_-(\nu, z_-)$ and $u_+'(\nu, z_+) = i k p_+(\nu ) u_+(\nu, z_+)$, as well as $u_-(\nu _j, z) \propto u_+(\nu _j, z)$, the right-hand side of (140) becomes

(141)$$\begin{array}{r} ik \left[ s^{{-}1} \left( p_+(\nu) - p_+(\nu_j) \right) u_-(\nu_j) u_+(\nu) \right]_{z=z_+}\\ - \left[s^{{-}1} \left( u_-(\nu_j) u_+'(\nu) - u_-' (\nu_j) u_+(\nu) \right)\right]_{z=z_-} \end{array}$$

Taking the derivative with respect to $\nu$ at $\nu = \nu _j$, the left-hand side of (140) becomes

(142)$$2 \nu_j k^2 \int_{z_-}^{z_+} s^{{-}1} u_-(\nu_j) u_+(\nu_j) \mathrm{d} z$$

and (141) evaluates to

(143)$$ik \left[ s^{{-}1} \frac{\partial p_+}{\partial \nu} u_- u_+ \right]_{z_+} - \partial_\nu W(\nu_j) - \left[s^{{-}1} \left( \frac{\partial u_-'}{\partial \nu} u_+{-} \frac{\partial u_-}{\partial \nu} u_+'\right)\right]_{z_-},$$

where we have substituted the expression for the Wronskian function (31). To simplify the last term in (143), we use that $u_-(z_-) \equiv 1$ and $u_-'(z_-) = - i k p_- u_-(z_-')$, which implies

(144)$$\left[s^{{-}1} \left( \frac{\partial u_-'}{\partial \nu} u_+{-} \frac{\partial u_-}{\partial \nu} u_+'\right)\right]_{z_-} ={-}ik \left[ s^{{-}1} \frac{\partial p_-}{\partial \nu} u_- u_+ \right]_{z_-}$$

Solving for the Wronskian function yields

(145)$$\begin{aligned}\partial_\nu W(\nu_j) &={-} 2 \nu_j k^2 \int_{z_-}^{z_+} s^{{-}1} u_- u_+ \mathrm{d} z\\ &+ ik \left[ s^{{-}1} \frac{\partial p_+}{\partial \nu} u_- u_+ \right]_{z_+} + ik \left[ s^{{-}1} \frac{\partial p_-}{\partial \nu} u_- u_+ \right]_{z_-}. \end{aligned}$$

Finally, for both the decaying choice of branch cut, $p = p_\mathrm {dec}$, and the vertical choice of branch cut, $p = p_\mathrm {vert}$, we get $\partial _\nu p_\pm = - \nu /p_\pm$, so that

(146)$$\begin{aligned}\partial_q W(\nu_j) &={-} 2 \nu_j k \int_{z_-}^{z_+} s^{{-}1} u_- u_+ \mathrm{d} z - i \nu_j \left[ \frac{u_- u_+}{s p_+} \right]_{z_+} - i \nu_j \left[ \frac{u_- u_+}{s p_-} \right]_{z_-}\\ &={-} 2 \nu_j k \left\langle \overline{s^{{-}1}u}_-, u_+\right\rangle_\partial. \end{aligned}$$

D. Computing the scattered field

A common problem is to evaluate the total field in the wave-guide given the specification of an incident beam in free space. This can be done using virtual sources.

The fields can be expressed in terms of the electric and magnetic currents. Substituting $\vec F$ in (7) and integrating by parts yields for the electric field,

(147)$$\vec E(\vec r) = \frac{i \omega \epsilon_0}{c^2}\int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') \cdot \vec J_\mathrm{e}(\vec r') - \int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') \times \overleftarrow{\nabla} \cdot \mu^{{-}1}(\vec r') \vec J_\mathrm{m}(\vec r'),$$

where we have assumed the boundary terms to vanish. Taking the curl of (9b) and inserting (8), it follows that

(148)$$\overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') \times \overleftarrow{\nabla}' \mu(\vec r')^{{-}1} = \epsilon(\vec r)^{{-}1} \nabla \times \overleftrightarrow{\mathbf{G}}_\mathrm{m}(\vec r, \vec r').$$

Equation (147) then becomes

(149)$$\vec E(\vec r) = \frac{i \omega \epsilon_0}{c^2} \int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_\mathrm{e}(\vec r, \vec r') \cdot \vec J_\mathrm{e}(\vec r') - \frac{1}{\epsilon(\vec r)} \nabla \times \int \mathrm{d}^3 \vec r' \overleftrightarrow{\mathbf{G}}_\mathrm{m}(\vec r, \vec r') \cdot \vec J_\mathrm{m}(\vec r').$$

A similar relation can be derived for the magnetic field.

Let us suppose that the free-space incident beam $\vec {E}_0$ is specified within some connected volume $V$ with boundary $\partial V$. We may then use Stratton-Chu equations to relate the electric field at a point $\vec {r}$ in the interior of $V$ to the electric and magnetic fields on the boundary via (see [28], ch. 1.4)

(150)$$\vec{E}_0(\vec{r}) = \int_{\partial V} \mathrm{d}^2 s \bigg( \frac{i \omega \epsilon_0}{c^2}\overleftrightarrow{\mathbf{G}}^{(0)}(\vec{r},\vec{s}) \cdot \vec{H}_0(\vec{s})\times \hat{n}(\vec{s}) - \overleftrightarrow{\mathbf{G}}^{(0)}(\vec{r},\vec{s}) \cdot \vec{E}_0(\vec{s}) \times \hat{n}(\vec{s}) \bigg),$$

where $\overleftrightarrow {\mathbf {G}}^{(0)}$ denotes a free-space Green’s function. Comparing with (149), we can see that this is equivalent to the field emitted by virtual electric and magnetic currents given by

(151)$$\vec{J}_{\mathrm{e},0}(\vec{r}) = \vec{H}_0(\vec{r}) \times \hat{S}(\vec{r})$$

(152)$$\vec{J}_{\mathrm{m},0}(\vec{r}) = \vec{E}_0(\vec{r}) \times \hat{S}(\vec{r}),$$

where $\vec {S}$ is the surface normal distribution for $\partial V$, equal to the product of the indicator distribution and the field of normal vectors of $\partial V$.

This can be used to obtain the scattered field as follows: suppose $V$ completely encloses all scatterers, and that the Green’s function taking into account the scatterers is known. Then the scattered field can be determined in terms of the incident field by treating the free-space field values along $\partial V$ as virtual sources via (151) and (152), and propagating their field using the full Green’s function, including the effects of scattering. Using (149), it is then straightforward to show that the total field $\vec {E}$, including the scattered field, is given in terms of the incident field $\vec {E}_0$, $\vec {B}_0$ by

(153)$$\begin{aligned} \vec{E}(\vec{r},\omega) &= \frac{i \omega \epsilon_0}{c^2} \int_{\partial V} \mathrm{d}^2{s} \overleftrightarrow{\mathbf{G}}^E(\vec{r},\vec{s},\omega) \cdot \vec{H}_0(\vec{s},\omega)\times \hat{n}(\vec{s})\\ &- \frac{1}{\epsilon(\vec r)} \nabla \times \int_{\partial V} \mathrm{d}^2{s} \overleftrightarrow{\mathbf{G}}^H(\vec{r},\vec{s},\omega) \cdot \vec{E}_0(\vec{s},\omega) \times \hat{n}(\vec{s}) , \end{aligned}$$

where $\overleftrightarrow {\mathbf {G}}$ refers to the full Green’s function.

Funding

Max Planck School of Photonics; Deutsche Forschungsgemeinschaft (432680300 - SFB 1456/C03, 429529648 - TRR 306/C04).

Acknowledgements

The authors have had fruitful discussions with Tim Salditt, Ralf Röhlsberger, Markus Osterhoff, and Adriana Pálffy during the preparation of this manuscript and are very grateful for their support.

Disclosures

The authors declare no conflict of interest.

Data Availability

Data underlying the results presented in this paper are available in Ref. [53].

References

1. E. Spiller and A. Segmüller, “Propagation of x rays in waveguides,” Appl. Phys. Lett. 24(2), 60–61 (1974). [CrossRef]

2. Y. P. Feng, S. K. Sinha, H. W. Deckman, J. B. Hastings, and D. P. Siddons, “X-ray flux enhancement in thin-film waveguides using resonant beam couplers,” Phys. Rev. Lett. 71(4), 537–540 (1993). [CrossRef]

3. M. J. Zwanenburg, J. F. Peters, J. H. H. Bongaerts, S. A. de Vries, D. L. Abernathy, and J. F. van der Veen, “Coherent propagation of x rays in a planar waveguide with a tunable air gap,” Phys. Rev. Lett. 82(8), 1696–1699 (1999). [CrossRef]

4. F. Pfeiffer, C. David, M. Burghammer, C. Riekel, and T. Salditt, “Two-dimensional x-ray waveguides and point sources,” Science 297(5579), 230–234 (2002). [CrossRef]

5. K. Giewekemeyer, H. Neubauer, S. Kalbfleisch, S. P. Krüger, and T. Salditt, “Holographic and diffractive x-ray imaging using waveguides as quasi-point sources,” New J. Phys. 12(3), 035008 (2010). [CrossRef]

6. T. Salditt, M. Osterhoff, M. Krenkel, R. N. Wilke, M. Priebe, M. Bartels, S. Kalbfleisch, and M. Sprung, “Compound focusing mirror and x-ray waveguide optics for coherent imaging and nano-diffraction,” J. Synchrotron Radiat. 22(4), 867–878 (2015). [CrossRef]

7. S. Hoffmann-Urlaub and T. Salditt, “Miniaturized beamsplitters realized by x-ray waveguides,” Acta Crystallogr. Sect. A Found. Adv. 72(5), 515–522 (2016). [CrossRef]

8. 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(20), 203902 (2015). [CrossRef]

9. 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]

10. J. Haber, J. Gollwitzer, S. Francoual, M. Tolkiehn, J. Strempfer, and R. Röhlsberger, “Spectral control of an x-ray L -edge transition via a thin-film cavity,” Phys. Rev. Lett. 122(12), 123608 (2019). [CrossRef]

11. Z.-R. Ma, X.-C. Huang, T.-J. Li, H.-C. Wang, G.-C. Liu, Z.-S. Wang, B. Li, W.-B. Li, and L.-F. Zhu, “First observation of new flat line fano profile via an x-ray planar cavity,” Phys. Rev. Lett. 129(21), 213602 (2022). [CrossRef]

12. R. Röhlsberger and J. Evers, “Quantum optical phenomena in nuclear resonant scattering,” in Modern Mössbauer Spectroscopy, G. L. Yutaka Yoshida, ed. (Springer Singapore, 2021), Topics in Applied Physics, pp. 105–171.

13. Y.-H. Chen, P.-H. Lin, G.-Y. Wang, A. Pálffy, and W.-T. Liao, “Transient nuclear inversion by x-ray free electron laser in a tapered x-ray waveguide,” Phys. Rev. Res. 4(3), L032007 (2022). [CrossRef]

14. 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 Radiat. 9(6), 383–393 (2002). [CrossRef]

15. T. Salditt, S. P. Krüger, C. Fuhse, and C. Bähtz, “High-transmission planar x-ray waveguides,” Phys. Rev. Lett. 100(18), 184801 (2008). [CrossRef]

16. M. Osterhoff and T. Salditt, “Coherence filtering of x-ray waveguides: analytical and numerical approach,” New J. Phys. 13(10), 103026 (2011). [CrossRef]

17. D. Pelliccia and D. M. Paganin, “Coherence filtering and revivals in x-ray waveguides: a communication-modes approach,” J. Opt. Soc. Am. A 31(8), 1754 (2014). [CrossRef]

18. T. Salditt and M. Osterhoff, X-ray Focusing and Optics in Nanoscale Photonic Imaging, T. Salditt, A. Egner, and D. R. Luke, eds. (Springer International Publishing, Cham, 2020), Topics in Applied Physics, pp. 71–124.

19. L. D. Caro, C. Giannini, D. Pelliccia, C. Mocuta, T. H. Metzger, A. Guagliardi, A. Cedola, I. Burkeeva, and S. Lagomarsino, “In-line holography and coherent diffractive imaging with x-ray waveguides,” Phys. Rev. B 77(8), 081408 (2008). [CrossRef]

20. 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]

21. M. Vassholz and T. Salditt, “Observation of electron-induced characteristic x-ray and bremsstrahlung radiation from a waveguide cavity,” Sci. Adv. 7(4), eabd5677 (2021). [CrossRef]

22. S.-Y. Lee, S. Ahrens, and W.-T. Liao, “Gravitationally sensitive structured x-ray optics using nuclear resonances,” arXivarXiv:2305.00613 (2023) .

23. P. Andrejic, L. M. Lohse, and A. Pálffy, “Waveguide QED with Mössbauer nuclei,” arXivarXiv:2305.11647 (2023).

24. D. Pelliccia, “Bending x-rays with nanochannels,” Nat. Photonics 10(2), 75–77 (2016). [CrossRef]

25. D. Marcuse, Theory of dielectric optical waveguides, Quantum electronics - principles and applications (Academic Press, 1991), 2nd ed.

26. C. Vassallo, “Radiating normal modes of lossy planar waveguides,” J. Opt. Soc. Am. 69(2), 311 (1979). [CrossRef]

27. J. J. Burke, “Propagation constants of resonant waves on homogeneous, isotropic slab waveguides,” Appl. Opt. 9(11), 2444 (1970). [CrossRef]

28. W. C. Chew, Waves and Fields in Inhomogenous Media (IEEE Press Series on Electromagnetic Wave Theory) (Wiley-IEEE Press, 1999).

29. M. J. N. van Stralen, K. F. I. Haak, and H. Blok, “On the classification of discrete modes in lossy planar waveguides: the modal analysis revisited,” Opt. Quantum Electron. 29(2), 243–262 (1997). [CrossRef]

30. W. H. Jark and S. D. Fonzo, “Transmission calculations for one- and two-dimensional x-ray waveguides,” in Design and Microfabrication of Novel X-Ray Optics II, vol. 5539 A. A. Snigirev and D. C. Mancini, eds. (SPIE, 2004).

31. D. Pelliccia, I. Bukreeva, A. Cedola, and S. Lagomarsino, “Dispersion properties of x-ray waveguides,” Appl. Opt. 45(12), 2821 (2006). [CrossRef]

32. S. Lagomarsino, I. Bukreeva, A. Cedola, D. Pelliccia, and W. Jark, “Theoretical analysis of x-ray waveguides,” in in Modern Developments in X-Ray and Neutron Optics, A. Erko, M. Idir, T. Krist, and A. G. Michette, eds. (Springer Berlin Heidelberg, 2008), Springer Series in optical science, (Springer Berlin Heidelberg, 2008), pp. 91–111.

33. I. Bukreeva, A. Cedola, D. Pellicia, W. Jark, and S. Lagomarsino, “X-ray waveguides,” in The X-Ray Standing Wave Technique, vol. 7 of Series on Synchrotron Radiation Techniques and Applications J. Zegenhagen and A. Kazimirov, eds. (World Scientific, 2013), pp. 143–162.

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

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

36. D. Pelliccia, I. Bukreeva, M. Ilie, W. Jark, A. Cedola, F. Scarinci, and S. Lagomarsino, “Computer simulations and experimental results on air-gap x-ray waveguides,” Spectrochim. Acta, Part B 62(6-7), 615–621 (2007). [CrossRef]

37. J. Soltau, L. M. Lohse, M. Osterhoff, and T. Salditt, “Finite-difference propagation for the simulation of x-ray multilayer optics,” Opt. Express 29(25), 41932 (2021). [CrossRef]

38. K. J. Tsanaktsidis, D. M. Paganin, and D. Pelliccia, “Analytical description of partially coherent propagation and absorption losses in x-ray planar waveguides,” Opt. Lett. 38(11), 1808–1810 (2013). [CrossRef]

39. F. Schenk, Optimization of resonances for multilayer x-ray resonators, vol. 3 of Göttingen series in x-ray physics (Universitätsverlag Göttingen, 2011).

40. I. Bukreeva, A. Popov, D. Pelliccia, A. Cedola, S. B. Dabagov, and S. Lagomarsino, “Wave-field formation in a hollow x-ray waveguide,” Phys. Rev. Lett. 97(18), 184801 (2006). [CrossRef]

41. C. Chen, P. Berini, D. Feng, S. Tanev, and V. Tzolov, “Efficient and accurate numerical analysis of multilayer planar optical waveguides in lossy anisotropic media,” Opt. Express 7(8), 260 (2000). [CrossRef]

42. I. Bukreeva, A. Sorrentino, A. Cedola, E. Giovine, A. Diaz, F. Scarinci, W. Jark, L. Ognev, and S. Lagomarsino, “Periodically structured x-ray waveguides,” J. Synchrotron Radiat. 20(5), 691–697 (2013). [CrossRef]

43. M. Osterhoff and T. Salditt, “Real structure effects in x-ray waveguide optics: The influence of interfacial roughness and refractive index profile on the near-field and far-field distribution,” Opt. Commun. 282(16), 3250–3256 (2009). [CrossRef]

44. F. Abelès, “Recherches sur la propagation des ondes électromagnétiques sinusoïdales dans les milieux stratifiés,” Annales de Physique 12(5), 596–640 (1950). [CrossRef]

45. M. Tolan, X-Ray Scattering from Soft-Matter Thin Films, Springer Tracts in Modern Physics (Springer Berlin Heidelberg, 1999).

46. D. Lentrodt, K. P. Heeg, C. H. Keitel, and J. Evers, “Ab initio quantum models for thin-film x-ray cavity QED,” Phys. Rev. Res. 2(2), 023396 (2020). [CrossRef]

47. P. Andrejić and A. Pálffy, “Superradiance and anomalous hyperfine splitting in inhomogeneous ensembles,” Phys. Rev. A 104(3), 033702 (2021). [CrossRef]

48. X. Kong, D. E. Chang, and A. Pálffy, “Green’s-function formalism for resonant interaction of x rays with nuclei in structured media,” Phys. Rev. A 102(3), 033710 (2020). [CrossRef]

49. S. Y. Buhmann, Dispersion Forces I, Springer Tracts in Modern Physics (Springer Berlin Heidelberg, 2013).

50. A. Asenjo-Garcia, J. D. Hood, D. E. Chang, and H. J. Kimble, “Atom-light interactions in quasi-one-dimensional nanostructures: A green’s-function perspective,” Phys. Rev. A 95(3), 033818 (2017). [CrossRef]

51. M. S. Tomaš, “Green function for multilayers: Light scattering in planar cavities,” Phys. Rev. A 51(3), 2545–2559 (1995). [CrossRef]

52. P. Johansson, “Electromagnetic green’s function for layered systems: Applications to nanohole interactions in thin metal films,” Phys. Rev. B 83(19), 195408 (2011). [CrossRef]

53. L. M. Lohse and P. Andrejic, “xwglib – a library for the simulation of x-ray waveguides,” https://gitlab.gwdg.de/irp/xwglib (2023).

54. G. W. Hanson and A. B. Yakovlev, Operator Theory for Electromagnetics (Springer New York, 2002).

55. A. Zettl, Sturm-Liouville Theory (American Mathematical Society, 2012).

56. A. Siegman, “Orthogonality properties of optical resonator eigenmodes,” Opt. Commun. 31(3), 369–373 (1979). [CrossRef]

57. P. T. Leung, S. Y. Liu, and K. Young, “Completeness and orthogonality of quasinormal modes in leaky optical cavities,” Phys. Rev. A 49(4), 3057–3067 (1994). [CrossRef]

58. I. Gohberg and S. Goldberg, “Spectral theory of integral operators,” in Basic Operator Theory, (Birkhäuser Boston, 2001), pp. 131–151.

59. M. Bôcher, “Certain cases in which the vanishing of the wronskian is a sufficient condition for linear dependence,” Trans. Amer. Math. Soc. 2(2), 139–149 (1901). [CrossRef]

60. R. E. Smith, G. W. Forbes, and S. N. Houde-Walter, “Unfolding the multivalued planar waveguide dispersion relation,” IEEE J. Quantum Electron. 29(4), 1031–1034 (1993). [CrossRef]

61. A. Bakhtazad, H. Abiri, and R. Ghayour, “A general transform for regularizing planar open waveguide dispersion relation,” J. Lightwave Technol. 15(2), 383–390 (1997). [CrossRef]

62. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photonics 1(1), 58–106 (2009). [CrossRef]

63. J. Chilwell and I. Hodgkinson, “Thin-films field-transfer matrix theory of planar multilayer waveguides and reflection from prism-loaded waveguides,” J. Opt. Soc. Am. A 1(7), 742 (1984). [CrossRef]

64. R. Parini, “cxroots: A Python module to find all the roots of a complex analytic function within a given contour,” (2018).

65. T. Schoonjans, A. Brunetti, B. Golosio, M. S. del Rio, V. A. Solé, C. Ferrero, and L. Vincze, “The xraylib library for x-ray–matter interactions. recent developments,” Spectrochim. Acta, Part B 66(11-12), 776–784 (2011). [CrossRef]

66. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University Press, 2012).

67. P. T. Kristensen and S. Hughes, “Modes and mode volumes of leaky optical cavities and plasmonic nanoresonators,” ACS Photonics 1(1), 2–10 (2014). [CrossRef]

68. O. Diekmann, D. Lentrodt, and J. Evers, “Inverse design approach to x-ray quantum optics with mössbauer nuclei in thin-film cavities,” Phys. Rev. A 105(1), 013715 (2022). [CrossRef]

polarization	$U$	$\partial _x U$	$\partial _z U$
TE	$E_y$	$i k \mu Z_0 H_z$	$-i k \mu Z_0 H_x$
TM	$H_y$	$-i k \epsilon Z_0^{-1} E_z$	$i k \epsilon Z_0^{-1} E_x$

Abstract

1. Introduction

2. Equations for the electromagnetic fields

2.1 Vector wave equations

2.2 Green’s functions

2.3 Transverse solutions in stratified media

2.4 Power flow

3. Free fields in stratified media

3.1 Separating the dimensions

3.2 Transversal eigenvalue problem

3.3 Resonant modes

3.4 Modified inner product

3.5 Wronskian resonance condition

3.6 Riemann sheets and branch cuts

3.7 Guided and leaky modes

3.8 Radiative modes

3.9 Mode dispersion

4. Line sources

4.1 Two-dimensional Green’s function

4.2 Homogeneous Green’s function

4.3 Semi-spectral representation

4.4 Derivation of real-space solutions

4.5 Computation of the residues

4.6 Dyadic Green’s functions for line sources

5. Point sources

5.1 Symmetry considerations

5.2 Derivation of the 3d real-space Green’s function

6. Transfer matrix formalism for layered media

6.1 Single layer

6.2 Multiple layers

6.3 Numerical solutions

7. Applications

7.1 Grazing incidence

7.2 Waveguide propagation

7.3 Dipole emission into waveguide modes

8. Conclusions

A. Integration by parts for the modified inner product

B. Derivation of mode dispersion

C. Derivation of the Wronskian derivative

D. Computing the scattered field

Funding

Acknowledgements

Disclosures

Data Availability

References

Data Availability

Cited By

Figures (8)

Tables (1)

Equations (164)

Optics Express