<figure><img src="image_1.png" /></figure>

# Resonant-state-expansion Born approximation for waveguides with dispersion

Mark Behzad Doost

► To cite this version:

Mark Behzad Doost. Resonant-state-expansion Born approximation for waveguides with dispersion. Physical Review A: Atomic, molecular, and optical physics [1990-2015], 2016, 93 (2), pp.023835. 10.1103/PhysRevA.93.023835. [https://arxiv.org/abs/1608.02385]

HAL Id: hal-01285305

https://hal.science/hal-01285305v2

Submitted on 14 Aug 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Resonant-state-expansion Born approximation for waveguides with dispersion

M. B. Doost

Independent Researcher

(Dated: August 14, 2016)

The resonant-state expansion (RSE) Born approximation, a rigorous perturbative method developed for electrodynamic and quantum mechanical open systems, is further developed to treat waveguides with a Sellmeier dispersion. For media that can be described by these types of dispersion over the relevant frequency range, such as optical glass, I show that the perturbed RSE problem can be solved by diagonalizing a second-order eigenvalue problem. In the case of a single resonance at zero frequency, this is simplified to a generalized eigenvalue problem. Results are presented using analytically solvable planar waveguides and parameters of borosilicate BK7 glass, for a perturbation in the waveguide width. The efficiency of using either an exact dispersion over all frequencies or an approximate dispersion over a narrow frequency range is compared. I included a derivation of the RSE Born approximation for waveguides to make use of the resonances calculated by the RSE, an RSE extension of the well-known Born approximation.

PACS numbers: 03.50.De, 42.25.-p, 03.65.Nk

# I. INTRODUCTION

Fundamental to scattering theory, the Born approximation consists of taking the incident field in place of the total field as the driving field at each point inside the scattering potential, it was first discovered by Born and presented in Ref. [1]. The Born approximation gave an expression for the differential scattering cross section in terms of the Fourier transform of the scattering potential. An important feature of this appearance of the Fourier transform is the availability of the inverse Fourier transform operation for the inverse scattering problem. The Born approximation is only valid for weak scatterers. In this paper I apply the RSE Born approximation Ref. [2, 3] to optical fibers or general open waveguide systems, which allows an arbitrary number of resonant-states (RSs) to be taken into account for scattering and transmission perturbative calculations.

Optical fibers provide a well-controlled optical path which can carry light over potentially very long distances of several hundred to thousands of kilometers since the light is contained within the fiber by total internal reflection. Fibers are clearly important for telecommunication since the light may be modulated to carry information. It has recently been reported that a hollow-core photonic-band gap fiber yielded a record combination of low loss and wide bandwidth [4]. Beyond telecommunications waveguides are used in integrated optical circuits [5], and terabit chip-to-chip interconnects [6].

Critical to understanding the response of a waveguide to being optically driven are its resonant states (RSs). The concept of RSs was first conceived and used by Gamow in 1928 in order to describe mathematically the process of radioactive decay, specifically the escape from the nuclear potential of an α particle by tunneling. Mathematically this corresponded to solving Schrödinger's equation for outgoing boundary conditions (BCs). These states have complex frequency ω with negative imaginary part meaning their time dependence $\exp(-i\omega t)$ decays

exponentially, thus giving an explanation for the exponential decay law of nuclear physics. The consequence of this exponential decay with time is that the further from the decaying system at a given instant of time, the greater the wave amplitude. An intuitive way of understanding this divergence of wave amplitude with distance is to notice that waves that are further away have left the system at an earlier time when less of the particle probability density had leaked out. There already exists numerical techniques for finding eigenmodes, such as the finite element method (FEM) and finite difference in time domain (FDTD) method to calculate resonances in open cavities. However, determining the effect of perturbations, which break the symmetry, presents a significant challenge as these popular computational techniques need large computational resources to model high-quality modes. Also these methods generate spurious solutions, which would damage the accuracy of the RSE Born approximation if included in the perturbation expansion basis.

In order to calculate the resonances of open systems we have recently developed the resonant-state expansion. Such an approach was not previously available due to the lack of a normalization for resonant states. I derived the normalization of resonant states as a contribution to Ref. [9], and the generalization of the normalisation to dispersive media I derived straight forwardly at the time of my rigorous derivation. The problem of normalising resonant states stems from the fact that RSs with complex frequencies have wave functions which are exponentially growing in space away from system, making a volume integral over all space as would be done for a hermitian system a meaningless exercise (see Appendix F). However the S-wave normalization was previously available and analytically correct but numerically unstable as I showed analytically in Ref [3].

So far the RSE has been applied to non dispersive systems of different dimensionality [7–11]. However, almost all realistic systems have a relevant frequency dispersion of the refractive index. I have recently found [3, 12] that
an Ohm's law dispersion, i.e. a term in the susceptibility scaling at the inverse frequency, can be introduced to the RSE while keeping its linearity. We found in Ref.[12] this dispersion can be a reasonable approximation for some materials over a limited range of wavelengths such as SHOTT BK7 glass over the optical range. In this work I generalize the RSE approach for waveguides detailed in Ref.[7] to systems constructed from dispersive media. Specifically I treat dispersion fully described by the Sellmeier equation, using a method similar to Muljarov et al in Ref.[13] for nano particles obeying Drude-Lorentz dispersion, however in this case generalised to inclined geometry. I also treat dispersion linear in wavelength squared, a generalisation of my work in Ref.[3, 12] to inclined geometry. I find the generalization of my Ohm's law approach to inclined geometry to be greatly superior to the generalization of the full dispersion treatment. This is most likely due to the unphysicality of the full dispersion at high and low frequencies, which are beyond the fitting range of the dispersion models. In order to make use of the RSs, I derive the RSE Born approximation for waveguides [2, 3], a theory for finding the field outside of a waveguide being internally or externally driven (see Appendix G).

The paper is organized as follows, Sec. II outlines the general recursive solution of Maxwell's equation using Green's function. Sec. III develops the solution in Sec. II into a perturbation theory for waveguides with linear dispersion in wavelength squared. Sec. IV develops the solution in Sec. II into a perturbation theory for waveguides with Sellmeier dispersion. Sec. VI demonstrates and compares the two approaches for including dispersion into the RSE. The perturbation considered corresponds to a narrowing of the waveguide and is discussed in Sec. VI C. In Sec. VI D I show the results using the simple dispersive approximation BK7. In Sec. VI E I use the full Sellmeier dispersion of BK7. Finally, we discuss the comparison and performance of the two methods in Sec. VI F, in Sec. VII, I derive the equations for the RSE Born approximation for planar waveguides. In Appendix C, D, E, and F I generalize my part of the proof of the normalization of RSs to open waveguide systems.

## II. RESONANT STATE EXPANSION FOR WAVEGUIDES AND NON-NORMAL INCIDENCE

In this section I develop the general recursive solution to the problem of calculating the resonant states of a perturbed waveguide. The recursive solution requires the Green's function (GF) of the unperturbed waveguide, and a perturbation which is retaining the translational invariance along the waveguide.

I first consider a waveguide of thickness $2a$ in vacuum with translational invariance in one direction, having the

dielectric constant

$$ \hat{\varepsilon}_{\omega}(\mathbf{r}) = \begin{cases} \varepsilon_{\omega} & \text{for } |\mathbf{r}| \le a, \\ 1 & \text{for } |\mathbf{r}| > a, \end{cases} \quad (1) $$

where $\varepsilon_{\omega}$ is the frequency dependent relative permittivity (RP) of the wave guide with the angular frequency $\omega$. $\dot{\mathbf{r}}$ are the coordinates normal to the waveguide, which for a cylindrical waveguide is given in polar coordinates $\mathbf{r} = (\rho, \theta)$ and for a planar waveguide $\dot{\mathbf{r}}$ is Cartesian $x$. I assume a relative permeability of $\mu = 1$ throughout this work. The electric field $\hat{\mathbf{E}}$ satisfies Maxwell's equation,

$$ \left[-\nabla \times \nabla \times -\hat{\varepsilon}_{\omega}(\mathbf{r}) \frac{1}{c^2} \frac{\partial^2}{\partial t^2}\right] \hat{\mathbf{E}}(\mathbf{r}, t) = 0, \quad (2) $$

For cylindrical and planar waveguides, due to translational invariance in the $z$ direction I first assume then prove that $\hat{\mathbf{E}}$ can be factorised as follows

$$ \hat{\mathbf{E}}(\mathbf{r}, t) = e^{i(pz - \omega t)} \mathbf{E}(\dot{\mathbf{r}}), \quad (3) $$

in which $p$ is the wave vector along the translationally invariant direction. For the component $\mathbf{E}(\dot{\mathbf{r}})$ of the electric field, Eq. (2) transforms to a one-dimensional (1D) wave equation in the planar case or a 2D wave-equation for the cylindrical waveguide to

$$ \left[-L + \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) \frac{\omega^2}{c^2}\right] \mathbf{E}(\dot{\mathbf{r}}) = 0. \quad (4) $$

Here $L$ is a linear operator not dependent on $z$ or $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$. The form of $L$ is derived in Appendix C, the derivation shows that the trial solution Eq. (3) is correct since it demonstrates that $\hat{\mathbf{E}}$ can be factorised by separation of variables. In principle $L$ could be any linear operator independent of both $z$ and $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$, and hence the treatment of waveguides in this paper is readily applicable to quantum mechanics or acoustics.

The non-normal incidence, characterized by $p \neq 0$, is treated here. The previously used spectral representation of the GF in the frequency domain contains a cut for $p \neq 0$, which can be removed by mapping the problem onto the complex normal wave-vector space $k$, as demonstrated in Ref.[7]. The relation between $k$, $p$, and $\omega$ is defined by us to be $\omega^2/c^2 = k^2 + p^2$, $k$ and $p$ are orthogonal by Pythagorean identity. I follow also here the approach used in Ref.[7] and formulate the RSE in the complex $k$-plane, for which the spectral representation of the GF of an infinite planar system with an in-plane momentum $p \neq 0$ written in the spectral form

$$ \hat{\mathbf{G}}_{k}^{(1)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_{n} \frac{\mathbf{E}_{n}(\dot{\mathbf{r}}) \otimes \mathbf{E}_{n}(\dot{\mathbf{r}}')}{2k_{n}(k - k_{n})}, \quad (5) $$

where $\mathbf{E}_n(\dot{\mathbf{r}})$ is the electric field of a RS, defined as an eigensolution of Eq. (4) with an arbitrary profile of $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ within a region $|\dot{\mathbf{r}}| < a$, satisfying the outgoing wave boundary conditions at infinity. The dyadic product of
the $\mathbf{E}_n(\dot{\mathbf{r}})$ fields is denoted by $\otimes$. The in-plane eigen-frequency $k_n$ are the poles in $\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$. By definition $\omega_n^2/c^2 = k_n^2 + p^2$ where $\omega_n^2$ is the eigenfrequency corresponding to the eigenstate $\mathbf{E}_n(\dot{\mathbf{r}})$. For a study of the derivation of such spectral representations see Appendix D.

The GF satisfies the equation

$$[-L + \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})(k^2 + p^2)] \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \hat{\mathbf{l}}\delta(\dot{\mathbf{r}} - \dot{\mathbf{r}}') \quad (6)$$

In the present work, we consider the refractive profile (RP) having the property

$$\lim_{\omega \to \infty} \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) = \hat{\varepsilon}(\dot{\mathbf{r}}), \quad (7)$$

where $\hat{\varepsilon}(\dot{\mathbf{r}})$ if the frequency independent part of $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$. If we assume $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ to be discontinuous through the boundary of the system then by substituting Eq. (5) into Eq. (6), convoluting with a finite function, and letting $k \to \infty$, we obtain the sum rule (see Appendix E)

$$\sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{k_n} = 0, \quad (8)$$

which allows us to re-write the Green's function a second way as

$$\hat{\mathbf{G}}_k^{(2)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n} \left[ \frac{1}{k-k_n} - \frac{k}{k^2+p^2} \right], \quad (9)$$

I now consider an arbitrary perturbation $\Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ of the dielectric constant inside the layer $|\dot{\mathbf{r}}| < a$. I use the Eq. (6) to solve the perturbed problem,

$$\mathcal{E}_{\nu}(\dot{\mathbf{r}}) = -\frac{\omega^2}{c^2} \int \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}') \mathcal{E}_{\nu}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}', \quad (10)$$

