%~Mouliné par MaN_auto v.0.32.2 2024-05-16 15:22:02 \documentclass[CRMATH,Unicode,XML,biblatex]{cedram} \usepackage{amssymb} \DeclareMathOperator{\curl}{curl} \DeclareMathOperator{\Curl}{Curl} \addbibresource{crmath20230376.bib} \newcommand*\bfx{\mathbf{x}} \newcommand*\bfE{\mathbf{E}} \newcommand*\bfH{\mathbf{H}} \newcommand*\bfu{\mathbf{u}} \newcommand*\bfU{\mathbf{U}} \newcommand*\bfT{\mathbf{T}} \newcommand*\bfX{\mathbf{X}} \newcommand*\bfe{\mathbf{e}} \newcommand*\bfnu{\boldsymbol{\nu}} \newcommand*\bftv{\boldsymbol{\tau}} \newcommand*\hboldE{\hat{\boldE}} \newcommand{\rext}{\mathrm{ext}} \newcommand{\rint}{\mathrm{int}} \newcommand{\SC}{\calK} \newcommand{\indSC}{p} \renewcommand{\Im}{\operatorname{Im}} \newcommand{\tv}{\boldtau} \newcommand{\ol}[1]{\overline{#1}} \newcommand{\bpm}{\begin{pmatrix}} \newcommand{\epm}{\end{pmatrix}} \renewcommand{\div}{\operatorname{div}} \newcommand{\boldtau}{\boldsymbol{\tau}} \newcommand{\setC}{\mathbb{C}} \newcommand{\setN}{\mathbb{N}} \newcommand{\setR}{\mathbb{R}} \newcommand{\boldu}{\mathbf{u}} \newcommand{\boldw}{\mathbf{w}} \newcommand{\boldx}{\mathbf{x}} \newcommand{\boldE}{\mathbf{E}} \newcommand{\boldH}{\mathbf{H}} \newcommand{\boldL}{\mathbf{L}} \newcommand{\calE}{\mathcal{E}} \newcommand{\calH}{\mathcal{H}} \newcommand{\calJ}{\mathcal{J}} \newcommand{\calK}{\mathcal{K}} \newcommand{\calU}{\mathcal{U}} \newcommand{\calX}{\mathcal{X}} \newcommand{\calY}{\mathcal{Y}} \graphicspath{{./figures/}} \newcommand*{\mk}{\mkern -1mu} \newcommand*{\Mk}{\mkern -2mu} \newcommand*{\mK}{\mkern 1mu} \newcommand*{\MK}{\mkern 2mu} \newcommand*{\relabel}{\renewcommand{\labelenumi}{(\theenumi)}} \newcommand*{\romanenumi}{\renewcommand*{\theenumi}{\roman{enumi}}\relabel} \newcommand*{\Romanenumi}{\renewcommand*{\theenumi}{\Roman{enumi}}\relabel} \newcommand*{\alphenumi}{\renewcommand*{\theenumi}{\alph{enumi}}\relabel} \newcommand*{\Alphenumi}{\renewcommand*{\theenumi}{\Alph{enumi}}\relabel} \let\oldtilde\tilde \renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}} \let\oldhat\hat \renewcommand*{\hat}[1]{\mathchoice{\widehat{#1}}{\widehat{#1}}{\oldhat{#1}}{\oldhat{#1}}} \TopicEN{Partial differential equations} \TopicFR{\'Equations aux dérivées partielles} \title{On the analysis of modes in a closed electromagnetic waveguide} \alttitle{Sur l'analyse des modes dans un guide d'ondes électromagnétiques fermé} \author{\firstname{Martin} \lastname{Halla}\CDRorcid{0000-0002-3010-3540}} \address{Institut f\"ur Numerische und Angewandte Mathematik, Georg-Augst Universit\"at G\"ottingen, Lotzestr.\ 16-18, 37083 G\"ottingen, Deutschland.} \email[M. Halla]{m.halla@math.uni-goettingen.de} \thanks{The first author was supported by DFG project 468728622 and DFG SFB 1456 project 432680300.} \author{\firstname{Peter} \lastname{Monk}\CDRorcid{0000-0003-4637-3059}\IsCorresp} \address{Department of Mathematical Sciences University of Delaware, Newark DE 19716, USA.} \email[P. Monk]{monk@udel.edu} \thanks{The research of the second author was partially supported by AFOSR grant FA9550-23-1-0256.} \subjclass{35Q60} \makeatletter \renewcommand*\cdr@preabstracthook{\vspace*{-5pt}} \def\@setafterauthor{% \vglue3mm% % \hspace*{0pt}% \begingroup\hsize=11cm\advance\hsize\abstractmarginL\raggedright \noindent %\hspace*{\abstractmarginL}\begin{minipage}[t]{10cm} \leftskip\abstractmarginL \normalfont\Small \@afterauthor\par \endgroup \vskip2pt plus 3pt minus 1pt } \makeatother \begin{abstract} Modal expansions are useful to understand wave propagation in an infinite closed electromagnetic waveguide. They can also be used to construct generalized Dirichlet-to-Neumann maps that provide artificial boundary conditions for truncating a computational domain when discretizing the field by finite elements. The modes of a waveguide arise as eigenfunctions of a non-symmetric eigenvalue problem, and the eigenvalues determine the propagation (or decay) of the modes along the waveguide. For the successful use of waveguide modes, it is necessary to know that the modes exist and form a dense set in a suitable function space containing the trace of the electric field in the waveguide. This paper is devoted to proving such a density result using the methods of Keldysh. We also show that the modes satisfy a useful orthogonality property, and show how the Dirichlet-to-Neumann map can be calculated. Our existence and density results are proved under realistic regularity assumptions on the cross section of the waveguide, and the electromagnetic properties of the materials in the waveguide, so generalizing existing results. \end{abstract} \begin{altabstract} Les expansions modales sont utiles pour comprendre la propagation des ondes dans un guide d'ondes électromagnétiques fermé et infini. Elles peuvent également être utilisées pour construire des cartes de Dirichlet à Neumann généralisées qui peuvent être utilisées pour fournir des conditions limites artificielles afin de tronquer un domaine de calcul lors de la discrétisation du champ par des éléments finis. Les modes d'un guide d'ondes apparaissent comme des fonctions propres d'un problème de valeurs propres non symétrique, et les valeurs propres déterminent la propagation (ou la décroissance) des modes le long du guide d'ondes. Pour une utilisation réussie des modes de guide d'ondes, il est nécessaire de savoir que les modes existent et forment un ensemble dense dans un espace de fonctions approprié contenant la trace du champ électrique dans le guide d'ondes. Cet article est consacré à la démonstration d'un tel résultat de densité en utilisant les méthodes de Keldysh. Nous montrons également que les modes satisfont une propriété d'orthogonalité utile, et nous montrons comment la carte de Dirichlet à Neumann peut être calculée. Nos résultats d'existence et de densité sont prouvés sous des hypothèses de régularité réalistes sur la section transversale du guide d'ondes et les propriétés électromagnétiques des matériaux dans le guide d'ondes, généralisant ainsi les résultats existants. \end{altabstract} \begin{document} \maketitle %\vspace*{5pt} \section{Introduction} Waveguides are devices used to transmit electromagnetic waves, for example connecting antenna elements in a device. There are several different types of transmission line including coaxial waveguides, microstrips, and strip lines~\cite{Pozar:11}. This paper focuses on closed waveguides. Out of the many possible designs, we shall focus on ones which have the following characteristics: they are long compared to the wavelength of the electromagnetic field with an invariant shape along the intended direction of propagation (assumed to be along the $z$-coordinate axis in this paper) and they are surrounded by a Perfectly Electrically Conducting (PEC or metallic) surface~\cite[Page~96]{Pozar:11}. They may involve magnetic or dielectric components. For these structures the modes of the waveguide are associated to a propagation constant $\beta$ and corresponding angular frequency $\omega$ of the radiation (see the next section for a precise definition of these terms). The presence of variable coefficients and perfect conductors in waveguides may produce solutions to Maxwell's equations with low regularity. We assume $\omega>0$ to be fixed and $\beta$ to be the eigenvalue. We also allow the presence of thin PEC structures (producing fields with corner singularities) as well as allowing both electric permittivity and magnetic permeability to vary. Our goal is then to prove the density of the eigen-modes in an appropriate space. To prove this density result, we apply~\cite[Theorem~2.1]{Halla:22MSteklovStab} (which in turn traces back to the work of Keldysh and is a reformulation of~\cite[Chapter~V, Theorem~8.3]{GohbergKrein:69} and~\cite[Theorem~4.2, Corollary~3.3]{Markus:88}) to a variational formulation of the modal eigenvalue problem suggested by Lee et al.~\cite{Lee91} and~\cite[Section~8.2]{Jin_book}. A related formulation is given by Vardepetyan\linebreak \& Demkowicz in~\cite{VardapetyanDemkowicz:03}, and we comment on the equivalence of these formulations of the problem. Such a density result is needed to give a theoretical underpinning to modal methods for solving waveguide problems. In addition, if we wish to compute the solution of a scattering problem in the waveguide in which there is a locally perturbed geometry or structure (for example a change in diameter of the waveguide, a corner or junction or a flaw) it is necessary to truncate the domain to allow finite elements to be used in a bounded region containing the perturbation, while a suitable truncation condition is used to account for the infinite translationally invariant parts of the structure. One way to do this is to use the electromagnetic version of the Dirichlet-to-Neumann map on a cross section of the waveguide. For this to be successful we need to know that, for a fixed $\omega$, the set of modes corresponding to all propagation constants are dense in the appropriate boundary space. Both Lee's method and the method of Vardepetyan \& Demkowicz can be used numerically. Numerical results for Lee et al.'s method can be found in~\cite{Lee91}. Vardepetyan \& Demkowicz~\cite{VardapetyanDemkowicz:03} prove that their method gives a convergent approximation of the the propagation constants, and numerical tests can be found in~\cite{VardapetyanDemkowiczNeikirk:03}. An alternative approach is to eliminate the four fields transverse to the direction of the transmission line. This has been studied for the problem of computing modes in a waveguide with piecewise constant electromagnetic properties~\cite{SheSmi13,SheSmi13a}. Assuming the electromagnetic coefficients are piecewise constant, an eigenvalue problem for a 4th order operator pencil can be derived. It is proved that eigenvalues exist~\cite{SheSmi13} and that only finitely many are real, and in~\cite{SheSmi13a} completeness properties are also derived. It seems to be critical for this approach that the coefficients are piecewise constant, whereas the approach in~\cite{Lee91,VardapetyanDemkowicz:03} allows for general piecewise smooth (precisely piecewise $W^{1,\infty}$) coefficients, so that our results extend the completeness theory of~\cite{SheSmi13,SheSmi13a} to allow for graded materials. The remainder of the paper proceeds as follows. In Section~\ref{deriv} we start by giving a derivation of Lee and Vardapetyan \& Demkowicz's eigenvalue problems from Maxwell's equations, and comment on conditions when they are equivalent. We then proceed to prove our main theorem, the existence and density result, in Section~\ref{density}. Then in Section~\ref{DtN} we show how the modes can be used to compute a DtN operator suitable for discretization. Concerning notation, we shall denote by $\Omega\subset\setR^2$ the cross section of the transmission line enclosed by a bounded PEC surface. Additional assumptions on $\Omega$ will be formulated in the theorems we state. We denote by $\bfnu$ the two dimensional outward normal to $\Omega$ and by $\bftv$ the unit tangential vector to the domain on the boundary $\partial\Omega$. For square integrable scalar or vector functions $f$ and $g$ we use the notation \[ \left\langle f,g\right\rangle=\int_{\Omega}f\cdot\overline{g}\,dx. \] The same notation is also used for duality pairings. \section{Derivation of the modal eigenvalue problem}\label{deriv} We now provide a brief derivation of the eigenvalue problems of Lee et al.~\cite{Lee91} and Vardepetyan \& Demkowicz~\cite{VardapetyanDemkowicz:03}. The axis of the waveguide is assumed to be parallel to the $z$ coordinate axis and we let $\Omega\subset \setR^2$ denote the cross-section of the waveguide in the $(x,y)$ plane. Critically the magnetic permeability $\mu$ and the electric permeability $\epsilon$ of the material in the waveguide are assumed to be independent of $z$. Usually $\mu$ is constant, but in case of magnetic components in the device it may differ from unity. For simplicity, we shall assume that $\mu$ is a scalar quantity. Both $\epsilon$ and $\mu$ are assumed to be real valued, bounded and uniformly positive scalar functions of position (precisely they are assumed to be piecewise $W^{1,\infty}$ functions). Let $\hat\bfx=(x,y)$ and $\bfx=(x,y,z)$. We search for propagating modes along the axis of the waveguide having a given angular frequency $\omega>0$. In particular if $(\calE,\calH)$ denotes the three-dimensional time harmonic electric and magnetic fields in the device, we assume that \begin{equation}\label{waveguide_form} \calE(\bfx)=\bfE(\hat\bfx)\exp(i\beta z)\quad\text{ and }\quad\calH(\bfx)=\bfH(\hat\bfx)\exp(i\beta z) \end{equation} where $\beta$ is the aforementioned propagation constant. Note that if ${\calU}(\bfx)=\bfU(\hat{\bfx})\exp(i\beta z)$ for some $\bfU=(U_1, U_2, U_3)^T\in \setR^3$ then \[ \nabla \times{\calU}=\left( \begin{array}{c}{U}_{3,y}-i\beta {U}_2\\-{U}_{3,x}+i\beta {U}_1\\{U}_{2,x}-{U}_{1,y} \end{array}\right)\exp(i\beta z):= \exp(i\beta z)\nabla_\beta \times \calU \] So the second order Maxwell system for $\calE$ (having eliminated the magnetic field $\calH$ in the usual way) can be rewritten as \begin{equation}\label{beta1} \nabla_\beta\times\mu^{-1}\nabla_\beta \times \bfE-\omega^2\epsilon\bfE=0. \end{equation} We now set \[ \hboldE=\left( \begin{array}{c} {E}_1\\{E}_2 \end{array}\right) \] and recall the standard scalar and vector curl in the plane given by \[ \curl\,\hat\boldw=w_{2,x}-w_{1,y}\quad\text{ and }\quad\Curl \,w=\left( \begin{array}{c}w_y\\-w_x \end{array}\right) \] where $\hat{\boldw}=(w_1,w_2)^T$ and $w$ is a scalar function. Then~\eqref{beta1} can be rewritten, using respectively the first two components, and then the third component, as \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\omega^2\epsilon\hboldE+\frac{\beta^2}{\mu}\hboldE+\frac{i\beta}{\mu}\nabla{E}_3 &=0.\label{eq:VD1} \\ \div\frac{1}{\mu}\nabla {E}_3+\omega^2\epsilon {E}_3-i\beta\div\frac{1}{\mu}\hboldE &=0.\label{eq:VD2} \end{align} The fact that the 3D field is divergence free gives \begin{equation}\label{eq:VD3} \div\epsilon\hboldE+i\beta \epsilon {E}_3=0. \end{equation} The above equations hold in $\Omega$. As mentioned in the introduction, we assume that the domain is enclosed in a perfectly conducting (PEC) boundary. In addition there may be PEC surfaces within the domain representing thin metallic components. On the PEC surface and the internal surfaces, we have $\nu_1{E}_2-\nu_2{E}_1=0$ where $(\nu_1,\nu_2)^T$ is the unit normal on these edges. In addition, for the PEC surfaces we know that ${E}_3=0$ so that the above equations are supplemented by the boundary conditions \begin{equation}\label{eq:VDBC} \bftv\cdot\hboldE=0\text{ and } {E}_3=0\text{ on }\partial \Omega, \end{equation} where $\boldsymbol{\tau}=(-\nu_2,\nu_1)^T$. Equations~\eqref{eq:VD2} and~\eqref{eq:VD3} are redundant. The formulation of the eigenvalue problem due to Lee~\cite{Lee91} is to use~\eqref{eq:VD1} and~\eqref{eq:VD2}, so seeking $\beta\in \setC$ and non-trivial $(\hboldE,E_3)\in H_0(\curl;\Omega)\times H^1_0(\Omega)$ such that \begin{subequations}\label{eq:systemVD2} \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\omega^2\epsilon\hboldE+\frac{\beta^2}{\mu}\hboldE+\frac{i\beta}{\mu}\nabla{E}_3&=0 \text{ in }\Omega,\\ \div\frac{1}{\mu}\nabla {E}_3+\omega^2\epsilon {E}_3-i\beta\div\frac{1}{\mu}\hboldE&=0\text{ in }\Omega.\label{eq:systemVD2b} \end{align} \end{subequations} This formulation will be used to prove the density of the eigenmodes in Section~\ref{density}. The eigenvalue problem that is the main focus of~\cite{VardapetyanDemkowicz:03} uses~\eqref{eq:VD1} and~\eqref{eq:VD3} and seeks $\beta\in \setC$ and non-trivial $(\boldE,E_3)\in H_0(\curl;\Omega)\times H^1_0(\Omega)$ such that \begin{subequations}\label{eq:systemVD1} \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\omega^2\epsilon\hboldE+\frac{\beta^2}{\mu}\hboldE+\frac{i\beta}{\mu}\nabla{E}_3&=0 \text{ in }\Omega,\label{eq:systemVD1a} \\ \div\epsilon\hboldE+i\beta \epsilon {E}_3&=0 \text{ in }\Omega.\label{eq:systemVD1b} \end{align} \end{subequations} together with the boundary conditions~\eqref{eq:VDBC}. The numerical analysis of this problem is the subject of~\cite{VardapetyanDemkowicz:03}, and this is the version they advocate for numerical purposes. In the discussion of the modal eigenvalue problem the cutoff frequencies sometimes require a special treatment. These are frequencies $\omega$ for which there exists a vanishing eigenvalue $\beta_j(\omega)=0$. Those cutoff frequencies are exactly the square roots of the eigenvalues $\lambda_n$ of the Dirichlet problems \begin{subequations}\label{eq:cutoff} \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\lambda\epsilon\hboldE&=0 \text{ in }\Omega,\\ \div\frac{1}{\mu}\nabla {E}_3+\lambda\epsilon {E}_3&=0\text{ in }\Omega. \end{align} \end{subequations} We note that the two formulations define the same modes provided $\beta\not=0$, so our theory also applies to the formulation of Vardapetyan \& Demkowicz. In particular, since $\omega>0$, if in addition we are not at a cutoff wave number so that $\beta\not=0$, equations~\eqref{eq:systemVD2} follow from~\eqref{eq:systemVD1} as can be seen by taking the divergence of~\eqref{eq:systemVD1a} and using the result to replace $\div\epsilon \hboldE$ in~\eqref{eq:systemVD1a} to obtain~\eqref{eq:systemVD2b}. Proceeding similarly, when $\beta\not=0$, we see that~\eqref{eq:systemVD2} implies~\eqref{eq:systemVD1}. Note that if $(\beta,(\hboldE,E_3))$ is an eigenpair, then, assuming $\epsilon$ and $\mu$ are real, so are $(\ol{\beta},(\ol{\hboldE},-\ol{E}_3))$ and $(-\beta,(\hboldE,-E_3))$. Hence the eigenvalues are symmetric about both the real and the imaginary axes. To obtain a linear eigenvalue problem, it is suggested in~\cite{Lee91,Jin_book,VardapetyanDemkowicz:03} to define $\tilde{E}_3=i\beta E_3$ so that system formulation due to Lee et al~\eqref{eq:systemVD2} becomes: seek $\beta\in \setC$ and non-trivial $(\hboldE,\tilde{E}_3)\in H_0(\curl;\Omega)\times H^1_0(\Omega)$ such that \begin{subequations}\label{eq:systemVD2m} \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\omega^2\epsilon\hboldE+\frac{\beta^2}{\mu}\hboldE+\frac{1}{\mu}\nabla{\tilde{E}}_3&=0 \text{ in }\Omega,\\ \div\frac{1}{\mu}\nabla {\tilde{E}}_3+\omega^2\epsilon \tilde{E}_3+\beta^2\div\frac{1}{\mu}\hboldE&=0,\text{ in }\Omega,\label{eq:systemVD2bm} \end{align} \end{subequations} together with the boundary conditions~\eqref{eq:VDBC}. In the same way, \eqref{eq:systemVD1} becomes \begin{subequations}\label{eq:systemVD1m} \begin{align} \Curl\frac{1}{\mu}\curl \hboldE-\omega^2\epsilon\hboldE+\frac{\beta^2}{\mu}\hboldE+\frac{1}{\mu}\nabla{\tilde{E}}_3&=0 \text{ in }\Omega,\label{eq:systemVD1ma} \\ \div\epsilon\hboldE+ \epsilon {\tilde{E}}_3&=0 \text{ in }\Omega,\label{eq:systemVD1mb} \end{align} \end{subequations} together with the boundary conditions~\eqref{eq:VDBC}. Now the two approaches define the same eigenvalues and eigenfunctions even if $\beta=0$. Note that in addition there exist several other possible ways to reformulation the eigenvalue problem~\eqref{eq:systemVD2m}. The main reason why we use~\eqref{eq:systemVD2m} is that after some technical computations it can be interpreted as a sufficiently weak perturbation of the linear eigenvalue problem for a selfadjoint operator. \section{Density of eigenmodes}\label{density} This section is devoted to analyzing the eigenmodes using~\eqref{eq:systemVD2} together with the boundary conditions~\eqref{eq:VDBC}. To this end lets recall that for Hilbert spaces $\calX,$ and $\calY$ a compact operator $K\in L(\calX, \calY)$ is in Schatten class $K_p(\calX,\calY)$ of order $p \in (0,\infty)$, if the sequence of singular values $s_n(K),$ $n\in \setN$, of $K$ is $\ell^p(\setN)$ summable. Then, a compact operator $K \in L(\calX,\calY)$ is said to be of finite order if there exists $p \in (0,\infty)$ such that $K\in K_p(\calX,\calY)$. The technique applied in this section is based on a perturbation analysis of selfadjoint operators due to Keldysh~\cite[Appendix]{Markus:88} which we use in the form as presented in~\cite[Theorem~2.2]{Halla:22MSteklovStab}. \begin{theo}[Keldysh]\label{thm:Keldysh} Let $X$ be an infinite dimensional separable Hilbert space and $L(X)$ be the space of bounded operators from $X$ to $X$. Let $T\in L(X)$ be compact and $K\in L(X)$ be injective, selfadjoint and compact. Let $T$ or $K$ be of finite order. Then the spectrum of the operator function $A(\lambda)=I-T-\lambda K$ consists of an infinite sequence of eigenvalues, each with finite algebraic multiplicity, which do not accumulate in $\setC$ and such that for every $\delta>0$ only finitely many eigenvalues lie outside the sectors $\{z\in\setC\colon |\arg z|<\delta\}$ and $\{z\in\setC\colon |\arg z-\pi|<\delta\}$. If $K$ is nonnegative, then for every $\delta>0$ only finitely many eigenvalues lie outside the sector $\{z\in\setC\colon |\arg z|<\delta\}$. Further the closure of the space spanned by the generalized eigenspaces of $A(\cdot)$ is dense in $X$. \end{theo} The second ingredient is a result on the Schatten class of Sobolev embeddings. \begin{theo}[{\cite[Theorem~2.2]{Halla:22MSteklovStab}}]\label{thm:SobolevFO} Let $n\in\setN$ and $\Omega\subset\setR^n$ be a bounded Lipschitz domain. Let $s_2>s_1\geq0$ and $k>\max\{s_2,n/2\}$. Then the embedding $E_{H^{s_2},H^{s_1}}$ from $H^{s_2}(\Omega)$ into $H^{s_1}(\Omega)$ is compact and in $K_{2k/(s_2-s_1)}(H^{s_2}(\Omega),H^{s_1}(\Omega))$. \end{theo} \begin{lemm}[{\cite[Lemma~2.3]{Halla:22MSteklovStab}}]\label{lem:SCproduct} Let $W_1, W_2$ and $W_3$ be separable Hilbert spaces. Let $B\in L(W_1,W_2)$ and $A\in L(W_2,W_3)$. If $s>0$ and either $B\in\SC_\indSC(W_1,W_2)$ or $A\in$ $\SC_\indSC(W_2,W_3)$, then $AB\in\SC_\indSC(W_1,W_3)$. and \begin{align*} \|AB\|_{\SC_\indSC(W_1,W_3)} &\leq \|A\|_{L(W_1,W_2)} \|B\|_{\SC_\indSC(W_1,W_2)}\\ \intertext{or} \|AB\|_{\SC_\indSC(W_1,W_3)} &\leq \|A\|_{\SC_\indSC(W_1,W_2)} \|B\|_{L(W_1,W_2)} \end{align*} respectively. \end{lemm} We assume that the equation $-\div\mu^{-1}\nabla E_3-\omega^2\epsilon E_3=0$ in $\Omega$, $E_3=0$ on $\partial\Omega$ is uniquely solvable which holds if we exclude those frequencies $\omega$ that support the existence of cutoff wave numbers $\beta=0$. To perform this analysis we define \[ X=H_0\left(\curl,\left(\div\mu^{-1}\right)^0;\Omega\right)\times\nabla H^1_0(\Omega) \] where, recalling that $\boldsymbol{\tau}$ denotes the unit tangent on $\partial \Omega$, \[ H_0\left(\curl,\left(\div\mu^{-1}\right)^0;\Omega\right)= \left\{\bfu\in H({\curl};\Omega)\;\middle|\; \boldsymbol{\tau}\cdot \bfu=0\text{ on }\partial \Omega\text{ and } \div\left(\mu^{-1}\bfu\right)=0\text{ in }\Omega\right\}. \] We will use the following scalar product on this space \[ \left\langle (\bfu,w),\left(\bfu',w'\right)\right\rangle_X=\left\langle \mu^{-1}\curl\boldu,\curl\boldu'\right\rangle +\left\langle\epsilon\nabla w,\nabla w'\right\rangle +\left\langle\epsilon\boldu,\boldu'\right\rangle, \] which simplifies some of the analysis. This section is devoted to proving the main theorem of this paper: \begin{theo}\label{th1} Assume that $\omega$ is not a cutoff frequency (i.e.\ $\omega^2$ is not an eigenvalue of the Dirichlet problem~\eqref{eq:cutoff}). Then the spectrum of the eigenvalue problem~\eqref{eq:systemVD2m} with the boundary conditions~\eqref{eq:VDBC} consists of an infinite sequence of eigenvalues, each with finite algebraic multiplicity, which do not accumulate in $\setC$ and such that for every $\delta>0$ only finitely many eigenvalues lie outside the sectors $\{z\in \setC\;|\; |\arg z-\pi|<\delta\}$ and $\{z \in\setC\;|\; |\arg z+\pi |<\delta\}$. Finally the closure of the space spanned by the $\hboldE$ components of the generalized eigenspaces of the problem is dense in $H_0(\curl;\Omega)$. \end{theo} \begin{proof} We start by building the Schur complement with respect to $E_3$ (this is done at the discrete, matrix, level in~\cite[Section~8.2]{Jin_book}) and obtain \begin{subequations}\label{eq:evp-Schur} \begin{align} \Curl\mu^{-1}\curl \hboldE - \omega^2\epsilon\hboldE+\beta^2\mu^{-1}\hboldE\label{Schur1} +\beta^2\mu^{-1}\nabla S \div\mu^{-1}\hboldE &=0,\text{ in }\Omega,\\ \tv\cdot\hboldE&=0,\text{ on }\partial\Omega,\label{Schur2} \end{align} \end{subequations} where formally $S:=(-\div\mu^{-1}\nabla-\omega^2\epsilon)^{-1}\in L(H^{-1}(\Omega),H^1_0(\Omega))$. Next we use the orthogonal decomposition $H_0(\curl;\Omega)=H_0(\curl,(\div\mu^{-1})^0;\Omega)\oplus^\bot \nabla H^1_0(\Omega)$ and write $\hboldE=\boldu+\nabla w$. Using this expansion, we obtain the variational equation for $(\bfu,w)\in X$: \begin{multline*} \left\langle \mu^{-1}\curl\boldu,\curl\boldu'\right\rangle-\omega^2\left\langle\epsilon\boldu,\boldu'\right\rangle+\beta^2\left\langle\mu^{-1}\boldu,\boldu'\right\rangle-\omega^2\big(\left\langle\epsilon\boldu,\nabla w'\right\rangle+\left\langle\epsilon\nabla w,\boldu'\right\rangle\big)\\ -\omega^2 \left\langle\epsilon\nabla w,\nabla w'\right\rangle+\beta^2 \left\langle\mu^{-1}\nabla w,\nabla w'\right\rangle +\beta^2 \left\langle \mu^{-1}\nabla S\div\mu^{-1}\nabla w,\nabla w'\right\rangle=0 \end{multline*} for all $(\bfu',w')\in X$. We compute, using the definition of $S$ \begin{align*} \left\langle \mu^{-1}\nabla S\div\mu^{-1}\nabla w,\nabla w'\right\rangle &=\left\langle \mu^{-1}\nabla S(\div\mu^{-1}\nabla w+\omega^2\epsilon w),\nabla w'\right\rangle -\omega^2\left\langle \mu^{-1}\nabla S(\epsilon w),\nabla w'\right\rangle \\ &=-\left\langle \mu^{-1}\nabla w,\nabla w'\right\rangle -\omega^2\left\langle \mu^{-1}\nabla S(\epsilon w),\nabla w'\right\rangle \end{align*} Then using integration by parts and the definition of $S$ again: \begin{align*} \left\langle \mu^{-1}\nabla S\div\mu^{-1}\nabla w,\nabla w'\right\rangle &=-\left\langle \mu^{-1}\nabla w,\nabla w'\right\rangle +\omega^2\left\langle \epsilon w,S\div\mu^{-1}\nabla w'\right\rangle\\ &=-\left\langle \mu^{-1}\nabla w,\nabla w'\right\rangle +\omega^2\left\langle \epsilon w,S\left(\div\mu^{-1}\nabla w'+\omega^2\epsilon w'\right)\right\rangle -\omega^4\left\langle \epsilon w,S \left(\epsilon w'\right)\right\rangle\\ &=-\left\langle \mu^{-1}\nabla w,\nabla w'\right\rangle -\omega^2\left\langle \epsilon w, w'\right\rangle -\omega^4\left\langle \epsilon w,S \left(\epsilon w'\right)\right\rangle. \end{align*} Hence $(\bfu,w)\in X$ satisfies \begin{multline*} \left\langle \mu^{-1}\curl\boldu,\curl\boldu'\right\rangle-\omega^2\left\langle\epsilon\boldu,\boldu'\right\rangle+\beta^2\left\langle\mu^{-1}\boldu,\boldu'\right\rangle -\omega^2\big(\left\langle\epsilon\boldu,\nabla w'\right\rangle+\left\langle\epsilon\nabla w,\boldu'\right\rangle\big)\\ -\omega^2 \left\langle\epsilon\nabla w,\nabla w'\right\rangle -\beta^2 \omega^2\big(\left\langle \epsilon w, w'\right\rangle +\omega^2\left\langle S(\epsilon w), \epsilon w'\right\rangle\big)=0 \end{multline*} for all $(\bfu',w')\in X$. We switch the sign of the test function $w'\to-w'$ and obtain \begin{equation}\label{opeq} (A-T-\lambda K)(\boldu,w)=0 \end{equation} with $\lambda:=-\beta^2$ and \begin{align*} \left\langle A(\boldu,w),\left(\boldu',w'\right)\right\rangle&:= \left\langle \mu^{-1}\curl\boldu,\curl\boldu'\right\rangle +\left\langle\epsilon\nabla w,\nabla w'\right\rangle +\left\langle\epsilon\boldu,\boldu'\right\rangle,\\ \left\langle T(\boldu,w),\left(\boldu',w'\right)\right\rangle&:= \omega^2\big(\left\langle\epsilon\boldu,\boldu'\right\rangle-\left\langle\epsilon\boldu,\nabla w'\right\rangle+\left\langle\epsilon\nabla w,\boldu'\right\rangle\big) +\left\langle\epsilon\boldu,\boldu'\right\rangle,\\ \left\langle K(\boldu,w),\left(\boldu',w'\right)\right\rangle&:= \left\langle\mu^{-1}\boldu,\boldu'\right\rangle -\omega^2\big(\left\langle \epsilon w, w'\right\rangle +\omega^2\left\langle S(\epsilon w), \epsilon w'\right\rangle\big). \end{align*} In view of the definition of the scalar product $\langle\cdot,\cdot\rangle_X$ we see that $A=I$. Both $T$ and $K$ are compact and of finite order. In particular, consider the embedding operator $E\in L(H_0(\curl,(\div\mu^{-1})^0;\Omega),\boldL^2(\Omega))$, the multiplication operator $M_\epsilon\in L(\boldL^2(\Omega))$ and the gradient operator $G\in L(H^1_0(\Omega),\boldL^2(\Omega))$, $Gw:=\nabla$. Thence $T=\omega^2(E^*M_\epsilon E- G^*M_\epsilon E+ E^*M_\epsilon G) + E^* M_\epsilon E$. Due to~\cite{Ciarlet:20} there exists $s>0$ such that $H_0(\curl,(\div\mu^{-1})^0;\Omega)$ embeds continuously into $\boldH^s(\Omega)$. Thus it follows with~\cite[Theorem~2.2, Lemma~2.3]{Halla:22MSteklovStab} that $E$ and hence $T$ is of finite order. With the same approach it can be seen that $K$ is of finite order too. In addition $K$ is selfadjoint. The operator $K$ is also injective. Indeed let $(\boldu,w)\in\ker K$. It easily follows that $\boldu=0$ and $w+\omega^2S(\epsilon w)=0$. Multiplying with $S^{-1}$ we obtain \begin{align*} 0=-\div\mu^{-1}\nabla w-\omega^2\epsilon w+\omega^2\epsilon w=-\div\mu^{-1}\nabla w, \end{align*} and thus $w=0$. Now we can simply apply~\cite[Theorem~2.1]{Halla:22MSteklovStab}. It remains to note that for an eigenpair $(\lambda,(\boldu,w))$ of~\eqref{opeq}, $(\pm\sqrt{-\lambda},\boldu+\nabla w)$ are eigenpairs of~\eqref{eq:evp-Schur} and vice-versa. \end{proof} Note that here we do not have to deal with the involved concept of double completeness, because the eigenfunctions for $\pm\beta$ are identical. Furthermore, for topologically nontrivial domains the operator $K$ remains injective. Note also that we could use this formulation also for the computation of modes. We can enforce the constraint $\div\mu^{-1}\boldu=0$ weakly by the introduction of a scalar auxiliary variable. The implementation of $S$ requires also the introduction of a scalar auxiliary variable. Thus we would obtain a problem in the space $H_0(\curl;\Omega)\times H^1_0(\Omega)\times H^1_0(\Omega)$, which is more expensive than using the space $H_0(\curl;\Omega)\times H^1_0(\Omega)$ which is used for~\eqref{eq:systemVD1}. \section{A transparent boundary condition}\label{DtN} Before we start the main derivation of this section we prove the following lemma. \begin{lemm}\label{lem:orthog} Assume that the equation $-\div\mu^{-1}\nabla E_3-\omega^2\epsilon E_3=0$ in $\Omega$, $E_3=0$ on $\partial\Omega$ is uniquely solvable. In addition, suppose $(\beta_j,(\hboldE_j,E_{3,j}))$ and $(\beta_k,(\hboldE_k,E_{3,k}))$ are both eigenpairs of~\eqref{eq:systemVD2m}. Assume $\beta_j^2\not=\overline{\beta}_k^2$, $\beta_j\not=0$, $\beta_k\not=0$, then \begin{align*} \left\langle \mu^{-1}\curl\hboldE_j,\curl\hboldE_k\right\rangle-\omega^2\left\langle\epsilon\hboldE_j,\hboldE_k\right\rangle=0. \end{align*} \end{lemm} \begin{proof} Using the Schur complement derived in the proof of the previous theorem (see~\eqref{Schur1}-\eqref{Schur2}) we see that the variational formulation of~\eqref{eq:systemVD2m} is to find $\beta$ and $\hboldE\in H_0(\curl;\Omega)$ with $\hboldE\not=0$ such that \begin{equation}\label{eq:junk} \left\langle \mu^{-1}\curl\hboldE,\curl \xi\right\rangle-\omega^2\left\langle \epsilon \hboldE,\xi\right\rangle=-\beta^2\left(\left\langle \mu^{-1}\hboldE,\xi\right\rangle-\left\langle S\left(\div\frac{1}{\mu}\hboldE\right),\div \frac{1}{\mu}\xi\right\rangle\right) \end{equation} for all $\xi\in H_0(\curl;\Omega)$. Equation~\eqref{eq:junk} is of the kind $A\hboldE=\beta^2B\hboldE$ with self-adjoint operators $A,B$. Hence $\langle A\hboldE_k,\hboldE_j \rangle = \beta_k^2 \langle B\hboldE_k,\hboldE_j \rangle$ and also $\langle A\hboldE_k,\hboldE_j \rangle = \langle \hboldE_k,A\hboldE_j \rangle = \ol{\beta_j}^2 \langle \hboldE_k,B\hboldE_j \rangle =\ol{\beta_j}^2 \langle B\hboldE_k,\hboldE_j \rangle$. Thus \[ \left(\beta_k^{-2}-\ol{\beta_j}^{-2}\right) \left\langle A\hboldE_k,\hboldE_j \right\rangle = 0. \] If $(\beta_k^{-2}-\ol{\beta_j}^{-2})$ does not vanish, $\langle A\hboldE_k,\hboldE_j \rangle$ has to vanish. \end{proof} \begin{figure} \centering \includegraphics[scale=0.3]{Terminate.jpg} \caption{Cartoon of a terminated wave guide. The semi-infinite waveguide $D^{\rext}$ (blue region) is terminated by the bounded domain $D^{\rint}$ (yellow region). If finite elements are used to discretize the termination region $D^{\rint}$, we can provide a boundary condition on the artificial boundary $\Omega$ by using a Dirichlet-to-Neumann map constructed from the modes propagating in the waveguide.}\label{Fig1} \end{figure} Now we can show how the eigen-modes can be used to provide a termination condition for waveguide calculations. When simulating waveguides by finite elements it is necessary to truncate the unbounded computational domain. One way to do this is using a transparent boundary condition that can be computed using the modes we have analyzed in the previous Section~\cite{Kim:17}. A cartoon example is shown in Figure~\ref{Fig1} where a vertical waveguide denoted $D^{\rext}$ is terminated by a bounded domain $D^{int}$. We assume a source $\calJ$ is located in the bounded domain $D^{\rint}$ (having compact support) and let $(\epsilon,\mu)$ denote the electromagnetic parameters of the three dimensional domain. In $D^{\rext}$, the functions $\mu$ and $\epsilon$ are assumed to satisfy the assumptions at the start of Section~\ref{deriv} so that $D^{\rext}$ is a waveguide and we can use the modes analyzed in the previous section to represent the solution on $\Omega$ (due to their verified density). In $D^{\rint}$ the electromagnetic parameters are positive real piecewise $W^{1,\infty}$ functions. Using finite element elements the termination region $D^{\rint}$ can be discretized using finite elements. Further up the waveguide we need to terminate the computational domain by cutting the waveguide at $\Omega$ which for simplicity we assume is in the plane $z=0$. Then the total electric field $\calE^{\rint}$ (now a general vector function not of the form~\eqref{waveguide_form}) satisfies \[ \nabla_3\times \mu^{-1}\nabla_3\times \calE^{\rint}-\omega^2\epsilon\calE^{\rint}=i\omega \calJ\quad\text{in}\quad D^{\rint}\subset \setR^3, \] together with the PEC boundary condition on the external boundary of $D^{\rint}$. Here $\nabla_3\times$ is the full three dimensional curl. Multiplying by a smooth test vector function $\calU'$ and integrating by parts, we obtain the equality: \[ \int_{D^{\rint}}\mu^{-1}\nabla_3\times \calE^{\rint}\cdot\nabla_3\times \overline{\calU'}-\omega^2\epsilon\calE^{\rint}\cdot\overline{\calU'}\,dV+\int_\Omega \bfe_3\times \mu^{-1}\nabla_3\times \calE^{\rint}\cdot\overline{\calU'}\,dA= \int_{D^{\rint}}i\omega \calJ\cdot\overline{\calU'}\,dV, \] where $\bfe_3=(0,0,1)^T$. By the continuity of the tangential trace of the electric and magnetic fields across $\Omega$ at $z=0$, \begin{align} \bfe_3\times\left(\calE^{\rext}\times\bfe_3\right)&=\bfe_3\times\left(\calE^{\rint}\times\bfe_3\right)\label{tanjE} \\ \bfe_3\times \mu^{-1}\nabla_3\times \calE^{\rext}&= \bfe_3\times \mu^{-1}\nabla_3\times \calE^{\rint}.\label{tanjH} \end{align} In addition, in the domain $D^{\rext}=\Omega\times (0,\infty)$, the field $\calE^\rext$ satisfies Maxwell's equations and is outgoing. By the assumed uniqueness of this problem in the semi-infinite domain we may define the analogue of the Dirichlet-to-Neumann map $\bfT$ such that \begin{align*} \bfT\left(\bfe_3\times\left(\calE^{\rext}\times\bfe_3\right)\right) &=\bfe_3\times \mu^{-1}\nabla_3\times \calE^{\rext} \\ \intertext{so that} \int_\Omega \bfe_3\times \mu^{-1}\nabla_3\times \calE^{\rint}\cdot\overline{\calU'}\,dA &=\int_\Omega\bfT\left(\bfe_3\times\left(\calE^{\rint}\times\bfe_3\right)\right)\cdot\overline{\calU'}\,dA. \end{align*} Thus we obtain a variational problem for $\calE^{\rint}\in\bfX$ that must satisfy \[ \int_{D^{\rint}}\mu^{-1}\nabla_3\times \calE^{\rint}\cdot\nabla_3\times \overline{\calU'}-\omega^2\epsilon\calE^{\rint}\cdot\overline{\calU'}\,dV+\int_\Omega \bfT\left(\bfe_3\times\left(\calE^{\rint}\times\bfe_3\right)\right)\cdot\overline{\calU'}\,dA= \int_{D^{\rint}}i\omega \calJ\cdot\overline{\calU'}\,dV, \] for all $\calU\in \bfX$ where \[ \bfX=\left\{\bfu\in H\left(\rm{curl};D^{\rint}\right)\;\middle|\;\bfnu\times(\bfu\times\bfnu)=0\text{ on }\partial D^{\rint}\setminus \Omega\right\} \] and $\bfnu$ is the unit outward normal to $D^{\rint}$. This problem is suitable for finite element discretization provided we can calculate an at least an accurate approximation to $\bfT$. We show how to do this next. In $D^{\rext}$ we set \begin{align*} \calE^\rext(\boldx,z)=\left( \begin{array}{c} \hboldE^\rext(\boldx,z)\\ E_3^\rext(\boldx,z) \end{array}\right). \end{align*} By~\eqref{tanjE} and the density of the modes on $\Omega$ we can approximate \begin{equation}\label{exp} \bfe_3\times\left(\calE^\rint(\bfx,0)\times\bfe_3\right)\approx \sum_{j=1}^N a_j\hboldE_j(\bfx)\text{ in } H_0({\curl};\Omega) \end{equation} where we can achieve any desired accuracy by taking $N$ large enough. The equality of the traces above and below $\Omega$ implies that \[ \calE^\rext(\bfx,z)\approx \sum_{j=1}^N a_j\left( \begin{array}{c} \hboldE_j(\bfx)\\ E_{3,j}(\bfx) \end{array}\right)\exp\left(i\beta_j z\right) \] where $\beta_j$ is chosen such that $\beta_j>0$ if $\beta_j$ is real, and $\Im(\beta_j)>0$ if $\beta_j$ has a non-zero imaginary part in order to obtain an outgoing or evanescent mode. Then, as we have seen, we encounter \begin{align*} \int_{\Omega}\bfT\left(\bfe_3\times\left(\calE^\rext\times\bfe_3\right)\right)\cdot\overline{\calU}'\,dA &=\int_{\Omega}\bfe_3\times\mu^{-1}\nabla_{3}\times\calE^\rext\cdot\overline{\calU}'\,dA\\ &=\left\langle \bfe_3\times\mu^{-1} \left( \begin{array}{c} E^\rext_{3,y}-E^\rext_{2,z}\\ -E^\rext_{3,x} + E^{\rint}_{1,z}\\ E^\rext_{2,x}-E^\rext_{1,y} \end{array}\right), \calU'\right\rangle\\ &=\left\langle \mu^{-1} \left( \begin{array}{c}E^\rext_{3,x}- E^\rext_{1,z}\\E^\rext_{3,y}-E^\rext_{2,z}\\0 \end{array}\right), \calU'\right\rangle =\left\langle\mu^{-1} \bpm \nabla E^{\rext}_3-\frac{\partial}{\partial z}\hboldE^{\rext}\\0\epm,\calU'\right\rangle, \end{align*} where we used the choce of $\bfe_3$ to compute the cross product. To construct a transparent boundary condition by means of a Dirichlet-to-Neumann operator we need to derive a meaningful expression for the right hand side above. For a mode of the cross sectional problem $(\hboldE,E_3)$, with propagation constant $\beta$ assuming that $\beta\not=0$, we note that the corresponding wave-guide mode is $(\hboldE,E_3)\exp(i\beta z)$ so for this mode we compute \begin{align*} \mu^{-1} \left( \begin{array}{c}E_{3,x}- E_{1,z}\\E_{3,y}-E_{2,z}\\0 \end{array}\right) &=\mu^{-1} \bpm \nabla E_3-i\beta\hboldE \\ 0 \epm\exp\left(i\beta z\right)\\ &=\mu^{-1} \bpm \frac{\mu}{i\beta}\left(-\Curl\frac{1}{\mu}\curl \hboldE+\omega^2\epsilon\hboldE-\frac{\beta^2}{\mu}\hboldE)\right)-i\beta\hboldE \\ 0 \epm\exp\left(i\beta z\right)\\ &= i\beta^{-1} \bpm \Curl\mu^{-1}\curl\hboldE-\omega^2\epsilon\hboldE \\ 0 \epm\exp\left(i\beta z\right). \end{align*} Note that $\calU'$ has a vanishing tangential trace on the waveguide boundary, and the tangential trace $\hboldE^\rint$ of $\calE^{\rint}$ on the interface $\Omega$ has a vanishing tangential trace on $\partial\Omega$. Thus we can integrate the right hand side of the former expression by parts and obtain, for a single mode, \begin{align*} \left\langle\mu^{-1} \left( \begin{array}{c}E_{3,x}- E_{1,z}\\E_{3,y}-E_{2,z}\\0 \end{array}\right),\calU'\right\rangle=i\beta^{-1} \left\langle \mu^{-1}\curl\hboldE,\curl\overline{\calU'}\right\rangle-\omega^2\left\langle\epsilon\hboldE,\overline{\calU'}\right\rangle. \end{align*} Assuming that $\beta=0$ is not an eigenvalue, we can normalize the eigenfunctions such that \begin{align*} \left\langle \mu^{-1}\curl\hboldE_j,\curl\hboldE_j\right\rangle-\omega^2\left\langle\epsilon\hboldE_j,\hboldE_j\right\rangle :=\gamma_j=\pm 1\text{ or }0. \end{align*} We then need to compute an expansion of the field in terms of modes on $\Omega$ (i.e. find the coefficients $a_n$ in the expansion~\eqref{exp}). To see how this can be done in a special case, suppose that $\beta_j^2$, $j=1,\cdots,\infty$ are real and the eigenvalues are simple. Then the propagation constants are either real or purely imaginary (i.e. no real part), so using Lemma~\ref{lem:orthog} the projection of $\hboldE^{\rint}$ onto mode $\hboldE_j$ is given by \begin{align*} \gamma_j\left(\left\langle \mu^{-1}\curl\hboldE^{\rint},\curl\hboldE_j\right\rangle-\omega^2\left\langle\epsilon\hboldE^{\rint},\hboldE_j\right\rangle\right)\hboldE_j. \end{align*} Thus, in this simple case, to construct an (approximate) DtN-operator, we have the following steps \begin{enumerate} \item Choose $N$ and compute the modes $(\beta_j,\boldE_j)$, ($j=1,\,\dots,\,N$). \item Compute the boundary term \begin{align*} \int_{\Omega}&\bfT(\bfe_3\times(\calE^\rext\times\bfe_3))\cdot\overline{\calU}'\,dA=\\ &\sum_j -i\beta_j^{-1} \gamma_j\left(\left\langle \mu^{-1}\curl\hboldE^{\rint},\curl\hboldE_j\right\rangle-\omega^2\left\langle\epsilon\hboldE^{\rint},\hboldE_j\right\rangle\right) \left(\left\langle \mu^{-1}\curl\hboldE_j,\curl\calU'\right\rangle-\omega^2\left\langle\epsilon\hboldE_j,\calU'\right\rangle\right). \end{align*} \end{enumerate} \section{Conclusion} Modal expansions are used to model closed waveguides. We have shown the most general denseness result to date for these modes, and this result allows for discontinuous electromagnetic parameters, geometric singularities in the solution and graded materials (spatially varying electromagnetic parameters). We have also shown a type of orthogonality result which my be used to simplify calculations with these modes. \section{Declaration of interests} The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations. \printbibliography \end{document}