which is the solution to the equation

$$L \mathcal{E}_{\nu}(\dot{\mathbf{r}}) = \frac{\omega^2}{c^2} [\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) + \Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})] \mathcal{E}_{\nu}(\dot{\mathbf{r}}) \quad (11)$$

Note that the perturbed modes $\mathcal{E}_{\nu}(\dot{\mathbf{r}})$ satisfy Eq. (4) with $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ replaced by $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) + \Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ and the outgoing boundary conditions with $\omega_n$ replaced by $\omega$. I show the solution of this equation perturbatively for two different types of dispersion in the following two sections.

### III. SINGLE SELLMEIER RESONANCE AT ZERO FREQUENCY

For a Sellmeier dispersion with a single resonance at zero frequency (SRZ), I write the perturbation in Eq. (10) as

$$\Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) = \Delta\hat{\varepsilon}(\dot{\mathbf{r}}) + \frac{c^2 \Delta\hat{\sigma}(\dot{\mathbf{r}})}{\omega^2}, \quad (12)$$

again present only inside the layer $|\dot{\mathbf{r}}| < a$. It is then possible to linearise the RSE by using the different forms

of the Green's function, similar to Ref. [13], Eq. (5) and Eq. (9) for the different components of the perturbation in Eq. (10),

$$\begin{aligned} \mathcal{E}(\dot{\mathbf{r}}) = & -\int \hat{\mathbf{G}}_k^{(1)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta\hat{\sigma}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \\ & -\frac{\omega^2}{c^2} \int \hat{\mathbf{G}}_k^{(2)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta\hat{\varepsilon}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}', \end{aligned} \quad (13)$$

which results in the following relationship between unperturbed and perturbed modes

$$\begin{aligned} \mathcal{E}(\dot{\mathbf{r}}) = & - (k^2 + p^2) \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}})}{2k_n} \left[ \frac{1}{k-k_n} - \frac{k}{k^2+p^2} \right] \\ & \times \int_{-a}^{a} \mathbf{E}_n(\dot{\mathbf{r}}') \Delta\hat{\varepsilon}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}} \\ & - \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}})}{2k_n(k-k_n)} \\ & \times \int_{-a}^{a} \mathbf{E}_n(\dot{\mathbf{r}}') \Delta\hat{\sigma}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}. \end{aligned} \quad (14)$$

The perturbed mode $\mathcal{E}(\dot{\mathbf{r}})$ satisfies Eq. (4) with $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ replaced by $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) + \Delta\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ and outgoing boundary conditions with the wave vector $k$. In the interior region $|\dot{\mathbf{r}}| < a$ which contains the perturbation, the perturbed RSs $\mathcal{E}_{\nu}$ of wavenumbers $\kappa_{\nu}$ can be expanded into the unperturbed ones, exploiting the completeness of the latter which follows from Eq. (A6) (see also Appendix E),

$$\mathcal{E}_{\nu}(\dot{\mathbf{r}}) = \sum_{n} b_{n\nu} \mathbf{E}_{n}(\dot{\mathbf{r}}). \quad (15)$$

Substituting this expansion into Eq. (13) and equating coefficients at the same basis functions $\mathbf{E}_n$ results in the matrix equation

$$b_{n\nu} = \sum_{m} -b_{m\nu} \frac{(p^2 + \kappa_{\nu} k_n) V_{nm} + S_{nm}}{2k_n(\kappa_{\nu} - k_n)}, \quad (16)$$

where

$$V_{nm} = \int \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \Delta\hat{\varepsilon}(\dot{\mathbf{r}}) \mathbf{E}_m(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (17)$$

and

$$S_{nm} = \int \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \Delta\hat{\sigma}(\dot{\mathbf{r}}) \mathbf{E}_m(\dot{\mathbf{r}}) d\dot{\mathbf{r}}. \quad (18)$$

With the substitution $c_{n\nu} = b_{n\nu}\sqrt{k_n}$, Eq. (16) can be rewritten as

$$0 = \sum_m c_{m\nu} \left\{ \kappa_{\nu} \left( \frac{\delta_{nm}}{k_n} + \frac{V_{nm}}{2\sqrt{k_n k_m}} \right) + \frac{p^2 V_{nm} + S_{nm}}{2k_n \sqrt{k_n k_m}} - \delta_{nm} \right\} \quad (19)$$

care should be taken to take the sign of $\sqrt{k_n}$ consistently between matrix elements. Eq. (19) is linear in $\kappa_{\nu}$ and
can be solved by libraries for generalized linear matrix eigenvalue problems. In the absence of dispersion, $S_{nm} = 0$, and Eq. (19) reverts back to the expression for non-dispersive waveguides [7]. In the absence of $p, pa = 0$, we see that Eq. (19) reverts to an expression for dispersive perturbation to nano-particles [12].

## IV. SELLMEIER DISPERSION

In this section I develop the existing RSE for wave guides into a perturbation theory for waveguides with Sellmeier dispersion, following a similar approach to Ref.[13]. The Sellmeier dispersion is the sum of a set of Lorentzians in frequency with poles $q_j$ on the real axis. The starting point for this derivation is Eq. (10), the recursive solution for the perturbed problem derived in Sec. II.

In this section, we consider the RP

$$ \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) + \sum_{j} \frac{c^2 \hat{\sigma}_j(\dot{\mathbf{r}})}{\omega^2 - \Omega_j^2}, \quad (20) $$

and perturbation,

$$ \Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) = \Delta \hat{\varepsilon}(\dot{\mathbf{r}}) + \sum_{j} \frac{c^2 \Delta \hat{\sigma}_j(\dot{\mathbf{r}})}{\omega^2 - \Omega_j^2}, \quad (21) $$

with $j$ numbering the resonances at frequencies $\Omega_j$ having oscillator strengths $\hat{\sigma}_j$. We introduce an effective resonance wave vector $\hat{q}_j$ as $\hat{q}_j^2 = \Omega_j^2/c^2 - p^2$ and re-write Eq. (21) as,

$$ \begin{aligned} \Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) &= \Delta \hat{\varepsilon}(\dot{\mathbf{r}}) + \sum_{j} \frac{\Delta \hat{\sigma}_{j}(\dot{\mathbf{r}})}{k^2 - \hat{q}_{j}^2} \\ &= \Delta \hat{\varepsilon}(\dot{\mathbf{r}}) + \sum_{j} \left( \frac{\Delta \hat{\sigma}_{j}(\dot{\mathbf{r}})}{2k(k - \hat{q}_{j})} + \frac{\Delta \hat{\sigma}_{j}(\dot{\mathbf{r}})}{2k(k + \hat{q}_{j})} \right), \end{aligned} \quad (22) $$

In the Appendix A I use my method of [3], which was first adapted for full dispersion of 3D nano-particles by E. A. Muljarov in Ref.[13], to derive for non-vanishing $\hat{\sigma}_j$ the sum rule

$$ \sum_{n} \frac{\mathbf{E}_{n}(\dot{\mathbf{r}}) \otimes \mathbf{E}_{n}(\dot{\mathbf{r}}')}{2k_{n}(k_{n} \pm \hat{q}_{j})} = 0. \quad (23) $$

which allows us to write the Green's function as

$$ \hat{\mathbf{G}}_{k}^{(\pm j)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_{n} \frac{\mathbf{E}_{n}(\dot{\mathbf{r}}) \otimes \mathbf{E}_{n}(\dot{\mathbf{r}}')(k \pm \hat{q}_{j})}{2k_{n}(k - k_{n})(k_{n} \pm \hat{q}_{j})}, \quad (24) $$

which has the useful $(k \pm \hat{q}_j)$ in the numerator for reducing the order of the perturbation matrix problem through cancellation with perturbation denominators which would otherwise have lead to high order polynomial eigen problems.

We now consider a perturbation $\Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ of the RP inside the layer $|\dot{\mathbf{r}}| < a$. Similar to the previous section we use the Eq. (6) to solve the perturbed problem,

$$ \mathcal{E}(\dot{\mathbf{r}}) = -\frac{\omega^2}{c^2} \int \hat{\mathbf{G}}_k^{(2)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}', \quad (25) $$

where we take $\Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ of the same form as Eq. (22). Using both forms of the Green's function [13], Eq. (9) and Eq. (24) in Eq. (25), yields

$$ \begin{aligned} \mathcal{E}(\dot{\mathbf{r}}) = & -\frac{\omega^2}{c^2} \int \hat{\mathbf{G}}_k^{(2)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta \hat{\varepsilon}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \\ & -\frac{\omega^2}{c^2} \sum_{\pm, j} \int G_k^{(\pm j)}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \frac{\Delta \hat{\sigma}_j(\dot{\mathbf{r}}')}{2k(k \pm \hat{q}_j)} \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}', \end{aligned} \quad (26) $$

which results in the following relationship between unperturbed and perturbed modes

$$ \begin{aligned} \mathcal{E}(\dot{\mathbf{r}}) = & -(k^2 + p^2) \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}})}{2k_n} \left[ \frac{1}{k - k_n} - \frac{k}{k^2 + p^2} \right] \\ & \times \int \mathbf{E}_n(\dot{\mathbf{r}}') \Delta \hat{\varepsilon}(\dot{\mathbf{r}}') \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \\ & -(k^2 + p^2) \sum_{n, \pm, j} \frac{\mathbf{E}_n(\dot{\mathbf{r}})(k \pm \hat{q}_j)}{2k_n(k_n \pm \hat{q}_j)(k - k_n)} \\ & \times \int \mathbf{E}_n(\dot{\mathbf{r}}') \frac{\Delta \hat{\sigma}_j(\dot{\mathbf{r}}')}{2k(k \pm \hat{q}_j)} \mathcal{E}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}'. \end{aligned} \quad (27) $$

The perturbed mode $\mathcal{E}(\dot{\mathbf{r}})$ satisfies Eq. (4) with $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ replaced by $\hat{\varepsilon}_{\omega}(\dot{\mathbf{r}}) + \Delta \hat{\varepsilon}_{\omega}(\dot{\mathbf{r}})$ and outgoing boundary conditions with the wave vector $k$. In the interior region $|\dot{\mathbf{r}}| < a$ which contains the perturbation, the perturbed RSs $\mathcal{E}_{\nu}$ of wavenumbers $\kappa_{\nu}$ can be expanded into the unperturbed ones, exploiting the completeness of the latter which follows from Eq. (A6) (see also Appendix E),

$$ \mathcal{E}_{\nu}(\dot{\mathbf{r}}) = \sum_{n} b_{n\nu} \mathbf{E}_{n}(\dot{\mathbf{r}}). \quad (28) $$

Substituting this expansion into Eq. (27) and equating coefficients at the same basis functions $\mathbf{E}_n$ results in the matrix equation

$$ \begin{aligned} b_{n\nu} = & -\frac{p^2 + \kappa_{\nu} k_n}{2k_n(\kappa_{\nu} - k_n)} \sum_m V_{nm} b_{m\nu} \\ & -\frac{(\kappa_{\nu}^2 + p^2)}{2\kappa_{\nu}(\kappa_{\nu} - k_n)} \sum_{m,j} \frac{A_{nm}^{(j)}}{k_n^2 - \hat{q}_j^2} b_{m\nu}, \end{aligned} \quad (29) $$

where

$$ V_{nm} = \int \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \Delta \hat{\varepsilon}(\dot{\mathbf{r}}) \mathbf{E}_m(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (30) $$

and

$$ A_{nm}^{(j)} = \int \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \Delta \hat{\sigma}_j(\dot{\mathbf{r}}) \mathbf{E}_m(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (31) $$
is the matrix of the perturbation in the basis of unperturbed RSs. Introducing the abbreviation

$$U_{nm} = \sum_{j} \frac{A_{nm}^{(j)}}{2(k_{n}^{2} - \bar{q}_{j}^{2})} \quad (32)$$

I arrive at

$$\sum_{m} b_{m\nu} \left\{ p^2 U_{nm} + \kappa_{\nu} \left( \frac{p^2 V_{nm}}{2k_n} - \delta_{nm} k_m \right) + \kappa_{\nu}^2 \left( \delta_{nm} + \frac{V_{nm}}{2} + U_{nm} \right) \right\} = 0, \quad (33)$$

which is of second order in $\kappa_{\nu}$. We can write this matrix problem compactly as $Q(\kappa_{\nu})b_{\nu} = 0$ with $Q(\kappa_{\nu}) = \kappa_{\nu}^2 M + \kappa_{\nu} C + K$. To solve this matrix problem for a basis of size $N$, I follow Ref.[14] and use the first companion linearization, defining

$$R(\kappa_{\nu}) = \kappa_{\nu} \begin{bmatrix} M & 0 \\ 0 & I \end{bmatrix} + \begin{bmatrix} C & K \\ -I & 0 \end{bmatrix} \quad (34)$$

where $I$ is the $N$-by-$N$ identity matrix. With the corresponding vector

$$z = \begin{bmatrix} \kappa_{\nu} b_{\nu} \\ b_{\nu} \end{bmatrix} \quad (35)$$

We solve $R(\kappa_{\nu})z = 0$ using the generalized eigenvalue solver from the numerical algorithms group (NAG) C++ library. We then take the first $N$ components of $z$ as the eigenvector $b_{\nu}$.

## V. POSSIBLE EXPLANATION FOR THE DIFFERENCE IN ORDER BETWEEN THE TWO METHODS

Here I go beyond the algebraic mathematics presented so far and discuss the fundamental topological differences between the two RSE approaches that I have developed.

The RSE perturbation theory presented here is an objective continuous mapping function except in the neighbourhood of the finite number of poles such that,

$$\hat{\varepsilon}_{\omega} \omega^2 = [\hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})} + \sum_{j} [\hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})} \quad (36)$$

is mapped to

$$\hat{\varepsilon}_{\omega} \omega^2 + \Delta \hat{\varepsilon}_{\omega} \omega^2 = [\hat{\varepsilon}_{\omega} \omega^2 + \Delta \hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})} + \sum_{j} [\hat{\varepsilon}_{\omega}^{j} \omega^2 + \Delta \hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})}, \quad (37)$$

in the following way

$$\text{RSE}: [\hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})} \rightarrow [\hat{\varepsilon}_{\omega} \omega^2 + \Delta \hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})} \quad (38)$$

$$\text{RSE}: [\hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})} \rightarrow [\hat{\varepsilon}_{\omega}^{j} \omega^2 + \Delta \hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})} \quad (39)$$

$$\text{RSE}: \{k_n\} \rightarrow \{\kappa_{\nu}\} \quad (40)$$

$$\text{RSE}: \{E_n\} \rightarrow \{\mathcal{E}_{\nu}\} \quad (41)$$

where $[\hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})}$ and $[\Delta \hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})}$ have poles in the same position in frequency space while $[\hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})}$ and $[\hat{\varepsilon}_{\omega} \omega^2 + \Delta \hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})}$ are regular functions.

The mapping in Eq. (38) is between two spaces which are topologically equivalent to non-dispersive spaces. Therefore, the mapping problem is mathematically equivalent to non-dispersive RSE perturbation theory and so the order of the eigenvalue problem is not increased upon that of the non-dispersive problem.

However, in the case of mapping by Eq. (39) $[\hat{\varepsilon}_{\omega}^{j} \omega^2]^{(\text{poles})}$ contains poles in frequency space and so cannot be continuously deformed into $[\hat{\varepsilon}_{\omega} \omega^2]^{(\text{reg})}$ hence the map given by Eq. (39) is not equivalent to a non-dispersive mapping and so the order of the RSE eigenvalue problem is increased upon that of the non-dispersive problem.

Formally $\mathbf{A}$ and $\mathbf{B}$ are topologically equivalent if there is a homeomorphism mapping orbits of $\mathbf{A}$ to orbits of $\mathbf{B}$ homeomorphically, and preserving orientation of the orbits.

## VI. APPLICATION TO A PLANAR WAVEGUIDE

In this section I discuss the application of the Sellmeier waveguide RSE to an effectively 1D planar waveguide system, translationally invariant in the Cartesian $z$ and $y$ directions, described by a scalar RP, i.e., $\hat{\varepsilon}_{\omega}(x) = \hat{1}_{\varepsilon_{\omega}}(x)$, $\Delta \hat{\varepsilon}_{\omega}(x) = \hat{1}_{\Delta \varepsilon_{\omega}}(x)$. As unperturbed system I use a homogeneous planar wave guide of half width $a$, so that

$$\varepsilon_{\omega}(x) = \begin{cases} \varepsilon_{\omega} & \text{for } |x| < a, \\ 1 & \text{elsewhere.} \end{cases} \quad (42)$$

### A. Unperturbed resonant states

The solutions of Eq. (4), which satisfy the outgoing-wave boundary conditions in TE polarization take the form [7]

$$E_n(x) = \begin{cases} (-1)^n A_n e^{-ik_n x}, & x \le -a, \\ B_n[e^{iq_n x} + (-1)^n e^{-iq_n x}], & |x| \le a, \\ A_n e^{ik_n x}, & x \ge a, \end{cases} \quad (43)$$

where the eigenvalues $k_n$ satisfy the secular equation

$$(k_n - q_n) e^{iq_n a} + (-1)^n (k_n + q_n) e^{-iq_n a} = 0, \quad (44)$$

with

$$q_n = \sqrt{\epsilon_{\omega} k_n^2 + (\epsilon_{\omega} - 1) p^2}. \quad (45)$$
I use here an integer index n which takes even (odd) values for symmetric (anti-symmetric) RSs, respectively. The normalization constants A_n and B_n are found from the continuity of E_n across the boundaries and the normalization condition found in Appendix F.

In order to arrive at the normalization condition I consider an effectively 2D system and I take $\bar{E}(\dot{\mathbf{r}}, k)$ as being the analytic continuation into k space of $E_n(\dot{\mathbf{r}})$ and A being an arbitrary cross-section of the translational invariant waveguide to arrive at (see Appendix F).

$$1 + \delta_{k_n, 0} = \int_A \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \left. \frac{\partial(\omega^2 \hat{\varepsilon}_\omega(\dot{\mathbf{r}}))}{\partial(\omega^2)} \right|_{\omega=\omega_n} \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (46) \\ + \lim_{k \to k_n} \int_A \frac{\mathbf{E}_n(\dot{\mathbf{r}}) L \bar{\mathbf{E}}(\dot{\mathbf{r}}, k) - \bar{\mathbf{E}}(\dot{\mathbf{r}}, k) L \mathbf{E}_n(\dot{\mathbf{r}})}{k^2 - k_n^2} d\dot{\mathbf{r}}$$

outside the system in free space Maxwell's equation simplifies (see Appendixes F and G), therefore by extending A outside the system and denoting its circumference $L_A$ I can write,

$$1 + \delta_{k_n, 0} = \int_A \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \left. \frac{\partial(\omega^2 \hat{\varepsilon}_\omega(\dot{\mathbf{r}}))}{\partial(\omega^2)} \right|_{\omega=\omega_n} \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} (47) \\ + \lim_{k \to k_n} \int_{L_A} \frac{\mathbf{E}_n \cdot \nabla \bar{\mathbf{E}} - \bar{\mathbf{E}} \cdot \nabla \mathbf{E}_n}{k^2 - k_n^2} d\mathbf{L}.$$

I made the necessary assumption that $\hat{\varepsilon}_\omega$ is a (real) symmetric matrix or a scalar so that (defined to be in this case) $\mathbf{E} \cdot \hat{\varepsilon}_{\omega_n} \mathbf{E}_n = \mathbf{E}_n \cdot \hat{\varepsilon}_{\omega_n} \mathbf{E}$ and non dispersive at high frequencies (see Appendix E and F).

Various schemes exist to evaluate the line integral limit in Eq. (47) such as analytic methods in Ref. [9, 11] or numerically extending the surface into a non-reflecting, absorbing, perfectly matched layer where it vanishes. Hence, I derive from Eq. (47) the relevant normalization condition for the planar waveguide systems, which I will use for the numerical demonstration,

$$\int_{-a}^{a} \left. \frac{\partial(\omega^2 \varepsilon_{\omega}(x))}{\partial(\omega^2)} \right|_{\omega=\omega_n} E_n(x) E_n(x) dx \\ - \frac{E_n(-a)E_n(-a) + E_n(a)E_n(a)}{i2k_n} = 1 + \delta_{k_n, 0}. (48)$$

Here in Eq. (48) the first integral is taken over an arbitrary simply connected line enclosing the inhomogeneity of the system and the center of the coordinates used, and the second term is evaluated at its end points.

The coefficients in Eq. (43) take the form

$$B_n = \frac{(-i)^n}{\sqrt{4(a\epsilon_\omega + ip^2/(k_n\omega_n^2/c^2)) + M_n\eta}}, \quad (49)$$

$$A_n = B_n \frac{e^{iq_n a} + (-1)^n e^{-iq_n a}}{e^{ik_n a}}, \quad (50)$$

with

$$\eta = \frac{1}{\epsilon_{\omega}} \frac{\partial(\omega^2 \epsilon_{\omega})}{\partial(\omega^2)} - 1, \quad (51)$$

and

$$M_n = \epsilon_\omega \left[ \frac{(-1)^n 2 \sin(2q_n a)}{q_n} + 4a \right], \quad (52)$$

where $\omega_n^2/c^2 = k_n^2 + p^2$.

Around each pole of $\epsilon_\omega$ the secular equation Eq. (44) has a countable infinite number of solutions, each one creating a RS. An analytic approximation to these solutions is given in Appendix B, which are used as starting values for the numerical solution of Eq. (44) to determine these RSs closely spaced in the complex frequency plane.

## B. BK7 Sellmeier dispersion

In our numerical examples I use dispersion parameters describing the common borosilicate glass SCHOTT BK7. Its refractive index $n_r$ in the optical frequency range is well described by the Sellmeier expression

$$\epsilon_{\omega} = 1 + \frac{1.03961212\lambda^2}{\lambda^2 - 6000.69867 \text{ nm}^2} + \frac{0.231792344\lambda^2}{\lambda^2 - 20017.9144 \text{ nm}^2} + \frac{1.01046945\lambda^2}{\lambda^2 - 103.560653 \text{ µm}^2}. \quad (53)$$

$$\epsilon_{\omega} = 1 - \frac{6839.577136}{(a\omega/c)^2 - 6578.970180} - \frac{457.1302872}{(a\omega/c)^2 - 1972.154382} - \frac{0.3852016551}{(a\omega/c)^2 - 0.3812105899} \quad (54)$$

The resulting RP is real and is shown in Fig. 1. This dispersion is of the form Eq. (20) and can be treated using the quadratic matrix equation Eq. (33). In order to compare this result with the one of the linear matrix equation Eq. (19), I have fitted by $\epsilon_\omega$ over the wavelength range from 1.25 µm to 1.75 µm with the form of Eq. (12), yielding

$$\epsilon_{\omega} = 2.28239 - 0.01262 \lambda^2 / \mu\text{m}^2 \quad (55)$$

The deviations of this fit over the fitted range are around $10^{-4}$, as shown in Fig. 1.

## C. Unperturbed and perturbed system

The unperturbed system I consider in this work as example is a planar waveguide of width $2a = 2$ µm. I consider a size perturbation, narrowing the waveguide by 10 percent, as shown in Fig. 2. For the in-plane wave vector component I have chosen $pa = 5$. I choose as basis of size N for the RSE all RSs with

$$|k_n \sqrt{\epsilon_{\omega_n}}| < k_{\max}(N), \quad (56)$$

i.e., with a wave vector in the medium below a suitably chosen maximum $k_{\max}(N)$, this follows the approach of Ref. [13].
<figure><img src="image_3.png" /><figcaption>FIG. 1: RP $\epsilon_{\omega}$ of SCHOTT BK7 glass given by the Sellmeier expression Eq. (53) (solid line), and fitted using a single resonance at zero frequency (dashed line). The absolute residual of the fit is given as dotted line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 2: Cross-section of the perturbed and unperturbed wave guide used in Sec. VID and Sec. VIE, shown at a particular frequency.</figcaption></figure>

The motivation for the choice of basis selection criteria is the analogy between the RSE and a Fourier expansion. As I increase $k_n \sqrt{\epsilon_{\omega_n}}$, I increase the number of oscillations in the field as can be seen from Eq. (43). These increasingly oscillating fields go into the basis and similar to Fourier series, increase the resolution of the composite field generated by the expansion.

## D. Single Sellmeier resonance at zero frequency

Here I show results of the SRZ yielding the linear eigenvalue problem Eq. (19). For the unperturbed and perturbed system I use a RP given by the SRZ Eq. (55). The refractive index of the unperturbed and perturbed system is compared in Fig. 3 with the Sellmeier expression Eq. (53). I can see that for the chosen waveguide width, the fitted wavelength range corresponds to $3.6 < a\omega/c < 5$. The RS frequencies are given in Fig. 3(b) for the unperturbed system $\omega_n = c\sqrt{k_n^2 + p^2}$ and for the perturbed system $\omega_{\nu} = c\sqrt{\kappa_{\nu}^2 + p^2}$. To compare $\omega_{\nu}$

<figure><img src="image_5.png" /><figcaption>FIG. 3: SRZ RSE results for a thickness perturbation of a planar waveguide as function of the real part of $\omega$. (a) Refractive index of the unperturbed and perturbed medium (dashed line) described by the SRZ Eq. (55) and of the BK7 Sellmeier dispersion (solid line) Eq. (53). (b) RS frequencies for $N = 200$. Shown are exact unperturbed (open circles), exact perturbed (open squares), and RSE perturbed (crosses) data. (c) Relative error $\Delta$ of the RSE perturbed RS frequencies for $N = 50, 100, 200$. Closed symbols give relative error to the SRZ Eq. (55) RS eigen-frequencies as shown. Open symbols give relative error to the Sellmeier dispersion RS eigen-frequencies as shown. I show in (b) states with both incoming and outgoing boundary conditions, i.e. having positive or negative imaginary parts of complex resonant frequency, only states with outgoing boundary conditions go into the basis.</figcaption></figure>

calculated using the RSE with the exact result $\omega_{\nu}^{(\text{exact})}$ obtained by the secular Eq. (44), I define the relative error $\Delta = |\omega_{\nu}/\omega_{\nu}^{(\text{exact})} - 1|$, and give the resulting values in Fig. 3(c). I find that as I increase $N$, $\Delta$ decreases proportional to $N^{-3}$, similar to the findings for the non-dispersive RSE [8–10], and values in the $10^{-6}$ range are reached for $N = 200$. This is actually smaller than the relative error due to the SRZ approximation of the Sellmeier dispersion (see Fig. 1).

The evolution of the perturbed RSs wavenumbers $\omega_{\nu}$
<figure><img src="image_6.png" /></figure>

FIG. 4: SRZ RSE results for a thickness perturbation of a planar waveguide and $N = 50$. Shown are $\mathbb{R}(\omega_\nu)$ as function of the in-plane wave vector $p$ with $\mathfrak{S}(\omega_\nu)$ given by the color according to the scale shown.

<figure><img src="image_7.png" /></figure>

FIG. 5: SRZ RSE results for a thickness perturbation of a planar waveguide and $N = 50$. Shown are $\mathbb{R}(\omega_\nu)$ as function of the in-plane wave vector $p$ with the relative error $\Delta$ given by the color according to the scale shown.

with the in-plane wave vector $p$ is shown in Fig. 4. I can distinguish [7] the waveguide (WG) and anti-waveguide (AWG) modes, which have real $\bar{\omega}_\nu$, and the Fabry-Pérot (FP) modes which have a finite imaginary part representing their losses. Increasing $p$, the FP modes split into WG and AWG mode at the bifurcation point at which $pc = \bar{\omega}_\nu$, at which $k=0$ i.e. at grazing incidence of the external field.

The relative error $\Delta$ of $\bar{\omega}_\nu$ is given in Fig. 5 as function of $p$. I can see that $\Delta$ has generally a weak dependence on $p$, except close to the bifurcation points, where the error of the FP and AWG mode is significantly increased. Under closer examination I see further smaller discontinuities as function of $p$, which correspond to a change of the basis states included according to Eq. (56) for fixed $N = 50$. These results indicate that the RSE is able to reproduce all relevant RSs with a good accuracy.

### E. Sellmeier Dispersion

Here I show results of the Sellmeier RSE yielding the quadratic eigenvalue problem Eq. (33).

The RS frequencies are given in Fig. 6(b) for the unperturbed system $\omega_n = c\sqrt{k_n^2 + p^2}$ and for the perturbed system $\bar{\omega}_\nu = c\sqrt{\kappa_\nu^2 + p^2}$. The relative error $\Delta$ is given in Fig. 6(c). Also here I find that as we increase $N$, $\Delta$ decreases proportional to $N^{-3}$, and values in the $10^{-5}$ range are reached for $N = 800$.

The evolution of the perturbed RSs wavenumbers $\bar{\omega}_\nu$ with the in-plane wave vector $p$ is shown in Fig. 7 including the relative error $\Delta$, and in Fig. 4 including the imaginary part $\mathfrak{S}(\bar{\omega}_\nu)$.

I can see that $\Delta$ has generally a weak dependence on $p$,

except for the bifurcation points, at which a Fabry-Perot mode splits into a waveguide and anti-waveguide mode.

These results indicate that with full Sellmeier dispersion the RSE is able to reproduce all relevant RSs with a good accuracy.

### F. Performance Comparison

Here I compare the computational complexity of the SRZ RSE and the Sellmeier RSE. The corresponding eigenvalue problems Eq. (19) and Eq. (33) show that for the same number of basis states $N$, the SRZ RSE is a generalized eigenvalue problem of size $N \times N$, while the Sellmeier RSE has a size $2N \times 2N$ due to the quadratic nature of Eq. (33). Also I see from Fig. 8 that for $pa = 5$ the resonances closely associated to the Sellmeier $\hat{q}_j$ poles, those eigen-frequencies found using the solutions to Eq. (B4) as a starting point for the Newton-Raphson search, contribute approximately 10% to the basis size thereby increasing the numerical complexity and reducing the efficiency.

In Fig. 9 I show a comparison between the average relative error in the perturbed resonances for $pa = 5$ found in the range $3 < a\omega < 5$ versus the number of seconds required to diagonalise the perturbation eigenvalue problems using an Intel Core 2 Duo, 6M Cache, 3.16GHz, 1333MHz FSB processor connected to 4GB of RAM and using the NAG generalized eigenvalue problem solver software.

For the case of SRZ dispersion the relative error reaches the limit of about $5 \times 10^{-5}$ due to the RP approximation at $N \approx 70$. Once the perturbation method using the full Sellmeier dispersion reaches this relative error I see that upon adding further basis states more noise in the relative
<figure><img src="image_8.png" /><figcaption>FIG. 6: Sellmeier RSE results for a thickness perturbation of a planar waveguide. (a) refractive index of the Sellmeier dispersion given by Eq. (53). (b) RS frequencies for $N = 800$ basis states. Shown are exact unperturbed (open circles), exact perturbed (open squares), and RSE perturbed (crosses) data. (c) relative error $\Delta$ of the RSE perturbed RS frequencies $\omega_{\nu}$ for $N = 50, 100, 200, 400, 800$ basis states, as labeled. I show in (b) states with both incoming and outgoing boundary conditions, i.e. having positive or negative imaginary parts of complex resonant frequency, only states with outgoing boundary conditions go into the basis.</figcaption></figure>

error plot is generated and that the average relative error changes little. Fig. 9 suggests that the SRZ RSE is several orders of magnitude more efficient than the full Sellmeier RSE when considering the resonances over a small range of frequencies such as the small range used for optical communications.

## VII. RSE BORN APPROXIMATION FOR PLANAR WAVEGUIDES

In this section I derive the RSE Born approximation for planar waveguide, the equivalent theory for effectively 2D systems is given in Appendix G.

In effectively one dimension Eq. (2) becomes in free

<figure><img src="image_9.png" /><figcaption>FIG. 7: Sellmeier RSE results for a thickness perturbation of a planar waveguide and $N = 200$. Shown are $\Re(\omega_{\nu})$ as function of the in-plane wave vector $p$ with the relative error $\Delta$ given by the colour according to the scale shown.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIG. 8: Graph showing quantitatively how the number of poles in the basis associated closely with the poles in the Sellmeier equation for $pa = 5$, those eigenfrequencies found using the solutions to Eq. (B4) as a starting point for the Newton-Raphson search, increase with basis size $N$ when the basis is selected according to Eq. (56). For reference I also show the number of ordinary poles, those found with the Newton-Raphson method without requiring the solutions of Eq. (B4) as a starting point for their search and discovery, in the basis versus basis size $N$.</figcaption></figure>

space (taking $c = 1$),

$$ \left[ \nabla^2 - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right] \theta_p E(x) = 0, \qquad (57) $$

However, $\nabla^2 \theta_p = -p^2 \theta_p$, therefore

$$ \left[ \frac{d^2}{dx^2} + k^2 \right] E(x) = 0, \qquad (58) $$
<figure><img src="image_11.png" /><figcaption>FIG. 9: Comparison of the average relative error of the wave guide modes occurring in the range $3 < a\omega < 5$ calculated using the two different perturbation schemes in this paper versus time required to diagonalise the perturbation matrix using an Intel Core 2 Duo, 6M Cache, 3.16GHz, 1333MHz FSB processor connected to 4GB of RAM and using the NAG generalised eigenvalue problem solver software</figcaption></figure>

Hence the free space GF equation is,

$$ \left[ \frac{d^2}{dx^2} + k^2 \right] G_k^{fs}(x, x') = \delta(x - x'), \quad (59) $$

which has the solution,

$$ G_k^{fs}(x, x') = -\frac{e^{ik|x-x'|}}{2ik} \quad (60) $$

The systems associated with $G_k^{fs}$ and $G_k$ of Eq. (6) are related by the Dyson Equations perturbing back and forth with $\Delta\varepsilon_\omega(x) = \varepsilon_\omega(x) - 1$ similar to Ref. [2],

$$ G_k(x, x'') = G_k^{fs}(x, x'') \quad (61) \\ -\omega^2 \int G_k^{fs}(x, x''') \Delta\varepsilon_\omega(x'') G_k(x'', x'') dx'', $$

$$ G_k(x'', x'') = G_k^{fs}(x'', x'') \quad (62) \\ -\omega^2 \int G_k(x'', x') \Delta\varepsilon_\omega(x') G_k^{fs}(x', x'') dx', $$

Combining Eq. (61) and Eq. (62) it is obtained similar to Ref. [2]

$$ G_k(x, x'') = G_k^{fs}(x, x'') \\ - \omega^2 \int G_k^{fs}(x, x') \Delta\varepsilon_\omega(x') G_k^{fs}(x', x'') dx' \\ + \omega^4 \iint G_k^{fs}(x, x') \Delta\varepsilon_\omega(x') G_k(x', x'') dx' \\ \times \Delta\varepsilon_\omega(x''') G_k^{fs}(x'', x'') dx'' dx'. \quad (63) $$

Hence, in one dimension the RSE Born approximation can be greatly simplified by using Eq. (60) and the spectral GF Eq. (E7) in Eq. (63) to arrive at,

$$ G_k(x, x'') = -\frac{e^{ik(x''-x)}}{2ik} + \frac{\omega^2 e^{ik(x''-x)}}{4k^2} \int_{-a}^{a} \Delta\varepsilon_\omega dx' \\ - \frac{\omega^4 e^{ik(x''-x)}}{4k^2} \sum_n \frac{A_n(x) A_n(x'')}{2k(k-k_n)}. \quad (64) $$

where $A_n$ is defined as the Fourier transform,

$$ A_n(x) = \int_{-a}^{a} e^{ik\hat{x}x'} \Delta\varepsilon_{\omega}(x') E_n(x') dx' \quad (65) $$

Interestingly in one dimension I do not require the far field approximation to make the simplification of the Green's function required to bring the RSE Born approximation to the form of Eq. (G9). Hence, in one dimension the RSE Born Approximation is valid everywhere outside of the slab and not just in the far field. I note that fast Fourier transforms are available for use upon Eq. (65).

I demonstrate the computational accuracy of the RSE Born approximation in Fig. 10 where I calculate the transmission defined as

$$ T(k, x') = |2kG(x', -a; k)|^2 \quad (66) $$

For comparison the analytic GF is found by solving Maxwell's wave equation in one dimension with a source of plane waves while making use of Maxwell's boundary conditions. The system treated is the unperturbed system in Sec. VI D with $pa = 5$. From Fig. 10 we can see that unlike the standard Born approximation the RSE Born approximation is valid over an arbitrarily wide range of $k$ depending only on the basis size $N$ used. Furthermore we see that as the basis size increases the RSE Born approximation converges to the exact solution. The absolute error in the RSE Born approximation is approximately reduced by an order of magnitude each time the basis size is doubled. Absolute errors of $10^{-8} - 10^{-5}$ are seen in the $k$ range shown for basis size $N = 401$.

VIII. SUMMARY

In this work I have extended the RSE to media having a simple dispersion linear in wavelength squared. This dispersion has a single pole at zero frequency and does not introduce an additional dynamic degree of freedom as it would be the case for a more general Sellmeier model of the material response. This property allows to keep the simplicity of the RSE formulation, therefore retaining the advantage of the RSE in computational efficiency discussed in Ref. 9.

Furthermore, in this work I have also extended the RSE to waveguides obeying the Sellmeier equation for glasses. In order to do this I have reduced the order
<figure><img src="image_12.png" /><figcaption>FIG. 10: I give transmission results for the unperturbed system in Sec. VID with $pa = 5$. I have used analytic modes as RSs for the RSE Born approximation. (a) Exact transmission as a function of frequency as defined by Eq. (66) (b) Absolute error in transmission calculated using the analytic form of $T(k, x')$ between $x = -a$ and $x = a$ versus frequency $\omega$ as comparison for the RSE Born approximation. Here $N = 51, 101, 201, 401$ as labeled.</figcaption></figure>

of the eigenvalue problem to second order by using sum rules as was done in [13].

To make use of the RSs generated by the RSE I derived the RSE Born approximation for effectively 1D and 2D waveguides.

My conclusion is that the simple dispersion treatment is more efficient than using the full Sellmeier dispersion. This is because the dispersion only has to be correctly reproduced over a narrow part of the optical frequency range and therefore it is inefficient to include the full and at some frequencies unphysical Sellmeier dispersion.

## Acknowledgments

I acknowledge support by the Cardiff University EP-SRC Doctoral Prize Fellowship EP/M50631X/1. I thank E. A. Muljarov for his positive and highly valuable contribution to this paper, which was to present me with his copy of his quite complete draft of Ref. [13]. I thank W. Langbein and T. Wood for their positive contributions to this paper during the early stages.

## Appendix A: Sum rule

Here I use my approach of Ref. [3], which was first adapted for full dispersion of 3D nano-particles by E. A. Muljarov in Ref. [13], to derive sum rules and completeness of modes for the waveguide GF.

By substituting Eq. (5) into Eq. (6) and using Eq. (4) I have,

$$ -\hat{\mathbf{1}}\delta(\hat{\mathbf{r}} - \hat{\mathbf{r}}') = \sum_{n} \frac{(\omega_n^2 \hat{\varepsilon}_{\omega_n} - \omega^2 \hat{\varepsilon}_{\omega})}{c^2 2k_n(k - k_n)} \mathbf{E}_n(\hat{\mathbf{r}}) \otimes \mathbf{E}_n(\hat{\mathbf{r}}'). \quad (A1) $$

I now consider the pole in $\hat{\varepsilon}_{\omega}$ at $k = \pm \hat{q}_j$

$$ \hat{\varepsilon}_{\omega} = \frac{\hat{\sigma}_{j}(\hat{\mathbf{r}})}{2k(k \pm \hat{q}_{j})} \quad (A2) $$

where I have been able to ignore other terms in the dispersion since they are constant in the limit $k \to \mp \hat{q}_j$. Using $\omega^2/c^2 = k^2 + p^2$ and $\omega_n^2/c^2 = k_n^2 + p^2$ and Eq. (A2) in Eq. (A1), I find

$$ \begin{aligned} & = \frac{(\omega_n^2 \hat{\varepsilon}_{\omega_n} - \omega^2 \hat{\varepsilon}_{\omega})}{c^2 2k_n(k - k_n)} \\ & = 2\hat{\sigma}_j(\hat{\mathbf{r}}) \frac{kp^2}{8k_n^2 k(k - k_n)(k_n \pm \hat{q}_j)} \\ & - \frac{2\hat{\sigma}_j(\hat{\mathbf{r}})}{(k \pm \hat{q}_j)} \left( \frac{p^2 k_n}{8k_n^2 k(k - k_n)} \pm \frac{k_n \hat{q}_j k}{8k_n^2 k(k_n \pm \hat{q}_j)} \right) \end{aligned} \quad (A3) $$

Then convoluting Eq. (A1) with arbitrary finite function and taking the limit $k \to \mp \hat{q}_j$, I see that the second term is diverging unless

$$ \sum_n \frac{\mathbf{E}_n(\hat{\mathbf{r}}) \otimes \mathbf{E}_n(\hat{\mathbf{r}}')}{k_n(k_n \pm \hat{q}_j)} = 0. \quad (A4) $$

since when $k \to \mp \hat{q}_j$

$$ \begin{aligned} & \left( \frac{p^2 k_n}{8k_n^2 k(k - k_n)} \pm \frac{k_n \hat{q}_j k}{8k_n^2 k(k_n \pm \hat{q}_j)} \right) \\ & = \pm \frac{p^2 + \hat{q}_j^2}{8k_n(k_n \pm \hat{q}_j)\hat{q}_j} \end{aligned} \quad (A5) $$

The closure and over completeness follow from letting $\omega \to \infty$ to give [3]

$$ \frac{\hat{\varepsilon}(\hat{\mathbf{r}})}{2} \sum_{n} \mathbf{E}_{n}(\hat{\mathbf{r}}) \otimes \mathbf{E}_{n}(\hat{\mathbf{r}}') = \hat{\mathbf{1}}\delta(\hat{\mathbf{r}} - \hat{\mathbf{r}}'), \quad (A6) $$

## Appendix B: Analytic solution of secular equation close to RP poles

Here I calculate an analytic approximation for the solutions of the secular equation Eq. (44) close to the poles of the RP, following the approach of Ref. [13]. In the limit
$\epsilon_{\omega} \to \infty$ I find from Eq. (45) that $q_n \to \sqrt{\epsilon_{\omega}}\omega_n/c$, so that Eq. (44) is given by $-e^{iq_n a} + (-1)^n e^{-iq_n a} = 0$, and thus

$$ \begin{aligned} \sin(\sqrt{\epsilon_{\omega}}\omega_n a/c) &= 0 & \text{for even } n \\ \cos(\sqrt{\epsilon_{\omega}}\omega_n a/c) &= 0 & \text{for odd } n \end{aligned} \quad (B1) $$

In the limit $k_n \to \pm \hat{q}_j$ to the poles of the RP given in Eq. (22) I find

$$ \epsilon_{\omega} \frac{\omega_n^2}{c^2} = \sigma_j \frac{k_n^2 + p^2}{2k_n(k_n \mp \hat{q}_j)} = \gamma_n \quad (B2) $$

and solutions of Eq. (B1) are given by

$$ \gamma_n a^2 = (\pi n/2)^2 \quad (B3) $$

The solutions of Eq. (B2) with $k_n$ close to $\pm \hat{q}_j$ are given by the analytic solution of the quadratic formula

$$ (\sigma_j - 2\gamma_n)k_n^2 \pm 2\gamma_n\hat{q}_j k_n + \sigma_j p^2 = 0 \quad (B4) $$

I use the solutions of Eq. (B4) as starting points for the Newton-Raphson search for RSs.

The RSs of the unperturbed system (see Sec. VI C) with Sellmeier dispersion Eq. (53) are given in Fig. 11. I note that close to the lowest resonance frequency of the dispersion at $\Omega_j a/c \approx 0.6174$, there are a large number of RSs approaching the pole from smaller frequencies. These RSs arise because the refractive index is diverging to positive infinity on the low frequency side, allowing for a countable infinite number of WG and AWG modes to form. $\sigma_j$ is negative for all the resonances in Eq. (53).

## Appendix C: Calculation of L

Here I derive the form of the $L$ operator occurring in the reduced Maxwell's wave equation for waveguide Eq. (4).

For planar systems $L$ is given in Ref. 7 as

$$ -L = \frac{d^2}{dx^2} - p^2 \quad (C1) $$

For effectively two-dimensional systems with one direction of translational invariance $L$ can be calculated as follows. Due to translational invariance in the direction $z$, the direction of propagation $I$ can write the full field as,

$$ \hat{\mathbf{E}}(\mathbf{r}) = e^{ipz - i\omega t} \mathbf{E}(\dot{\mathbf{r}}) = \theta_p \mathbf{E}(\dot{\mathbf{r}}), \quad (C2) $$

where $\dot{\mathbf{r}}$ give the coordinates in the plane perpendicular to $z$.

Therefore using standard vector identities I can write

$$ \begin{aligned} \nabla \times \nabla \times \hat{\mathbf{E}} &= \nabla \times \nabla \times (\theta_p \mathbf{E}) \\ &= \nabla \times (\nabla \theta_p \times \mathbf{E} + \theta_p \nabla \times \mathbf{E}) \\ &= \nabla \theta_p (\nabla \cdot \mathbf{E}) - \mathbf{E} (\nabla \cdot (\nabla \theta_p)) \\ &\quad + (\mathbf{E} \cdot \nabla) \nabla \theta_p - (\nabla \theta_p \cdot \nabla) \mathbf{E} \\ &\quad + \nabla \theta_p \times \nabla \times \mathbf{E} \\ &\quad + \theta_p \nabla \times \nabla \times \mathbf{E} \end{aligned} \quad (C3) $$

<figure><img src="image_13.png" /><figcaption>FIG. 11: (a) RP as a function of frequency given by Eq. (53). (b) Resonant frequencies $\omega_n$ of the homogeneous dielectric waveguide with $pa = 5$ and this RP, forming the basis states for the Sellmeier RSE discussed in Sec. VI E. (c) Zoom of (b) showing the series of poles on the low frequency side of the singularity in RP. I show in (b) and (c) states with both incoming and outgoing boundary conditions, i.e. having positive or negative imaginary parts of complex resonant frequency, only states with outgoing boundary conditions go into the basis.</figcaption></figure>

I make use of the following identities,

$$ \nabla \theta_p = ip \theta_p e_z \quad (C4) $$

$$ \nabla \cdot \nabla \theta_p = -p^2 \theta_p \quad (C5) $$

$$ \nabla \theta_p (\nabla \cdot \mathbf{E}) = ip \theta_p (\nabla \cdot \mathbf{E}) e_z \quad (C6) $$

$$ (\mathbf{E} \cdot \nabla) \nabla \theta_p = (e_z \cdot \mathbf{E}) p^2 \theta_p e_z \quad (C7) $$

$$ \partial_z \mathbf{E} = 0 \quad (C8) $$

where $e_z$ is the unit vector in the $z$ direction so as to simplify Eq. (C3) to

$$ \begin{aligned} \nabla \times \nabla \times \hat{\mathbf{E}} &= \nabla \times \nabla \times (\mathbf{E} \theta_p) \\ &= ip \theta_p e_z (\nabla \cdot \mathbf{E}) + \mathbf{E} p^2 \theta_p \\ &\quad + (e_z \cdot \mathbf{E}) p^2 \theta_p e_z - ip \theta_p \partial_z \mathbf{E} \\ &\quad + ip \theta_p e_z \times \nabla \times \mathbf{E} + \theta_p \nabla \times \nabla \times \mathbf{E} \\ &= \theta_p L(\mathbf{E}) \end{aligned} \quad (C9) $$
Hence I see that,

$$L(\mathbf{E}) = ipe_z(\nabla \cdot \mathbf{E}) + \mathbf{E}p^2 + (e_z \cdot \mathbf{E})p^2 e_z + ipe_z \times \nabla \times \mathbf{E} + \nabla \times \nabla \times \mathbf{E} \quad (\text{C10})$$

So $L$ is a linear operator independent of the RP and $z$ as required. Therefore the eigenvalue problems derived in Sec. III and Sec. IV are valid for planar and cylindrical wave guides. Also the derivation of $L$ demonstrates I can use separation of variables in Eq. (2) to arrive at Eq. (C2).

## Appendix D: Spectral representation of the GF of an open system

Here I almost exactly repeat my derivations which I contributed to Ref. [8] using exactly the same method but with increased mathematical rigour, in order to prove in this section the spectral representation of the Green's function (GF) of a general wave equations.

The GF of an open waveguide equation system is a tensor function $\hat{\mathbf{G}}_k$ which satisfies the outgoing wave BCs and the waveguide wave equation Eq. (4) with a delta function source term,

$$-L\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') + \hat{\varepsilon}(\dot{\mathbf{r}}, k)k^2\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \hat{\mathbf{1}}\delta(\dot{\mathbf{r}} - \dot{\mathbf{r}}'), \quad (\text{D1})$$

In this Appendix the effective permittivity $\hat{\varepsilon}(\dot{\mathbf{r}}, k) = \hat{\varepsilon}_\omega(\dot{\mathbf{r}})(1 + p^2/k^2)$ with fixed $p$. Assuming a simple-pole structure of the GF inside the scatterer with poles at $k = q_n$ and taking into account its large-$k$ vanishing asymptotics, the Mittag-Leffler theorem allows us to express the GF only inside the scatterer as the convergent

$$\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}')}{k - q_n}. \quad (\text{D2})$$

It will probably be that Eq. (D2) will need to be experimentally verified by comparing scattering predicted by the RSE Born approximation with experimental results of scattering from well defined scatterers. It may be that Eq. (D2) is a fundamental law of Physics. My justification for this form of the GF is the superposition of Lorentzians which make up the scattering profile of resonators and the numerical verification of this form of GF made in Ref. [7–10].

Assuming no degeneracy with the mode $n$, the definition of the residue tensor $\hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$ at a simple pole of the function $\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$ is,

$$\lim_{k \to q_n} (k - q_n) \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \quad (\text{D3})$$

I have again assumed $\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$ to be holomorphic in this neighbourhood of $q_n$ except for at the poles $q_n$ so that it

has a Laurent series at $q_n$. Substituting the expression Eq. (D2) into Eq. (D3) gives

$$\lim_{k \to q_n} (k - q_n) \sum_m \frac{\hat{\mathbf{Q}}_m(\dot{\mathbf{r}}, \dot{\mathbf{r}}')}{k - q_m} = \hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \quad (\text{D4})$$

so that

$$\lim_{k \to q_n} (k - q_n) \sum_{m \neq n} \frac{\hat{\mathbf{Q}}_m(\dot{\mathbf{r}}, \dot{\mathbf{r}}')}{k - q_m} = 0 \quad (\text{D5})$$

Substituting the expression Eq. (D2) into Eq. (D1) and convoluting with an arbitrary finite vector function $\mathbf{D}(\dot{\mathbf{r}})$ over a finite volume $V$ we obtain

$$\sum_n \frac{-L\mathbf{F}_n(\dot{\mathbf{r}}) + \hat{\varepsilon}(\dot{\mathbf{r}}, k)k^2\mathbf{F}_n(\dot{\mathbf{r}})}{k - q_n} = \mathbf{D}(\dot{\mathbf{r}}), \quad (\text{D6})$$

where $\mathbf{F}_n(\dot{\mathbf{r}}) = \int_V \hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}')\mathbf{D}(\dot{\mathbf{r}}')d\dot{\mathbf{r}}'$. Multiplying by $(k-q_n)$ and taking the limit $k \to q_n$ yields

$$\lim_{k \to q_n} (k - q_n) \sum_m \frac{-L\mathbf{F}_m(\dot{\mathbf{r}}) + \hat{\varepsilon}(\dot{\mathbf{r}}, k)k^2\mathbf{F}_m(\dot{\mathbf{r}})}{k - q_m} \quad (\text{D7})$$

$$= \lim_{k \to q_n} (k - q_n) \mathbf{D}(\dot{\mathbf{r}}) = 0. \quad (\text{D8})$$

From Eq. (D5) we can see,

$$\lim_{k \to q_n} (k - q_n) \sum_{m \neq n} \frac{-L\mathbf{F}_m(\dot{\mathbf{r}}) + \hat{\varepsilon}(\dot{\mathbf{r}}, k)k^2\mathbf{F}_m(\dot{\mathbf{r}})}{k - q_m} = 0, \quad (\text{D9})$$

so we can drop terms $n \neq m$ from the summation in Eq. (D8) to give

$$\lim_{k \to q_n} (k - q_n) \frac{-L\mathbf{F}_n(\dot{\mathbf{r}}) + \hat{\varepsilon}(\dot{\mathbf{r}}, k)k^2\mathbf{F}_n(\dot{\mathbf{r}})}{k - q_n} = 0. \quad (\text{D10})$$

or

$$-L\mathbf{F}_n(\dot{\mathbf{r}}) + \hat{\varepsilon}(\dot{\mathbf{r}}, q_n)q_n^2\mathbf{F}_n(\dot{\mathbf{r}}) = 0. \quad (\text{D11})$$

Due to the convolution with the GF, $\mathbf{F}_n(\dot{\mathbf{r}})$ satisfies the same outgoing wave BCs. Then, according to Eq. (4), $\mathbf{F}_n(\dot{\mathbf{r}}) \propto \mathbf{E}_n(\dot{\mathbf{r}})$ and $q_n = q_n$, i.e.

$$L\mathbf{E}_n(\dot{\mathbf{r}}) = \hat{\varepsilon}(\dot{\mathbf{r}}, k_n)k_n^2\mathbf{E}_n(\dot{\mathbf{r}}), \quad (\text{D12})$$

Note that the convolution of the kernel $\hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$ with different vector functions $\mathbf{D}(\dot{\mathbf{r}})$ can be proportional to one and the same vector function $\mathbf{E}_n(\dot{\mathbf{r}})$ only if the kernel has the form of a product:

$$\hat{\mathbf{Q}}_n(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}') / 2k_n, \quad (\text{D13})$$

where $\otimes$ is the dyadic product operator.

The symmetry in Eq. (D13) follows from the reciprocity theorem, described mathematically by the relation

$$s_1 \hat{\mathbf{G}}_k(\dot{\mathbf{r}}_1, \dot{\mathbf{r}}_2) s_2 = s_2 \hat{\mathbf{G}}_k(\dot{\mathbf{r}}_2, \dot{\mathbf{r}}_1) s_1, \quad (\text{D14})$$
which holds for any two point sources $s_{1,2}$ at points $\dot{\mathbf{r}}_{1,2}$ emitting at the same frequency. Hence $\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}')$ is symmetric.

In the case of a GF made up of degenerate modes the proof of Eq. (D13) is modified by making use of orthogonality of the degenerate modes to choose $\mathbf{D}(\dot{\mathbf{r}})$ such that,

$$ \int_V \mathbf{E}_m(\dot{\mathbf{r}}) \cdot \mathbf{D}(\dot{\mathbf{r}}) d\dot{\mathbf{r}} = 0, \qquad (\text{D15}) $$

for $m \neq n$ and where state $m$ is degenerate with $n$. Hence I obtain

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n(k-k_n)}. \qquad (\text{D16}) $$

## Appendix E: Derivation of sum rule and completeness

Here I exactly repeat my derivation of the sum rule and completeness of the GF which I made for Ref. [3] but in more detail.

In order to simplify the RSE Born approximation we require an appropriate spectral form of the GF which is different from the one already proven in Appendix D. To obtain this correct form I start with the GF valid inside the scatterer only,

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n(k-k_n)}. \qquad (\text{E1}) $$

Substituting Eq. (E1) in

$$ -L\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') + (k^2+p^2)\hat{\varepsilon}_\omega(\dot{\mathbf{r}})\hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \hat{\mathbf{1}}\delta(\dot{\mathbf{r}}-\dot{\mathbf{r}}'), \qquad (\text{E2}) $$

gives for $k \to \infty$,

$$ \hat{\varepsilon}(\dot{\mathbf{r}}) \sum_n \frac{(k+k_n)\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n} = \hat{\mathbf{1}}\delta(\dot{\mathbf{r}}-\dot{\mathbf{r}}'). \qquad (\text{E3}) $$

since throughout the derivation in this appendix we are considering the limit where $k \to \infty$ at which $\hat{\varepsilon}_k(\dot{\mathbf{r}}) = \hat{\varepsilon}(\dot{\mathbf{r}})$, i.e. the system is non-dispersive at high frequencies.

Convoluting Eq. (E3) with arbitrary finite functions $\mathbf{D}$ which are only non-zero inside the resonator, gives

$$ \int \left( \hat{\varepsilon}(\dot{\mathbf{r}}) \sum_n \frac{(k+k_n)\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n} \right) \mathbf{D}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' = \mathbf{D}(\dot{\mathbf{r}}). \qquad (\text{E4}) $$

and assuming the series are convergent as $k \to \infty$ we get,

$$ \int \left( \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n} \right) \mathbf{D}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' = 0. \qquad (\text{E5}) $$

and since this is true for all $\mathbf{D}$ it follows that

$$ \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k_n} = 0. \qquad (\text{E6}) $$

Combining Eq. (E1) and Eq. (E6) yields

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k(k-k_n)}. \qquad (\text{E7}) $$

Combining Eq. (E3) and Eq. (E6) leads to the closure relation

$$ \frac{\hat{\varepsilon}(\dot{\mathbf{r}})}{2} \sum_n \mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}') = \hat{\mathbf{1}}\delta(\dot{\mathbf{r}} - \dot{\mathbf{r}}'), \qquad (\text{E8}) $$

which expresses the completeness of the RSs, so that any function can be written as a superposition of RSs. If in the perturbed system some of the series are not convergent or are instead conditionally convergent then we will not arrive at the sum rule and completeness, in which case I expect that the RSE Born approximation will still give convergence to the exact solution but only if a valid spectral Green's function is used, such as Eq. (E1).

## Appendix F: Derivation of the correct normalisation for effectively 2D waveguides

Here I almost exactly repeat my derivations which I contributed to Ref. [9] using exactly the same method as for Maxwell's equations in order to prove in this section that the spectral representation, only valid inside the waveguide (similar to Ref. [3] for nano-particles), rigorously derived in Appendix C, D, and E

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{E}_n(\dot{\mathbf{r}}')}{2k(k-k_n)}, \qquad (\text{F1}) $$

leads to the RS normalization condition Eq. (F8) for general waveguide equations. To do so, I consider an analytic continuation $\mathbf{E}(\dot{\mathbf{r}}, k)$ of the vector wave function $\mathbf{E}_n(\dot{\mathbf{r}})$ around the point $k = k_n$ in the complex $k$-plane ($k_n$ is the wavenumber of the given RS). I select the analytic continuation such that it satisfies the outgoing wave boundary condition and my waveguide equation (taking $c = 1$)

$$ -L\mathbf{E}(\dot{\mathbf{r}}, k) + \hat{\varepsilon}(\dot{\mathbf{r}}, \omega)(k^2+p^2)\mathbf{E}(\dot{\mathbf{r}}, k) = (k^2-k_n^2)\boldsymbol{\sigma}(\dot{\mathbf{r}}) \qquad (\text{F2}) $$

with an arbitrary source term.

The source $\boldsymbol{\sigma}(\dot{\mathbf{r}})$ has to be zero outside the cross-section of the inhomogeneity of $\hat{\varepsilon}_\omega(\dot{\mathbf{r}})$ for the field $\mathbf{E}(\dot{\mathbf{r}}, k)$ to satisfy the outgoing wave boundary condition. It also has to be non-zero somewhere inside that cross-section, as otherwise $\mathbf{E}(\dot{\mathbf{r}}, k)$ would be identical to $\mathbf{E}_n(\dot{\mathbf{r}})$. It is further required that $\boldsymbol{\sigma}(\dot{\mathbf{r}})$ is normalized according to

$$ \int_A \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \boldsymbol{\sigma}(\dot{\mathbf{r}}) d\dot{\mathbf{r}} = 1 + \delta_{k_n, 0}, \qquad (\text{F3}) $$

The integral in Eq. (F3) is taken over an arbitrary cross-section $A$ which includes all system inhomogeneity of $\hat{\varepsilon}_\omega(\dot{\mathbf{r}})$. I do not derive Eq. (F3), it is simply a convenient condition to put on the otherwise arbitrary spatial
dependence of the source which lies only inside the waveguide scatterer. If I had made the condition anything else (which in my earliest proof I did) then algebra and cancellation would lead us back to the same result, however with more mathematical complexity and operations. The condition make the mathematics easier because it causes $\mathbf{E} \rightarrow \mathbf{E}_n$ exactly and without proportionality constant appearing, equation (F3) ensures that the analytic continuation reproduces $\mathbf{E}_n(\dot{\mathbf{r}})$ in the limit $k \rightarrow k_n$. Solving Eq. (F2) with the help of the GF and using its spectral representation Eq. (F1), I find:

$$ 
\begin{aligned}
\mathbf{E}(\dot{\mathbf{r}}, k) &= \int_A \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}')(k^2 - k_n^2)\boldsymbol{\sigma}(\dot{\mathbf{r}}')d\dot{\mathbf{r}}' \\
&= \sum_m \mathbf{E}_m(\dot{\mathbf{r}}) \frac{k^2 - k_n^2}{2k(k - k_m)} \int_A \mathbf{E}_m(\dot{\mathbf{r}}') \cdot \boldsymbol{\sigma}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}', \quad (\text{F4})
\end{aligned}
$$

and then, using Eq. (F3), obtain

$$ \lim_{k \to k_n} \mathbf{E}(\dot{\mathbf{r}}, k) = \mathbf{E}_n(\dot{\mathbf{r}}), $$

for any $\dot{\mathbf{r}}$ inside the system. Outside the system, the analytic continuation $\mathbf{E}(\dot{\mathbf{r}}, k)$ is defined as a solution of the waveguide equation wave equation in free space. This solution is connected to the field inside the system [given by Eq. (F4)] through the boundary conditions. Note that in the case of degenerate modes, $k_m = k_n$ for $m \neq n$, the current $\boldsymbol{\sigma}(\dot{\mathbf{r}})$ has to be chosen in such a way that it satisfies Eq. (F3) and, additionally,

$$ \int_A \mathbf{E}_m(\dot{\mathbf{r}}) \cdot \boldsymbol{\sigma}(\dot{\mathbf{r}}) d\dot{\mathbf{r}} = 0, $$

in order that the degenerate modes can be normalised separately.

I now consider the integral

$$ I_n(k) = \frac{\int_A (\mathbf{E} \cdot L\mathbf{E}_n - \mathbf{E}_n \cdot L\mathbf{E}) d\dot{\mathbf{r}}}{k^2 - k_n^2} \quad (\text{F5}) $$

and evaluate it by using the waveguide equations Eqs. (4) and (F2) for $\mathbf{E}_n$ and $\mathbf{E}$, respectively, and the source term normalization Eq. (F3):

$$ I_n(k) = \frac{\int_A (\mathbf{E} \cdot \omega_n^2 \hat{\mathbf{e}}_{\omega_n} \mathbf{E}_n - \mathbf{E}_n \cdot \omega^2 \hat{\mathbf{e}}_{\omega} \mathbf{E}) d\dot{\mathbf{r}}}{k^2 - k_n^2} + 1 + \delta_{k_n, 0}, \quad (\text{F6}) $$

where I assume that $\hat{\mathbf{e}}_{\omega}$ is a a (real) symmetric matrix or a scalar so that (defined to be in this case) $\mathbf{E} \cdot \hat{\mathbf{e}}_{\omega_n} \mathbf{E}_n = \mathbf{E}_n \cdot \hat{\mathbf{e}}_{\omega_n} \mathbf{E}$, and so I obtain also by commutation of $\mathbf{E}$ and $\mathbf{E}_n$ and simple calculus in the integral shown Eq. (F6) the dispersion factor and normalization becomes,

$$ 
\begin{aligned}
1 + \delta_{k_n, 0} = & \int_A \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \frac{\partial(\omega^2 \hat{\mathbf{e}}_{\omega}(\dot{\mathbf{r}}))}{\partial(\omega^2)} \Big|_{\omega=\omega_n} \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (\text{F7}) \\
& + \lim_{k \to k_n} \int_A \frac{\mathbf{E}_n(\dot{\mathbf{r}}) L \mathbf{E}(\dot{\mathbf{r}}, \mathbf{k}) - \mathbf{E}(\dot{\mathbf{r}}, \mathbf{k}) L \mathbf{E}_n(\dot{\mathbf{r}})}{k^2 - k_n^2} d\dot{\mathbf{r}}
\end{aligned}
$$

and then finally with $L_A$ the circumference of A,

$$ 
\begin{aligned}
1 + \delta_{k_n, 0} = & \int_A \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \left. \frac{\partial \omega^2 \hat{\mathbf{e}}(\dot{\mathbf{r}}, \omega)}{\partial \omega^2} \right|_{\omega=\omega_n} \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} \quad (\text{F8}) \\
& + \lim_{k \to k_n} \oint_{L_A} \frac{\mathbf{E}_n \cdot \nabla \mathbf{E} - \mathbf{E} \cdot \nabla \mathbf{E}_n}{k^2 - k_n^2} d\mathbf{L}.
\end{aligned}
$$

I have also assumed a simpler form for the waveguide equation in free space, see Appendix G so that I can use the vector and scalar Green's theorem to simplify Eq. (F7). Various schemes exist to evaluate the line integral limit in Eq. (F8) such as analytic methods in Ref. [11] or numerically extending the closed arc into a non-reflecting, absorbing, perfectly matched layer where it vanishes.

Due to the over completeness of the basis demonstrated in Appendix E we do not have the relation we find for Hermitian systems

$$ \int_{allspace} \frac{\omega_n^2 \hat{\mathbf{e}}_{\omega_n}(\dot{\mathbf{r}}) - \omega_m^2 \hat{\mathbf{e}}_{\omega_m}(\dot{\mathbf{r}})}{\omega_m^2 - \omega_n^2} \mathbf{E}_m(\dot{\mathbf{r}}) \cdot \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} = 0, \quad (\text{F9}) $$

for all $m \neq n$. The basis states of non-Hermitian systems are not all mutually orthogonal although some are orthogonal to one another. Hence we cannot derive for non-Hermitian systems the normalization relation, which is well known for Hermitian systems and follows from Eq. (F9),

$$ \int_{allspace} \left. \frac{\partial(\omega^2 \hat{\mathbf{e}}_{\omega}(\dot{\mathbf{r}}))}{\partial(\omega^2)} \right|_{\omega=\omega_n} \mathbf{E}_n(\dot{\mathbf{r}}) \cdot \mathbf{E}_n(\dot{\mathbf{r}}) d\dot{\mathbf{r}} = 1 + \delta_{k_n, 0}, \quad (\text{F10}) $$

and, hence, Eq. (F8) and Eq. (F10) are different.

## Appendix G: RSE Born Approximation for effectively 2D waveguides

In effectively two dimensions Eq. (2) becomes in free space (taking $c = 1$),

$$ \left[ \nabla^2 - \frac{1}{c^2} \frac{\partial^2}{\partial t^2} \right] \theta_p \mathbf{E}(\dot{\mathbf{r}}) = 0, \quad (\text{G1}) $$

However $\nabla^2 \theta_p = -p^2 \theta_p$, therefore

$$ [\nabla^2 + k^2] \mathbf{E} = 0, \quad (\text{G2}) $$

Hence the free space GF equation is,

$$ [\nabla^2 + k^2] \hat{\mathbf{G}}_k^{fs}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = \hat{\mathbf{1}} \delta(\dot{\mathbf{r}} - \dot{\mathbf{r}}'), \quad (\text{G3}) $$

which has the solution, also assuming $\rho \gg \rho'$,

$$ \hat{\mathbf{G}}_k^{fs}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}'|) \hat{\mathbf{1}} \approx Q \sqrt{\frac{1}{\rho\pi}} e^{ik|\dot{\mathbf{r}} - \dot{\mathbf{r}}'|} \hat{\mathbf{1}} \quad (\text{G4}) $$

The systems associated with $G_k^{fs}$ and $G_k$ of Eq. (6) are related by the Dyson Equations perturbing back and forth with $\Delta \hat{\mathbf{e}}_{\omega}(\dot{\mathbf{r}}) = \hat{\mathbf{e}}_{\omega}(\dot{\mathbf{r}}) - \hat{\mathbf{1}}$ similar to Ref. [2],
$$ \hat{\mathbf{G}}_k(\mathbf{r}, \dot{\mathbf{r}}'') = \hat{\mathbf{G}}_k^{\mathrm{fs}}(\mathbf{r}, \dot{\mathbf{r}}'') \qquad (\mathrm{G5}) \\ -\omega^2 \int \hat{\mathbf{G}}_k^{\mathrm{fs}}(\mathbf{r}, \dot{\mathbf{r}}'') \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}'') \hat{\mathbf{G}}_k(\dot{\mathbf{r}}''', \dot{\mathbf{r}}'') d\dot{\mathbf{r}}'', $$

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}''', \dot{\mathbf{r}}'') = \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}''', \dot{\mathbf{r}}'') \qquad (\mathrm{G6}) \\ -\omega^2 \int \hat{\mathbf{G}}_k(\dot{\mathbf{r}}''', \dot{\mathbf{r}}') \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}', \dot{\mathbf{r}}'') d\dot{\mathbf{r}}', $$

Combining Eq. (G5) and Eq. (G6) it is obtained similar to Ref. [2]

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') \\ - \omega^2 \int \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}', \dot{\mathbf{r}}'') d\dot{\mathbf{r}}' \\ + \omega^4 \int \int \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}, \dot{\mathbf{r}}') \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') \hat{\mathbf{G}}_k(\dot{\mathbf{r}}', \dot{\mathbf{r}}''') \\ \times \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}'') \hat{\mathbf{G}}_k^{\mathrm{fs}}(\dot{\mathbf{r}}''', \dot{\mathbf{r}}'') d\dot{\mathbf{r}}''' d\dot{\mathbf{r}}'. \quad (\mathrm{G7}) $$

In order to improve the numerical performance further I make a final few steps as in the original Born approximation [1], I define unit vector $\hat{\mathbf{r}}$ such that $\hat{\mathbf{r}} = \rho\hat{\mathbf{r}}$ and $\mathbf{k}_s = k\hat{\mathbf{r}}$. Then for $\rho \gg \rho'$,

$$ k|\dot{\mathbf{r}} - \dot{\mathbf{r}}'| = k\rho - \mathbf{k}_s \cdot \dot{\mathbf{r}}' + O\left(\frac{1}{\rho^2}\right) + \dots \quad (\mathrm{G8}) $$

Therefore substituting Eq. (5) and Eq. (G4) in to Eq. (G7) and using Eq. (G8) because both $\dot{\mathbf{r}}, \dot{\mathbf{r}}''$ are far from the scatterer I arrive at the RSE Born approximation

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q^2 e^{ik(\rho+\rho'')}}{\pi \sqrt{\rho \rho''}} \int e^{i(\mathbf{k}_s - \mathbf{k}_s'') \cdot \dot{\mathbf{r}}'} \Delta \hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \\ + \omega^4 \frac{Q^2 e^{ik(\rho+\rho'')}}{\pi \sqrt{\rho \rho''}} \sum_n \frac{\mathbf{A}_n(\mathbf{k}_s) \otimes \mathbf{A}_n(-\mathbf{k}_s'')}{2k(k-k_n)} \quad (\mathrm{G9}) $$

or using Eq. (E1) instead

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q^2 e^{ik(\rho+\rho'')}}{\pi \sqrt{\rho \rho''}} \int e^{i(\mathbf{k}_s - \mathbf{k}_s'') \cdot \dot{\mathbf{r}}'} \Delta \hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \\ + \omega^4 \frac{Q^2 e^{ik(\rho+\rho'')}}{\pi \sqrt{\rho \rho''}} \sum_n \frac{\mathbf{A}_n(\mathbf{k}_s) \otimes \mathbf{A}_n(-\mathbf{k}_s'')}{2k_n(k-k_n)} \quad (\mathrm{G10}) $$

The vector $\mathbf{A}_n$ is defined as a Fourier transform of the RSs,

$$ \mathbf{A}_n(\mathbf{k}_s) = \int e^{i\mathbf{k}_s \cdot \dot{\mathbf{r}}'} \Delta\hat{\mathbf{\epsilon}}_{\omega}(\dot{\mathbf{r}}') \mathbf{E}_n(\dot{\mathbf{r}}') d\dot{\mathbf{r}}' \quad (\mathrm{G11}) $$

I note that the fast Fourier transform method is available. The first two terms in Eq. (G9) correspond to the standard Born approximation, the final summation term corresponds to the RSE correction to the Born approximation.

A simple corollary of this theory is as follows, I can see from the arguments just stated that from Eq. (G5) if $\dot{\mathbf{r}}''$ is inside the resonator and $\rho \gg \rho''$ then

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q e^{ik\rho}}{\sqrt{\rho\pi}} \sum_n \frac{\mathbf{A}_n(\mathbf{k}_s) \otimes \mathbf{E}_n(\dot{\mathbf{r}}'')}{2k(k-k_n)}, \quad (\mathrm{G12}) $$

or using Eq. (E1) instead

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q e^{ik\rho}}{\sqrt{\rho\pi}} \sum_n \frac{\mathbf{A}_n(\mathbf{k}_s) \otimes \mathbf{E}_n(\dot{\mathbf{r}}'')}{2k_n(k-k_n)}, \quad (\mathrm{G13}) $$

similarly from Eq. (G6) if $\dot{\mathbf{r}}$ is inside the resonator and $\rho'' \gg \rho$ then

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q e^{ik\rho''}}{\sqrt{\rho''\pi}} \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{A}_n(-\mathbf{k}_s'')}{2k(k-k_n)}, \quad (\mathrm{G14}) $$

or using Eq. (E1) instead

$$ \hat{\mathbf{G}}_k(\dot{\mathbf{r}}, \dot{\mathbf{r}}'') = -\frac{i}{4} H_0(k|\dot{\mathbf{r}} - \dot{\mathbf{r}}''|) \hat{\mathbf{1}} \\ - \omega^2 \frac{Q e^{ik\rho''}}{\sqrt{\rho''\pi}} \sum_n \frac{\mathbf{E}_n(\dot{\mathbf{r}}) \otimes \mathbf{A}_n(-\mathbf{k}_s'')}{2k_n(k-k_n)}, \quad (\mathrm{G15}) $$

other permutations are possible.

It may be possible to use these formulas to develope a hollow glass cable with light and fluid passing along it, the escaping light being measured in the far field for analysis with the RSE Born approximation for the determination of the content of the fluid via inverse emission problem. To elaborate on this point, it was shown in Ref. [3] that

$$ \lim_{\rho \to \infty} \mathbf{E}_n(\dot{\mathbf{r}}) \propto \mathbf{A}_n(\dot{\mathbf{r}}k_n) \frac{e^{ik_n \rho}}{\rho^{1/2}} \quad (\mathrm{G16}) $$

However close to a sharp resonance, when $[k \approx \Re(k_n), \rho$ is fixed], scattering is dominated by a single resonance and so the scattered $\mathbf{E}$-field $\mathbf{E}^{\mathrm{scattered}}$ is approximately

$$ \lim_{\rho \to \infty} \mathbf{E}^{\mathrm{scattered}}(\dot{\mathbf{r}}, k) \approx C \mathbf{A}_n(\dot{\mathbf{r}}k) \frac{e^{ik\rho}}{\rho^{1/2}} \quad (\mathrm{G17}) $$

Hence $\mathbf{A}_n(\dot{\mathbf{r}}k)$ can be partial inverse Fourier transformed with respect to angle to find information about the internal structure of the cable, $C$ is some constant. The
same approach can be taken in three dimensions to aid the determination of resonator structure.

The reason Eq. (G16) takes the form it does can be seen more clearly by substituting the GF of Eq. (G12) into the first line of Eq. (F4) and taking the limit $k \to k_n$, while following the arguments in that section that,

$$ \lim_{k \to k_n} \mathbf{E}(\dot{\mathbf{r}}, k) = \mathbf{E}_n(\dot{\mathbf{r}}). $$

The 3D equivalent [2, 3] of Eq. (G16) one day might be valuable information for solving the inverse emission problem from such things as black hole gravitational wave emitters, i.e. calculations of the potential

from the emission (decay) via fast inverse Fourier transform methods upon the set of $A_n(\hat{\mathbf{r}}k_n)$, particularly if the potentials of interest are rotating about a fixed axis so we know their orientation to some extent such as occurs for decaying magnetic nuclei as part of a non-magnetic crystaline compound placed inside a NMR (nuclear-magnetic-resonance) machine. Because $k_n$ are discrete values $A_n(\hat{\mathbf{r}}k_n)$ only give angular information when inverse Fourier transformed and so the inverse Fourier methods might have to be used self-consistently in conjunction with the RSE perturbation theory and the values of $k_n$. This is a highly speculative aside and might be a possible topic for future research.

[1] Max Born, Zeitschrift für Physik 38 802 (1926)

[2] R. -C. Ge, P. T. Kristensen, Jeff. Young, S. Huges, New J. Phys. 16 113048 (2014).

[3] M. B. Doost arXiv:1508.04103.

[4] F. Poletti et al., Nature Photonics 7, 279 (2013).

[5] S. John, Nature (London) 460, 337 (2009).

[6] J. C. Cervantes-Gonzalez, D. Ahn, X. Zheng, S. K. Banerjee, A. T. Jacome, J. C. Campbell, and I. E. Zaldivar-Huerta, Appl. Phys. Lett. 101, 261109 (2012).

[7] Armitage, L. J. and Doost, M. B. and Langbein, W. and Muljarov, E. A., Phys. Rev. A, 89, (2014).

[8] Doost, M. B. and Langbein, W. and Muljarov, E. A., Phys. Rev. A, 87, 043827, (2013).

[9] Doost, M. B. and Langbein, W. and Muljarov, E. A., Phys. Rev. A, 90, 013834, (2014).

[10] Doost, M. B. and Langbein, W. and Muljarov, E. A., Phys. Rev. A, 85, 023835, (2012).

[11] E. A. Muljarov and W. Langbein and R. Zimmermann, Europhys. Lett., 5, 50010, (2010).

[12] arXiv:1508.03851 Doost, M. B. and Langbein, W. and Muljarov, E. A.

[13] E. A. Muljarov, W. Langbein arXiv:1510.01182

[14] F. Tisseur and K. Meerbergen, The quadratic eigenvalue problem, SIAM Rev., 43 (2001), pp. 235-286