On spurious solutions encountered in Helmholtz scattering resonance computations in R^d with applications to nano-photonics and acoustics

TL;DR

This paper introduces a new low-cost sorting scheme based on volume integral equations to identify potentially spurious scattering resonances in electromagnetic and acoustic problems, applicable to complex material structures.

Contribution

A novel sorting method using Lippmann-Schwinger volume integral equations for detecting spurious solutions in scattering resonance computations.

Findings

Effective identification of spurious solutions in test cases

Applicable to graded and piecewise constant materials

Works for electromagnetic and acoustic resonances

Abstract

In this paper, we consider a sorting scheme for potentially spurious scattering resonant pairs in one- and two-dimensional electromagnetic problems and in three-dimensional acoustic problems. The novel sorting scheme is based on a Lippmann-Schwinger type of volume integral equation and can, therefore, be applied to structures with graded materials as well as to configurations including piece-wise constant material properties. For TM/TE polarized electromagnetic waves and for acoustic waves, we compute first approximations of scattering resonances with finite elements. Then, we apply the novel sorting scheme to the computed eigenpairs and use it to mark potentially spurious solutions in electromagnetic and acoustic scattering resonances computations at a low computational cost. Several test cases with Drude-Lorentz dielectric resonators as well as with graded material properties are considered.

Tables

Table 1. Table 1: Memory consumption estimation for matrices in ( 29 ), and ( 26 ) for d = 1 , 2 , 3 𝑑 1 2 3 d=1,2,3 .

$d$	$p$	$m$	$N_K$	$N_z^S$	$N_z^D$	$W^S\left(bytes\right)$	$W^D\left(bytes\right)$
$1$	$2$	$10$	${10}^1$	$9.0\times {10}^1$	$4.0\times {10}^2$	$1.4\times {10}^3$	$6.4\times {10}^3$
$2$	$2$	$10$	${10}^2$	$8.1\times {10}^3$	$1.6\times {10}^5$	$1.3\times {10}^5$	$2.6\times {10}^6$
$3$	$2$	$10$	${10}^3$	$7.3\times {10}^5$	$6.4\times {10}^7$	$1.2\times {10}^7$	$1.0\times {10}^9$
$1$	$2$	${10}^2$	${10}^2$	$9.0\times {10}^2$	$4.0\times {10}^4$	$1.4\times {10}^4$	$6.4\times {10}^5$
$2$	$2$	${10}^2$	${10}^4$	$8.1\times {10}^5$	$1.6\times {10}^9$	$1.3\times {10}^7$	$2.5\times {10}^{10}$
$3$	$2$	${10}^2$	${10}^6$	$7.3\times {10}^8$	$6.4\times {10}^{13}$	$1.2\times {10}^{10}$	$1.0\times {10}^{15}$
$1$	$2$	${10}^3$	${10}^3$	$9.0\times {10}^3$	$4.0\times {10}^6$	$1.4\times {10}^5$	$6.4\times {10}^7$
$2$	$2$	${10}^3$	${10}^6$	$8.1\times {10}^7$	$1.6\times {10}^{13}$	$1.3\times {10}^9$	$2.6\times {10}^{14}$
$3$	$2$	${10}^3$	${10}^9$	$7.3\times {10}^{11}$	$6.4\times {10}^{19}$	$1.2\times {10}^{13}$	$1.0\times {10}^{21}$

Table 2. Table 2: Cost and time estimation for computing ( 42 ) for d = 1 , 2 , 3 𝑑 1 2 3 d=1,2,3 .

$d$	$p$	$m$	$N_r$	cost	$t_i\left(s\right)$	$t\left(s\right)$
$1$	$2$	${10}^1$	${10}^1$	$1.8\times {10}^3$	$9.0\times {10}^{-5}$	$1.8\times {10}^{-3}$
$2$	$2$	${10}^1$	${10}^2$	$3.2\times {10}^6$	$8.1\times {10}^{-3}$	$3.2\times {10}^0$
$3$	$2$	${10}^1$	${10}^3$	$5.8\times {10}^9$	$7.3\times {10}^{-1}$	$5.8\times {10}^3$
$1$	$2$	${10}^2$	${10}^2$	$1.8\times {10}^5$	$9.0\times {10}^{-4}$	$1.8\times {10}^{-1}$
$2$	$2$	${10}^2$	${10}^4$	$3.2\times {10}^{10}$	$8.1\times {10}^{-1}$	$3.2\times {10}^4$
$3$	$2$	${10}^2$	${10}^6$	$5.8\times {10}^{15}$	$7.3\times {10}^2$	$5.8\times {10}^9$
$1$	$2$	${10}^3$	${10}^3$	$1.8\times {10}^7$	$9.0\times {10}^{-2}$	$1.8\times {10}^2$
$2$	$2$	${10}^3$	${10}^6$	$3.2\times {10}^{14}$	$8.1\times {10}^1$	$3.2\times {10}^8$
$3$	$2$	${10}^3$	${10}^9$	$5.8\times {10}^{21}$	$7.3\times {10}^5$	$5.8\times {10}^{15}$

Table 3. Table 3: Drude Lorentz data for Gold, taken from [ 52 ] , with time convention e − i ​ ω ​ t superscript 𝑒 𝑖 𝜔 𝑡 e^{-i\omega t} .

${ϵ}_{\infty }=1$	${\omega }_p=$ 9.03	-
$f_0=$ 0.76	${\omega }_0=$ 0	${\gamma }_0=$ 0.053
$f_1=$ 0.024	${\omega }_1=$ 0.415	${\gamma }_1=$ 0.241
$f_2=$ 0.01	${\omega }_2=$ 0.83	${\gamma }_2=$ 0.345
$f_3=$ 0.071	${\omega }_3=$ 2.969	${\gamma }_3=$ 0.87
$f_4=$ 0.601	${\omega }_4=$ 4.304	${\gamma }_4=$ 2.494
$f_5=$ 4.384	${\omega }_5=$13.32	${\gamma }_5=$ 2.214

Equations

$$-\Delta E_3-{\omega }^2ϵE_3=0\text{and}-\nabla \cdot \left(\frac{1}{ϵ}\nabla H_3\right)-{\omega }^2{H}_3=0,$$

$$T(\omega )u:=u-K(\omega )u=0$$

$$\begin{array}[]{l}T\!M\!:\,K(\omega)u(x):=\omega^{2}\int_{\Omega_{r}}\Phi(x,y)\left(\epsilon(y)-1\right)u(y)\,dy\\[5.69054pt] T\!E:\,\,K(\omega)u(x):=\nabla\cdot\int_{\Omega_{r}}\Phi(x,y)\left(\tfrac{1}{\epsilon(y)}-1\right)\nabla u(y)\,dy\end{array},\,\,\Phi(x,y):=\left\{\begin{array}[]{ll}\frac{i}{2\omega}e^{i\omega|x-y|},&d=1\\[5.69054pt] \tfrac{i}{4}H_{0}^{(1)}(\omega|x-y|),&d=2.\end{array}\right.$$

$$-\Delta u-\frac{{\omega }^2}{c(x{)}^2}u=0,$$

$$K(\omega )u(x):={\omega }^2{\int }_{{\Omega }_r}\Phi (x,y)\left(\frac{1}{c(y{)}^2}-1\right)u(y)dy,\Phi (x,y)=\frac{e^{i\omega |x-y|}}{4\pi |x-y|}.$$

$$\Vert (A-{\omega }^2)u\Vert \le \delta \Vert u\Vert$$

$${\sigma }_{\delta }(T)=\sigma (T)\cup \{\omega \in \rho (T):\Vert T^{-1}(\omega )\Vert >{\delta }^{-1}\}.$$

$${ϵ}_{metal}(\omega ):={ϵ}_{\infty }+\sum _{j=0}^{N_p}\frac{f_j{\omega }_p^2}{{\omega }_j^2-{\omega }^2-i\omega {\gamma }_j},$$

$$ϵ(x,\omega ):=\sum _{m=0}^{N_r}{ϵ}_m(\omega ){\chi }_{{\Omega }_m}(x),x\in {\Omega }_a,\omega \in \mathcal{D},$$

$$\begin{array}[]{rcl}-\nabla\cdot\left(\rho\nabla u\right)-\omega^{2}\eta u&=&0,\end{array}$$

$$\mathfrak{a}(\omega )[u,v]:={\int }_{{\Omega }_a}\rho \nabla u\cdot \nabla \overline{v}\mathrm{dx},\mathfrak{b}(\omega )[u,v]:={\int }_{{\Omega }_a}\eta u\overline{v}\mathrm{dx},$$

$$\frac{\partial u}{\partial n}=\mathcal{G}(\omega )u\text{on }\Gamma :=\partial {\Omega }_a,$$

$$-{\left(\rho u^{\prime }\right)}^{\prime }-{\omega }^2\eta u=0forx\in {\Omega }_a,$$

$$u^{\prime }(-a)=-i\omega u(-a),u^{\prime }(a)=i\omega u(a).$$

$$g_1(\omega )[u,v]:=i\omega (u(a)\overline{v}(a)+u(-a)\overline{v}(-a)).$$

$$q_1(\omega )[u,v]:=\mathfrak{a}(\omega )[u,v]-{\omega }^2\mathfrak{b}(\omega )[u,v]-g_1(\omega )[u,v]=0,$$

$$\mathcal{G}(\omega )u:=\frac{1}{2\pi }\sum _{\nu =-\infty }^{\infty }\omega \frac{H_{\nu }^{(1)\prime }(\omega a)}{H_{\nu }^{(1)}(\omega a)}{e}^{i\nu \theta }{\int }_0^{2\pi }u(a,{\theta }^{\prime })e^{-i\nu {\theta }^{\prime }}d{\theta }^{\prime }$$

$$q(\omega )[u,v]:=\mathfrak{a}(\omega )[u,v]-{\omega }^2\mathfrak{b}(\omega )[u,v]-g(\omega )[u,v]=0,$$

$$g(\omega )[u,v]:=({\mathcal{G}}_{{\nu }_{\max }}(\omega )u,v{)}_{{\Gamma }_a}=\sum _{\nu =-{\nu }_{\max }}^{{\nu }_{\max }}\omega a\frac{H_{\nu }^{(1)\prime }(\omega a)}{H_{\nu }^{(1)}(\omega a)}{\hat{u}}_{\nu }{\bar{\hat{v}}}_{\nu },{\hat{\phi }}_{\nu }=\frac{1}{\sqrt{2\pi }}{\int }_0^{2\pi }\phi (a,\theta )e^{-i\nu \theta }d\theta .$$

$$\tilde{\sigma}(r):=\left\{\begin{array}[]{ll}0,&\hbox{if}\,\,r<a\\ P(r),&\hbox{if}\,\,a\leq r\leq b\\ \sigma_{0},&\hbox{if}\,\,r>b\end{array}\right.$$

$$\frac{d}{d\tilde{x}}=\frac{1}{\tilde{\alpha }(x)}\frac{d}{dx},with\tilde{\alpha }(x)=1+i\tilde{\sigma }(x).$$

$$-\frac{d}{dx}\left(\frac{\rho }{\tilde{\alpha }}\frac{du}{dx}\right)-{\omega }^2\eta \tilde{\alpha }u=0,x\in {\Omega }_{\ell },u(r_{\ell })=0andu(-r_{\ell })=0.$$

$$t_1(\omega )[u,v]:=\mathfrak{a}(\omega )[u,v]-{\omega }^2\mathfrak{b}(\omega )[u,v]+{\hat{t}}_1(\omega )[u,v]=0,$$

$$\begin{array}[]{lllll}\tilde{\alpha}(r)\!\!\!&:=1+i\tilde{\sigma}(r),&&\tilde{r}(r)\!\!\!&:=(1+i\tilde{\sigma})r=\tilde{\alpha}(r)\,r,\\[5.69054pt] \sigma(r)\!\!\!&:=\tilde{\sigma}(r)+r\displaystyle\frac{\partial\tilde{\sigma}}{\partial r},&&\alpha(r)\!\!\!&:=\displaystyle\frac{\partial\tilde{r}}{\partial r}=1+i\sigma(r),\end{array}$$

$$t(\omega )[u,v]:=\mathfrak{a}(\omega )[u,v]-{\omega }^2\mathfrak{b}(\omega )[u,v]+\hat{t}(\omega )[u,v]=0,$$

$$\mathcal{A}=\left(\begin{array}[]{cc}\tfrac{\tilde{\alpha}}{\alpha}\cos^{2}\theta+\tfrac{\alpha}{\tilde{\alpha}}\sin^{2}\theta&\bigg{(}\tfrac{\tilde{\alpha}}{\alpha}-\tfrac{\alpha}{\tilde{\alpha}}\bigg{)}\sin\theta\cos\theta\\ \bigg{(}\tfrac{\tilde{\alpha}}{\alpha}-\tfrac{\alpha}{\tilde{\alpha}}\bigg{)}\sin\theta\cos\theta&\tfrac{\tilde{\alpha}}{\alpha}\sin^{2}\theta+\tfrac{\alpha}{\tilde{\alpha}}\cos^{2}\theta\end{array}\right),\quad\mathcal{B}:=\alpha\tilde{\alpha}.$$

$$\begin{array}[]{ll}&T^{\gamma}(\omega^{\gamma})\xi=(I-K(\omega^{\gamma}))\xi=0,\,\,\,\hbox{with}\\[8.53581pt] T\!M:&K_{ij}(\omega):=\omega^{2}\int_{\Omega_{r}}\Phi(x_{i},y)\left(\epsilon(y)-1\right)\varphi_{j}(y)\,dy,\\[5.69054pt] T\!E:&K_{ij}(\omega):=\nabla\cdot\int_{\Omega_{r}}\Phi(x_{i},y)\left(\tfrac{1}{\epsilon(y)}-1\right)\nabla\varphi_{j}(y)\,dy.\end{array}$$

$$S^{\gamma }(\Omega ):=\{u\in H^1({\Omega }_a):{u|}_{K_j}\in {\mathcal{P}}_p(K_j)for{K}_j\in \mathcal{T}\}.$$

$$u^{\gamma }=\sum _{j=1}^N{\xi }_j{\phi }_j,A_{ij}=(\rho \nabla {\phi }_j,\nabla {\phi }_i{)}_{{\Omega }_a},M_{ij}=(\eta {\phi }_j,{\phi }_i{)}_{{\Omega }_a}.$$

$$F(\omega )\xi :=\left(A-{\omega }^{2}M+Q\right)(\omega )\xi =0,$$

Full text

††journal: JOURNAL

1 Introduction

The most common approach to approximate scattering resonances is to truncate the domain with a perfectly matched layer (PML) and discretize the differential equations with a finite element method. This results in approximations of the true resonances but in practice also a large number of solutions that are unrelated to the true resonances. One approach to lessen the problem with nonphysical solutions in finite element (FE) computations and decrease the number of spurious solutions is to design an appropriate $hp$-refinement strategy [1]. However, for large problems spurious eigenvalues remains a major problem in resonance computations [2] and in quasimodal expansions [3, 4].

The identification of spurious pairs from resonance computations with PML has been attempted using a sensitivity approach [5, 6, 2]. The approach is based on the observation that spurious eigenvalues are sensitive to parameter perturbations, while well-converged resonances are not. However, it has been observed that this approach is able to identify approximations to true resonances only in small regions of the complex plane and only if the finite element space is large. Additionally, the method is for a fixed discretization expensive as it requires the computation of all eigenpairs in a spectral window several times, from where re-meshing and re-assembling of the FE must be performed.

The origin of nonphysical eigenvalues in scattering resonance computations is spectral instability, which is common for non-normal operators [7, 8]. Spectral instability is known to be less problematic with volume integral equations compared with formulations based on differential operators. Therefore, we proposed in [9] to use a Lippmann-Schwinger type of integral equation for marking potentially spurious solutions in a one-dimensional setting. However, the direct application of the ideas presented in [9] to higher dimensions turns out to be more challenging and computationally demanding.

In this paper, we develop a scheme for marking potentially spurious solutions in higher dimensions when the computational domain is truncated by a PML or a Dirichlet-to-Neumann (DtN) map. In particular, we show that the developed scheme is computationally cheap and provides valuable information on the location of potentially spurious eigenvalues. For example, this information is useful in quasimodal expansions of the field and when adaptive finite element methods are employed.

In the new sorting strategy, a particular approximation of a Lippmann-Schwinger based residual is computed for each eigenpair. The sorting strategy is computationally efficient since (i) the computation re-uses the pre-computed FE environment (ii) the testing has low memory requirements, and (iii) the algorithm is fully parallelizable. The method is successfully applied to metal-dielectric nanostructures in ${\mathbb{R}}^{2}$, where the metal is modeled by a sum of Drude-Lorentz terms, and to an acoustic benchmark problem in ${\mathbb{R}}^{3}$. Several additional benchmarks confirm the efficiency of the proposed sorting scheme in ${\mathbb{R}}^{d}$ with $d=1,2,3$.

2 Electromagnetic and acoustic scattering resonances

Assume that $\epsilon(x,\omega)=\epsilon(x_{1},x_{2},\omega)$ is independent of $x_{3}$ and consider electromagnetic waves propagating in the $(x_{1},x_{2})$-plane. The $x_{3}$-independent electromagnetic field $({\bf E},{\bf H})$ is then decomposed into transverse electric (TE) polarized waves $(E_{1},E_{2},0,0,0,H_{3})$ and transverse magnetic (TM) polarized waves $(0,0,E_{3},H_{1},H_{2},0)$ [10]. This decomposition reduces Maxwell’s equations to one scalar equation for $H_{3}$ and one scalar equation for $E_{3}$. The TM-polarized waves and the TE-polarized waves satisfy formally

[TABLE]

respectively. For the scattering resonance problems, $E_{3}$ and $H_{3}$ are assumed to be locally $L_{2}$-integrable functions that satisfy an outgoing condition [11, 12].

Let the physical domain $\Omega_{a}\subset\displaystyle\mathbb{R}^{d}$ be an open ball of radius $a$ with boundary $\Gamma_{a}$, and let $\Omega_{r}:=\text{supp}\,(\epsilon-1)\subset\Omega_{a}$ be the bounded domain defining the resonators. Hence, we assume that the relative permittivity $\epsilon$ for $x\in{\mathbb{R}}^{d}\setminus\Omega_{r}$ is one. Furthermore, let $\Omega_{r}:=\cup_{i=1}^{N}\Omega_{i}$ denote the union of disjoint resonators $\Omega_{1},\,\Omega_{2},...,\Omega_{N}$ as shown in Figure 1. A scattering resonance was in [5] formally defined as a complex number $\omega$ for which the Lippmann-Schwinger equation

[TABLE]

has a non-zero solution $u$.

The integral operator $K$ in (2) is for TM/TE waves given by

[TABLE]

Here, $\Phi(x,y)$ is known as the outgoing Green function in free space for the Helmholtz equation [11, 12]. Notice that while we are interested in $x\in\Omega_{a}$, the integration in (2) is only performed over $\Omega_{r}$, since the integration over the air region $\Omega_{0}:=\Omega_{a}\setminus\Omega_{r}$ is zero.

The scattering resonance problem (2) is a highly non-linear eigenvalue problem, where the matrices after discretization are large and full. It is possible to accurately solve the non-linear eigenvalue problem $T(\omega)u=0$ in ${\mathbb{R}}$ using a standard laptop; see e.g. [13, 9]. However, accurate computations of eigenpairs of (2) in higher dimensions would require huge computer resources; See [14] and the discussion in Section 7.5.

Acoustic scattering resonances in ${\mathbb{R}}^{3}$

Sound-soft materials are characterized by the speed of sound $c(x)$ and we assume that acoustic resonators are defined by $\Omega_{r}:=\text{supp}\,(c^{-2}-1)\subset\Omega_{a}$. Then, the acoustic pressure $u$ satisfies formally the Helmholtz equation

[TABLE]

where the outgoing condition can be expressed as $u$ satisfying an expansion in spherical harmonics outside the open ball $\Omega_{a}$ [15].

Moreover, $(\omega,u)$ is a scattering resonance pair if (2) holds with

[TABLE]

Note that the acoustic problem in ${\mathbb{R}}^{d}$, $d=1,2,3$ is analogous to TM-polarized electromagnetic waves with $\epsilon=1/c^{2}$.

2.1 Alternative formulations

Understanding the resonance behavior of structures in unbounded domains is important and many different approaches have been proposed. Graded material properties are increasingly popular in applications [16] and we will therefore not consider boundary integral equations. However, boundary integral equation based methods are a good alternative for cases with piecewise constant coefficients and not too complicated geometries [17]. The most popular method to compute resonances in ${\mathbb{R}}^{d}$, $d>1$ is the finite element (FE) method with a perfectly matched layer (PML). In recent years, finite element methods based on Hardy space infinite elements (HIF) [18] and DtN maps [19] have also been proposed as strong alternatives to compute resonances in higher dimensions. For the DtN map, recent developments in computational linear algebra are a key to the high performance of the method [20, 21]. Discretization with FE of the PML and HIF formulations result in sparse matrices, and a formulation in terms of a DtN result in sparse matrices except a small dense block corresponding to the DtN map. Moreover, the PML and HIF formulations result in a standard generalized eigenvalue problem if $\epsilon$ is $\omega$-independent and in the general case the non-linearity in $\omega$ is completely determined by $\epsilon(x,\omega)$. Hence, the PML and HIF formulations seem to have the most attractive properties of the considered methods. However, it is very important to also take into account the so called spectral instability. Then, the picture changes completely, as discussed in the next section.

2.2 Spectral instability and pseudospectra

Let $A$ denote an unbounded closed linear operator in a Hilbert space with domain $\text{dom}\,A$, spectrum $\sigma(A)$, and resolvent set $\rho(A)$. Then $A$ exhibits high spectral instability if for a very small $\delta>0$ there exist many $\omega^{2}\in{\mathbb{C}}$ and $u\in\text{dom}\,A$ such that

[TABLE]

even though $\omega^{2}$ is not close to $\sigma(A)$ [22]. This is closely related to the pseudospectrum $\sigma_{\delta}(A)$ which is defined as the union of $\sigma(A)$ and all $\omega^{2}$ in the resolvent set $\rho(A)$ for which it exists an $u\in\text{dom}\,A$ such that (6) holds. The generalization of those results to an operator function $T$ is straightforward and we will in some of the numerical computations rely on the following alternative characterization of the pseudospectrum:

[TABLE]

It is well known that PML and HIF based methods encounter high spectral instability [18, 9]. Methods based on a DtN map encounter medium spectral instability [19, 9] and integral equation based methods encounter low spectral instability [9]. Hence, the Lippmann-Schwinger equation is in our setting the preferred method in terms of spectral stability. This will be further discussed in the paper.

2.3 Domain and material properties

In optics, the material properties of non-magnetic metals are characterized by the complex relative permittivity function $\epsilon$, which changes rapidly at optical frequencies $\omega$. The most common accurate material model is then the Drude-Lorentz model

[TABLE]

where $\epsilon_{\infty}\geq 1$ and $f_{j}$, $\omega_{p},\,\omega_{j}$, $\gamma_{j}$ are non-negative [10]. Hence, the Maxwell eigenvalue and scattering resonance problems in $\omega$ are nonlinear for metal-dielectric nanostructures. Assume that the domain of the resonators can be written in the form $\Omega_{r}:=\cup_{i=1}^{N}\Omega_{i}$ and let $\chi_{\Omega_{m}}$ denote the characteristic function of the subset $\Omega_{m}$. For material properties that are piecewise constant in $\Omega_{a}$, we assume a permittivity function in the form

[TABLE]

where the dependencies on $\omega\in\mathcal{D}\subset\displaystyle\mathbb{C}$ in $\epsilon_{m}$ for $m=0,1,\dots$ are of Drude-Lorentz type (7). In addition, we will consider graded material properties, meaning that $\epsilon$ is a continuous function in $x$. In linear acoustics, the speed of sound $c$ is assumed to be independent of the frequency.

3 DtN and PML based methods

In the next sections, we will describe two common approaches to compute scattering resonances and the restriction of resonance modes to a compact subset of ${\mathbb{R}}^{d}$. In the following, we use the notation

[TABLE]

where $u:=E_{3},\,\rho:=1,\,\eta:=\epsilon$ for the TM-case and $u:=H_{3},\,\rho:=1/\epsilon,\,\eta:=1$ for the TE-case.

We define for $u,v\in H^{1}(\Omega_{a})$ the forms

[TABLE]

where in (10), $\rho$ and $\eta$ are functions of $\omega\in\mathcal{D}$. Let $\mathcal{Z}$ denote the set of values $\omega$ that are zeros or poles of $\epsilon$ and set $\mathcal{D}:=\displaystyle\mathbb{C}\setminus\mathcal{Z}$.

3.1 DtN based methods

Scattering resonances $\omega$ and quasi-normal modes $u$ restricted to $\Omega_{a}$ can be determined from a problem that utilizes a Dirichlet-to-Neumann (DtN) map operator denoted $\mathcal{G}(\omega)$ [23, 24, 19]. Below we present variational formulations for ${\mathbb{R}}^{d}$, $d=1,2$. Formally, $(\omega,u)$ is a scattering resonance pair if (9) holds in $\Omega_{a}$ and

[TABLE]

where $\partial u/\partial n$ is the normal derivative.

3.1.1 DtN formulation in ${\mathbb{R}}$

In one space dimension the scattering resonance problem restricted to $\Omega_{a}:=(-a,a)$ is formally: Find a non-zero $u$ and a complex $\omega$ such that

[TABLE]

where the DtN-map at $x=\pm a$ is

[TABLE]

Define for $u,v\in H^{1}(\Omega_{a})$ and $\omega\in\mathcal{D}\subset{\mathbb{C}}$ the forms $\mathfrak{a},\,\mathfrak{b}$ as in (10), and

[TABLE]

The nonlinear eigenvalue problem is then as follows: Find $u\in H^{1}(\Omega_{a})\backslash\{0\}$ and $\omega\in\mathcal{D}$ satisfying

[TABLE]

for all $v\in H^{1}(\Omega_{a})$. Note that (15) is a quadratic eigenvalue problem if $\epsilon$ is independent of $\omega$ and a rational eigenvalue problem for Drude-Lorentz type of materials (7).

3.1.2 DtN formulation in ${\mathbb{R}}^{2}$

In this subsection, we present a DtN formulation in polar coordinates $(r,\theta)$. Let $H^{(1)}_{\nu}(z)$ denote the Hankel function of first kind, then the DtN operator (11) on the circle $\Gamma_{a}$ has the explicit form

[TABLE]

and $\mathcal{G}(\omega):H^{1/2}(\Gamma_{a})\rightarrow H^{-1/2}(\Gamma_{a})$ is bounded [23].

The resonance problem restricted to $\Omega_{a}$ is formally to find non-trivial solutions $(\omega,u)$ such that (9) and (11) with (16) holds. The theory presented in [23] can with minor changes be used in the present case to derive properties of a variational formulation of the problem.

Variational formulation: Let $S$ denote the union of the set of zeros of $H^{(1)}_{\nu}(\omega a),\,\nu\in\displaystyle\mathbb{Z}$, and $\mathcal{G}_{\nu_{\rm max}}(\omega)$ denote the operator (11) truncated after $|\nu|={\nu_{\rm max}}$. The eigenvalues of the truncated version of (9)-(11) are determined by the following variational problem: Find $u\in H^{1}(\Omega_{a})\setminus\{0\}$ and $\omega\in\mathcal{D}:=\displaystyle\mathbb{C}\setminus\{\mathcal{Z}\cup\displaystyle\mathbb{R}^{-}\cup S\}$ such that for all $v\in H^{1}(\Omega_{a})$

[TABLE]

where the forms $\mathfrak{a},\,\mathfrak{b}$ are defined as in (10), and

[TABLE]

3.2 PML based methods

In the previous section, a DtN-map was used to reduce the exterior Helmholtz problem to a bounded domain. In this section, we consider an alternative approach based on a complex coordinate stretching (the PML method), which results in a linear eigenvalue problem for non dispersive material coefficients [25]. The method consists on attaching to $\Omega_{a}$ a buffer layer of thickness $\ell$, where outgoing solutions decay rapidly. The buffer domain is referred to as $\Omega_{{P\!M\!L}}$ and the full computational domain $\Omega:=\Omega_{a}\cup\Omega_{P\!M\!L}$ is enlarged as shown in the Figure 1.

3.2.1 PML formulation in ${\mathbb{R}}$

Let $\ell>0$, and $0<a<b<r_{\ell}$, with $r_{\ell}=b+\ell$. The action of the PML is defined through the stretch function

[TABLE]

with $r=|x|$. The polynomial $P$ in (19) is required to be increasing in $[a,b]$ and $\tilde{\sigma}(r)$ is sufficiently smooth $C^{2}(0,r_{\ell})$. For this we introduce the fifth order polynomial $P$ satisfying: $P(a)=P^{\prime}(a)=P^{\prime\prime}(a)=P^{\prime}(b)=P^{\prime\prime}(b)=0$ and $P(b)=\sigma_{0}$.

The PML problem is restricted to $(-r_{\ell},r_{\ell})$ and the PML strength function has then the profile shown in Figure 1. In the following sections, we consider the transformation rule

[TABLE]

For finite element computations we restrict the domain to $\Omega_{\ell}:=(-r_{\ell},r_{\ell})$ and choose as in [25] homogeneous Dirichlet boundary conditions. Formally, the finite PML problem is then: Find the eigenpairs $(\omega,u)$ such that

[TABLE]

Let $\Omega_{P\!M\!L}=(-r_{\ell},-a)\cup(a,r_{\ell})$ and denote by $\mathfrak{a},\,\mathfrak{b}$ the forms in (10). As variational formulation of (21), we consider:

Find $u\in H_{0}^{1}(\Omega_{\ell})\setminus\{0\}$ and $\omega\in\mathcal{D}:=\displaystyle\mathbb{C}\setminus\mathcal{Z}$ such that for all $v\in H_{0}^{1}(\Omega_{\ell})$

[TABLE]

where $\hat{t}_{1}(\omega)[u,v]=(\tfrac{1}{\tilde{\alpha}}u^{\prime},v^{\prime})_{\Omega_{{P\!M\!L}}}-\omega^{2}(\tilde{\alpha}u,v)_{\Omega_{{P\!M\!L}}}$.

3.2.2 PML formulation for ${\mathbb{R}}^{d}$

Approximation of resonances using a radial PML was analyzed in [25] and we will here only consider the PML problem truncated to the ball $\Omega\subset{\mathbb{R}}^{d},\,d=2,3$. The used complex stretching functions are, as in one-dimension, of the form (19). Let $\mathcal{Z}$ denote the set of values $\omega$ that are zeros or poles of $\epsilon$ and let $\Omega:=\Omega_{a}\cup\Omega_{P\!M\!L}$ denote a partition into the PML-region and the part of the domain $\Omega_{a}$ containing the resonators.

In the sequel we need the following definitions

[TABLE]

with the properties $\sigma(r)=\partial(r\tilde{\sigma})/\partial r$ and $\alpha(r)=\tilde{\alpha}(r)=1+i\sigma_{0}\quad\hbox{for}\quad r>b$.

It is clear that the PML coefficients are designed such that for $r\leq a$, we obtain $\tilde{\sigma}(r)=0$ and $\alpha(r)=1$, which is the same as $\tilde{r}=r$. Hence, the PML operator restricted to $r\leq a$ corresponds to the original operator in problem (9), where there is no PML effect.

Variational formulation: The eigenvalues of (9) are then determined by the following variational problem: Find $u\in H_{0}^{1}(\Omega)\setminus\{0\}$ and $\omega\in\mathcal{D}:=\displaystyle\mathbb{C}\setminus\mathcal{Z}$ such that for all $v\in H_{0}^{1}(\Omega)$

[TABLE]

where $\hat{t}(\omega)[u,v]:=(\mathcal{A}\nabla u,\nabla v)_{\Omega_{P\!M\!L}}-\omega^{2}(\mathcal{B}u,v)_{\Omega_{P\!M\!L}}$, and the forms $\mathfrak{a},\,\mathfrak{b}$ are defined as in (10).

As an example for $d=2$, direct transformation of (24) from polar to Cartesian coordinates results in

[TABLE]

Note that even though $\mathcal{A}$, and $\mathcal{B}$ are defined in $\Omega$, their action takes place only in $\Omega_{P\!M\!L}$.

4 Discretization of the Lippmann-Schwinger equation

In this section we present a collocation method for the discretization of the Lippmann-Schwinger equation (2), which will be used to compute resonances in one-dimension and is the base for the numerical sorting algorithm in Section 6. Further computational details are given in Section 7.5.

4.1 A Galerkin-Nyström method

We present a Galerkin-Nyström discretization method for linear Fredholm integral equations of the second kind. In [26], this method is referred to as case (A) of the Galerkin methods, and convergence for the problem with sources is also discussed. In [5, Sec. 3.2] the method was used for resonance computations for $d=2$.

Let $\{\varphi_{j}\}_{j=1}^{N}$ be piecewise polynomial functions with the property $\varphi_{j}(x_{i})=\delta_{ji}$, $\{x_{i}\}_{i=1}^{N}\in\Omega_{a}$. We introduce the representation $u^{\gamma}=\sum^{N}_{j}\xi_{j}\varphi_{j}$, and with the use of (2) we obtain the nonlinear eigenvalue problem: Find $\xi\in\displaystyle\mathbb{C}^{N}$ and $\omega^{\gamma}\in\displaystyle\mathbb{C}$ such that

[TABLE]

The Nyström method consists in choosing the collocation points $x_{i}$ as the nodes of a high order quadrature rule. By doing this, the convergence of the scheme is considerably improved. In Section 7.5 we describe some of the implementation details of the Galerkin-Nyström discretization.

The resulting nonlinear matrix eigenvalue problem (26) is solved by using a contour integration based method [27, 28, 29].

Remark** 1****.**

The formulation in (2) uses information of the exact solution of the problem at every discretization node $x_{i}$, through its fundamental solution $\Phi(x,y)$. From where a numerical scheme based on (26) is flexible in the way that it can be posed in the smallest domain $\Omega_{r}$, as well as in larger domains $\Omega\supset\Omega_{r}$, without taking any special treatment for boundary conditions.

5 FE discretization of the DtN and PML based formulations

In this section we discuss briefly the details involved in the assembly of the matrices corresponding to the discretization of the formulations given in (9), and (2).

5.1 Discretization with the finite element method

Let the domain $\Omega\subset\displaystyle\mathbb{R}^{d}$ be covered with a regular and quasi uniform finite element mesh $\mathcal{T}(\Omega)$ consisting of elements $\{K_{j}\}^{N_{K}}_{j=1}$. The mesh is designed such that the permittivity function $\epsilon$ is continuous in each $K_{j}$. Let $h_{j}$ be the length of the largest diagonal of the non-curved primitive $K_{j}$ and denote by $h$ the maximum mesh size $h:=\max_{j}{h_{j}}$.

Let $\mathcal{P}_{p}$ denote the space of polynomials on $\displaystyle\mathbb{R}^{d}$ of degree $\leq p$ in each space coordinate and define the $N$ dimensional finite element space

[TABLE]

The computations of discrete resonance pairs $(u^{\gamma},\omega^{\gamma})$ are for $d=2,3$ performed in the approximated domain $\Omega^{\gamma}$ using curvilinear elements [30]. The meshes used are shape regular in the sense of [31, Sec. 4.3], and consist of quadrilateral/brick elements with curvilinear edges/surfaces that deviate slightly from their non-curved primitives.

5.2 Assembly of the FE matrices

In this section we refer to domains $\Omega\in{\mathbb{R}}^{d}$. Let $\{\varphi_{1},\dots,\varphi_{N}\}$ be a basis of $S^{{\gamma}}(\Omega)$. Then $u^{\gamma}\in S^{{\gamma}}(\Omega)$, and the entries in the finite element matrices are of the form

[TABLE]

The matrix eigenvalue problem is then: Find the eigenpars $(\omega,{\xi})\in\mathcal{D}\times{\mathbb{C}}^{N}\setminus\{0\}$ such that

[TABLE]

where the corresponding matrix valued function is

[TABLE]

In the case where $\epsilon(\omega,x)$ is given as piecewise smooth function of space, we write (29) as

[TABLE]

with matrices $\tilde{A}^{m}_{ij}=(\nabla\varphi_{j},\nabla\varphi_{i})_{\Omega_{m}},\,\,\tilde{M}^{m}_{ij}=(\varphi_{j},\varphi_{i})_{\Omega_{m}},\,\,m=0,1,\ldots,N_{r}.$

Remark** 2****.**

Truncation of the DtN:* Let $\left\lceil z\right\rceil$ be the smallest integer greater than or equal to $z$. We use the rule ${\nu_{\rm max}}=\left\lceil a\,\omega_{M}\right\rceil$ according to the results in [19], where from the considered spectral window, $\omega_{M}$ is the largest real part allowed for computations of eigenvalues.*

Remark** 3****.**

Truncation of the PML:* The PML is set up following the discussions in [25, 9], which accounts for large enough $\ell$ and $\sigma_{0}$ such that the search region is feasible. Additionally, we use the space $S_{0}^{{\gamma}}(\Omega):=\{u\in S^{{\gamma}}(\Omega):u=0\,\,\,\hbox{for}\,\,\,x\in\partial\Omega\}$ for computations with the PML formulation.*

Finally, we mention that all discretization methods presented in this work use the same FE space over $\Omega_{a}$.

Remark** 4****.**

All formulations (LS, DtN, PML) use the FE triangulation $\mathcal{T}(\Omega_{a})$, which is the restriction of $\mathcal{T}(\Omega)$ to $\Omega_{a}$.

By using the FE mesh suggested in Remark 4, we ensure that the approximation properties in the physical domain are the same for all formulations.

6 Numerical sorting of potentially spurious solutions

In this section we derive a discrete indicator from equation (2) that allow us to identify potentially spurious solutions once we have computed FE solutions $(\omega^{\gamma}_{m},u^{\gamma}_{m})$ to (29) or to (31). The identification of potentially spurious solutions is important since adaptive finite element methods and quasimodal expansions in practice will fail if spurious solutions are included in the expansion or influence the mesh-refinement. The idea of sorting potentially spurious solutions is in the spirit of the standard residual error estimator and marking strategy used in adaptive finite element methods [32]. Adaptive FE can be applied when the coarsest mesh is sufficiently fine. However, the pre-asymptotic regime is very large in resonance computations, in particular for cases with large air regions in the computational domain. A goal of the paper is therefore to supplement the information given by the PDE based residual estimator with information from an integral equation. The sorting scheme is based on a computationally cheap approximation of the condition $\|\chi_{a}T(\omega)\chi_{a}u\|<\delta$.

Let $\{\varphi_{j}\}$ be a basis for $S^{\gamma}(\Omega_{a})$. Then, the discrete Lippmann-Schwinger equation (26) is written in the form

[TABLE]

Definition 5**.**

Pseudospectrum indicator:* The computed eigenvalue $\omega^{\gamma}$ belongs, for given $\delta>0$, to the $\delta$-psudospectrum $\sigma_{\delta}(T^{\gamma})$ if the pair $(\omega^{\gamma},u^{\gamma})$ satisfies $\|T^{\gamma}(\omega^{\gamma})u^{\gamma}\|_{\Omega_{a}}<\delta$. Then, for a given domain $\Omega\supseteq\Omega_{r}$, we define the pseudospectrum indicator as*

[TABLE]

Particularly, we want to be able to measure whether the computed eigenpair $(\omega^{\gamma},u^{\gamma})$ is converging to a physical pair. Naturally, non convergent pairs exhibit large $\delta^{\gamma}$ values. Additionally, we can identify and remove certain spurious eigenpairs in the PML based formulation by considering the following definition.

Definition 6**.**

PML added eigenpairs:* The use of the coordinate stretching technique in formulations (22), (24) introduces new eigenpairs to problem (9). These new eigenvalues accumulate close to the critical line of the modified PML problem [9], and the corresponding eigenfunctions $v_{m}$ exhibit oscillations in $\Omega_{P\!M\!L}$, but decay in the physical region $\Omega_{a}$. Then, by using the normalization $\|v_{m}\|_{L^{2}(\Omega)}=1$, the PML critical eigenvalues exhibit*

[TABLE]

which can be successfully used as a filtering criterion for removing PML added eigenpairs.

A similar idea was in [2] used to study photoacoustic resonators.

6.1 Numerical pseudospectra computation

Computations of the pseudospectra provide insight into the behavior of the resolvent of the discretized operator function $F^{\gamma}$, allowing us to evaluate its spectral stability. In our computations, we use that $\sigma_{\delta}(F^{\gamma})$ is the set of all $z\in\displaystyle\mathbb{C}$ such that

[TABLE]

where $s_{\scriptsize\hbox{min}}\,F^{\gamma}(z)$ denotes the smallest singular value of $F^{\gamma}(z)$ [7, Def. 2.10]. For the singular value computations we used SLEPc [33].

6.2 On spurious solutions in ${\mathbb{R}}$

In [9] we introduced a sorting strategy based on Definition 5 but performing computations only inside $\Omega_{r}$. As suggested by (32), computations can be also performed in $\Omega_{0}$, and it is natural to ask whether or not the sorting scheme performs worse if air is included in the evaluation. Then, we numerically test the solutions to the Lippmann-Schwinger (LS) formulation by enlarging the computational domain to include air. If the LS operator function would exhibit undesired spurious eigenvalues, this would render the method unreliable for detection of spurious pairs.

The following problem has been considered by several authors including [34, 25]. Define for $n_{1}\neq 1$ the piecewise constant function $n$ as

[TABLE]

The corresponding exact resonances of (12)-(13) for TM polarization are given by

[TABLE]

We start by discussing the stability of the spectrum of the DtN, PML and LS formulations. Namely, the discretization (26) inherits the property described in Remark 1, from where it is expected that the spectrum of the LS formulation would not be sensitive to perturbations of $\Omega_{a}=(-a,a)$ in the sense discussed in Section 2.2. For verification, we propose a similar experiment to the one presented in [9, Figs. 4.1,5.2,5.5], but now we allow larger domains $\Omega_{a}\supseteq\Omega_{r}$ in the computation of the Slab problem (34).

The FE approximations to eigenvalues and the pseudospectrum of the problem are for $a=1,2,3$ obtained using the tools introduced in Section 6.1. We keep track of the number of eigenvalues in a fixed region of the complex plane, and the location of eigenvalues that remain for all perturbations.

As it can be seen in Figure 2 for a coarse mesh, and in Figure 3 with one mesh refinement, the number of eigenvalues in the Lippmann Schwinger formulation remained constant for $a=1,2,3$. This means that no spurious eigenvalues are added due to the inclusion of air. Furthermore, from the plots we see that the computed eigenvalues remain unperturbed when increasing $a$, which shows that the eigenvalues are not sensitive to perturbations of the domain.

In the same Figure we follow the last discussion but for computations using the DtN and PML formulations. It is observed that by increasing $a$, there is an increase of the number of computed eigenvalues in a fixed region of the complex plane. Additionally, it can be seen that the location of eigenvalues is slightly modified when perturbing $a$. These observations lead us to conclude that eigenvalues from the DtN and PML formulations are very sensitive to perturbations of the domain. Note that we for many 1D problems can use any of the methods to obtain very accurate approximations of the resonances without using large computer resources. However, the sensitivity to perturbations is in practice very important for problems in higher dimensions.

Before the next discussion, we briefly introduce quadrature rules for $d=1$, but postpone specific details until Section 7.5. In the one dimensional case, the integration of $f\in C^{2m}(K),\,m\in\displaystyle\mathbb{N}$ over an element $K$ is approximated by formulas of the form $\int_{K}f(x)\,dx=\sum_{i=1}^{m}{\rm w}_{i}f({\rm x}_{i})+E$, where ${\rm w}_{i}$ are quadrature weights, ${\rm x}_{i}$ quadrature nodes, and $E$ is the remainder [35, Ch. 8]. By employing $m$-point quadrature rules with large enough $m$ we can ensure that the remainder $|E|$ is smaller than the FE discretization error. Finally, we gather the quadrature points defined over $\Omega_{0}$ and $\Omega_{r}$, in the sets $\mathcal{X}_{0},\mathcal{X}_{r}$ respectively. Let $\mathcal{X}_{0}:=\cup_{i}\{\mathrm{x}_{j}\in K_{i}|K_{i}\subset\Omega_{0}\}_{j=1}^{m}$, $\mathcal{X}_{r}:=\cup_{i}\{\mathrm{x}_{j}\in K_{i}|K_{i}\subset\Omega_{r}\}_{j=1}^{m}$, and $\mathcal{X}_{a}:=\mathcal{X}_{0}\cup\mathcal{X}_{r}$.

In the discretization of the Lippmann-Schwinger formulation (26), we encounter the situation $x_{i},y\in K_{l}\subset\Omega_{r}$, where the evaluation of the kernel $\Phi(x_{i},y)$ becomes problematic at the evaluation point $y=x_{i}$. In one-dimensional problems ($d=1$), the kernel $\Phi({x},y)$ is continuous, but has a jump in the derivative at points $x=y$. In the troublesome element $K_{l}:=(x_{l},x_{l+1})$, we can always split the integration interval $K_{l}\rightarrow(x_{l},x_{j})\cup(x_{j},x_{l+1})$ and perform two separate quadrature integrations. Then, by using Gauss-type of quadratures, it is possible to avoid the evaluation of $\Phi(x_{j},x_{j})$ [9].

In Figure 4 we show convergence of various quantities of interest for the development of the sorting strategy to be presented in Section 7.6.3. The plots are computed for the slab problem (34) with $p=2$ and featuring five levels of mesh refinement ref $=1,2,\ldots,6$. Horizontally, we show results for the eigenpairs corresponding to $m=1,5,10$ from equation (35). In black solid line we show $|(\omega_{m}-\omega^{\gamma})/\omega_{m}|$ where $\omega^{\gamma}$ corresponds to the closest computed eigenvalue to $\omega_{m}$, and $u^{\gamma}$ its corresponding computed eigenfunction. In blue solid line we show the pseudospectrum indicator $\delta^{\gamma}(\Omega_{a})$, and in dashed lines we show $\min_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ and $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$. Finally, in red circles we add $\max_{{\rm x}\in\mathcal{X}_{0}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$.

From the figure it is evident that as the mesh is refined, there is convergence of the eigenvalue error $|(\omega_{m}-\omega^{\gamma})/\omega_{m}|$, and residuals $\min_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ and $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$. Additionally, the rate of convergence of $\min_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ and $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ is the same as the rate of convergence of $|(\omega_{m}-\omega^{\gamma})/\omega_{m}|$ and the residuals are around two order of magnitude higher than the eigenvalue error. Moreover, the values $\min_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ and $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ are tight bounds for $\delta^{\gamma}(\Omega_{a})$. Finally, the value $\max_{{\rm x}\in\mathcal{X}_{0}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ coincides with $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ for the one-dimensional case.

The fact that the residuals $\min_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ and $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ exhibit the same rate of convergence of the eigenvalue error, can be used as an indication of the eigenvalue convergence for problems where an exact solution is not available.

7 On spurious solutions in ${\mathbb{R}}^{d}$

From Figure 4 it is, for the studied case, evident that the quantity $\max_{{\rm x}\in\mathcal{X}_{a}}|T^{\gamma}(\omega^{\gamma})u^{\gamma}({\rm x})|$ can be used to mark potentially spurious scattering resonance pairs. The proposed sorting strategy for higher dimensions is based on this observation. In the following subsections, we discuss computational details and computational cost before the new sorting indicator is introduced in formula (45).

7.1 Computational details for the evaluation of integrals

In this section we outline the computational details for the FE discretization of the DtN and PML based formulations in Section 5. Furthermore, we present results on integration of weakly singular kernels that are used in the new numerical sorting scheme in Subsection 7.6. For convenience of the reader we provide a summary of the used standard techniques.

7.1.1 Curved elements

Consider a physical element $K$ and let $\mathcal{K}:=(-1,1)^{d}$ denote the master element. Numerical quadrature is used to integrate a function over $K$ and when high order polynomial spaces are used, it is convenient to compute information from the shape functions $\varphi_{j},\,\nabla\varphi_{j}$ in the master element $\mathcal{K}$ and then store it. In this way we gain in performance as computations involving high order polynomials are expensive. Consequently, functions defined over a, possibly curved, physical element $K$ are mapped to act over $\mathcal{K}$, where numerical quadratures become simpler. Then, the mapping $X_{K}:{\mathbb{R}}^{d}\rightarrow{\mathbb{R}}^{d}$ transforms coordinates as $K=X_{K}(\mathcal{K})$. The action of the mapping is enforced by the Jacobian’s determinant $J=\det(DX_{K})$ [36, Sec. 3.3], [37, Sec. 3.4], where integrals transform as

[TABLE]

In the case where $K$ is a line, quadrilateral or a brick element, the explicit expression for $X_{K}$ is a known linear transformation. When $K$ has curved edges, then $X_{K}$ can be described by the so-called theory of Transfinite Interpolation [38], and the implementation and computational details can be found in [39], [37, Sec. 3.2]. A general rule of thumb is that the bending of the edges must be small compared to the diameter of the element, and that the angles at the element corners should be close to $\pi/2$. For further details and explicit error estimates on curved elements the reader is referred to [36, Sec. 3.3], [40, Sec. 6.7]. For the description on how $\varphi_{j},\,\nabla\varphi_{j}$ transform from $K$ to $\mathcal{K}$, and other related details, the reader is referred to [37, Sec. 3.3].

7.1.2 Numerical quadrature

In this subsection, we briefly revise numerical integration by Gauss-Legendre quadratures [41, Ch. 4]. In the one dimensional case, integration over the master element $\mathcal{K}$ is approximated by formulas of the form $\int_{\mathcal{K}}f(x)\,dx=\sum_{i=1}^{m}w_{i}f(x_{i})+E$, with $w_{i}$ the quadrature weights and $x_{i}$ the quadrature nodes. The coefficients $w_{i}$ are all positive [35, Sec. 8.4]. The so-called quadrature error or remainder is denoted $E$, and under the assumption that $f\in C^{2m}(\mathcal{K}),\,m\in\displaystyle\mathbb{N}$, then the $m$-point Gauss quadrature’s remainder satisfies

[TABLE]

Then, the Weierstrass approximation theorem [35, Sec. 1.2] guarantees the existence of a polynomial $P(x)$ such that $\sup_{x\in\mathcal{K}}|f(x)-P(x)|\leq\delta$, for an specified $\delta>0$. In this way $w_{i},x_{i}$ can be set to minimize $|E|$, and $P(x)$ is integrated exactly. Particularly, if $f$ is a polynomial of order $p<2m$, the remainder vanishes and the quadrature gives the exact value of the integral.

An effective way of reducing $|E|$ is by increasing the polynomial degree $p$ until the residual is below $\delta$. In the quadrature formula, increasing $p$ is equivalent to increasing the number of evaluation points $m$.

In ${\mathbb{R}}^{d}$, $d>1$, it is possible to derive similar quadrature formulas for the integration of functions that are in $C^{2m}(\mathcal{K})$. The resulting formulas are the so-called composite Gauss quadratures [41, Sec. 4.2.1], which can be constructed in our case by using tensor products of one dimensional quadrature rules. Particularly, when a physical element $K$ is allowed to be curved, we approximate integration with the formula

[TABLE]

where ${\rm w}_{j}=w_{j}J(x_{j})$ and ${\mathrm{x}}_{j}=X_{K}(x_{j})$ corresponding to a composite Gauss quadrature with $m^{d}$ weights and nodes $w_{j},x_{j}$. Finally, we discuss quadrature rules when integration is performed over a domain $\Omega\supseteq\cup_{i}K_{i}$ defined as the union of several elements $K_{i}$. The integrand is now required to be piecewise smooth $f\in C^{2m}(K_{i})$, for $i=1,2,\ldots,N_{\text{elements}}$. Then, we obtain

[TABLE]

where for each element $K_{i}$, we have the quadrature pairs $\mathrm{x}_{j}(K_{i}),\mathrm{w}_{j}(K_{i})$, similarly as in expression (38). The polynomial spaces that we use for $d=2,3$ are based on the tensor product of one dimensional finite element spaces [41, Sec. 2.2], [42]. As we work with piecewise smooth coefficients, we make sure that the jumps of the integrand coincide with the possibly curved element edges $\partial K_{i}$, such that $f\in C^{2m}(K_{i})$. This allow us to use quadrature rules in each individual element and guarantee convergence of the error of the numerical integration.

Further details on Gaussian quadratures can be revised in e.g. [35, Ch. 8], and implementation details are provided in [41, Ch. 4].

7.1.3 Integrating weakly singular kernels

In the discretization of the Lippmann-Schwinger formulation (26), we encounter the problematic situation where the kernel $\Phi(x_{i},y)$ is required to be evaluated at $y=x_{i}$. In the one-dimensional case ($d=1$), this obstacle is easily overcome by splitting the integration interval and performing two separate integrations by using Gauss-type of quadratures rules.

In higher dimensions ($d=2,3$) the kernel is weakly singular [43, Sec. 2.3], what makes the integration in (26) more demanding. This difficulty can be overcome by specializing the quadratures [44, 45, 46]. Usually, extra effort is spent in refining adaptively elements $K_{l}$ containing the singularity. Then, a Nyström type of high order quadratures, combined with interpolation in polar coordinates along with other techniques are used in order to keep $|E(\Omega_{r})|$ small to desired order. As expected, the challenge becomes more pronounced in higher dimension [45, 46]. These techniques where successfully tested [45, 46] for the solution of scattering problems for a given incoming wave employing the Lippmann-Schwinger formulation. However, our case is very different as we look for scattering resonances where the corresponding eigensolver is computationally more demanding than a linear solve.

7.2 Solution of the nonlinear eigenvalue problems

The approximation of resonances based on the DtN formulation for $d=2$ leads to the matrix problem in (29). The solution of this NEP is based on the solution strategy presented by Araújo et al. [19], where we use a specialization of the Infinite Arnoldi method [47, 48] called the tensor infinite Arnoldi method (TIAR). In particular we introduce a pole cancellation technique in order to increase the radius of convergence for computation of eigenvalues that lie close to the poles of the matrix-valued function.

For the approximation of resonances when the permittivity function is described by the rational model (7), we prefer to solve the corresponding matrix NEP by the techniques presented by Araújo et al. [1], which is a specialization of the solver by Güttel et al. [49] implemented in the SLEPc library [33].

7.3 Properties of volume integral equations for resonance computation

We discretize volume integral equations by using the scheme presented in (26), which does not require any boundary conditions as mentioned in Remark 1. In Section 6.2 numerical computations illustrate that the spectrum of the resulting discrete operator exhibits desired stability properties for perturbations of $a$.

Remark** 7****.**

The resulting system matrices has dimensions comparable to FE matrices: $N_{LS}\times N_{LS}$, with $N_{LS}\leq ch^{-d}$. However, the matrices are dense ($\mathcal{O}(N_{LS}^{2})$ storage), non-symmetric, and the elements of $T(\omega)$ are transcendental functions of $\omega$.

The consequences of Remark 7 in a NEP solution strategy is that the matrix $T(\omega_{j})$ from (26) must be re-assembled for each new iteration $\omega_{j+1}=\omega_{j}+\delta\omega$, which is computationally expensive, especially for problems of dimension $d>1$.

7.4 Memory requirements

Let $N_{K}$ be the number of cells in a triangulation in space dimension $d$, $p$ the polynomial degree of the basis functions in use, and $w=16$ bytes is the memory required to store a complex number in double precision. Given the number of non zeros elements $N_{z}$ in a matrix, the memory required to store it is $W=N_{z}\times w$.

In the collocation method given in Section 4.1, matrices are dense and we get $N^{D}_{z}\approx[N_{K}\times p^{d}]^{2}$. In turn, the FE matrices from Section 5.2 are sparse, and each cell in the triangulation contributes with a block of size $[(p+1)\times(p+1)]^{d}$ support points. Additionally, we have scattered connections of order $(p+1)^{d}$ with neighboring cells that we omit for simplicity. A simple estimation gives $N_{z}^{S}\approx N_{K}\times(p+1)^{2d}$. Furthermore, by assuming $m$ the division in a one dimensional partition, it is reasonable to have $N_{K}=m^{d}$, and $m=10,10^{2},10^{3}$ for a small, moderate and large problem respectively.

We list our simple estimations in table 1 for $p=2$, for both FE discretization methods in Section 5.2 and collocation methods for volume integral equations in Section 4.1. As expected, dealing with FE matrices is a standard way of discretizing wave problems and results in manageable sparse matrices for current computer memory constraints. It is evident that the load becomes larger with higher space dimension and number of cells, from where matrices for $d=3$ can be stored only for small and moderate problems. In the case of matrices from the collocation strategy in (26), we conclude that storage becomes computationally unfeasible in higher space dimensions, for moderate and large problems. Additionally, working with higher polynomial degree rules out even small problems.

7.5 Computational platform and details

All numerical experiments have been carried out using the finite element library deal.II [42] with Gauss-Lobatto shape functions [41, Sec. 1.2.3]. For fast assembly and computations with complex numbers the package PETSc [50] is used.

The computational platform was provided by the High-Performance Computing Center North (HPC2N) at Umeå  University, and all experiments were run on the distributed memory system Abisko. The jobs were run in serial on an exclusive node: during the process, no other jobs were running on the same node. Node specifications: four AMD Opteron 6238 processors with a total of 48 cores per node.

7.6 Computational details of the sorting scheme

In order to evaluate the sorting scheme, we are interested in computing the integrals from (26) as accurately as possible. The available FE machinery for computing integrals over $\Omega_{0}$ facilitates the numerical integration, which is done similarly as described in Section 7.1.2.

Due to the growth of most resonant modes, the point-wise residual $|(T^{\gamma}(\omega^{\gamma})u^{\gamma})(x_{j})|$ is expected to be larger for $x_{j}\in\Omega_{0}$. Additionally, as discussed in Section 7.1.3, due to the unboundedness of $\Phi(x_{j},x_{j})$, the computation of (26), and (32) for $x_{j}\in\Omega_{r}$ requires considerable more effort compared to its evaluation for $x_{j}\in\Omega_{0}$. The apparent reason for this is that for $x_{j}\in\Omega_{r}$ and $d>1$, we have to numerically compute an integral with a weak singularity as discussed in Section 7.1.3.

Below, we write explicitly the steps involved in computing $\delta^{\gamma}(\Omega_{a})$ from Definition 5. First, we split the integration into separate parts over $\Omega_{0},\Omega_{r}$ and use the composite quadrature rules (39) for evaluating the integrals. We need the following definitions.

Definition 8**.**

Let $\mathcal{K}_{0}:=\{i:K_{i}\subset\Omega_{0}\}$ and $\mathcal{K}_{r}:=\{i:K_{i}\subset\Omega_{r}\}$ be index sets defined over $\Omega_{0}$ and $\Omega_{r}$, respectively. We define the sets $\mathcal{X}_{0}:=\cup_{i\in\mathcal{K}_{0}}\{\mathrm{x}_{j}\in K_{i}\}_{j=1}^{m^{d}}$, $\mathcal{X}_{r}:=\cup_{i\in\mathcal{K}_{r}}\{\mathrm{x}_{j}\in K_{i}\}_{j=1}^{m^{d}}$, and denote by $\mathcal{I}_{0}$, $\mathcal{I}_{r}$ the resulting extracted index sets from the new ordering. Additionally, let $\mathcal{X}_{a}:=\mathcal{X}_{0}\cup\mathcal{X}_{r}$ and $\mathcal{I}_{a}$ its corresponding ordering. In this way we obtain the new quadrature rules

[TABLE]

Then, we have

[TABLE]

where $\alpha_{k}:=|(T^{\gamma}(\omega^{\gamma})u^{\gamma})(\mathrm{x}_{k})|,\,\,\beta_{l}:=|(T^{\gamma}(\omega^{\gamma})u^{\gamma})(\mathrm{y}_{l})|$, for $\mathrm{x}_{k}\in\Omega_{r},\,\mathrm{y}_{l}\in\Omega_{0}$. From (26), the evaluation of $T(\omega)u$ involves an integration over $\Omega_{r}$, which we refer to as inner loop. Then, for each $k$ in $\alpha_{k}$, and each $l$ in $\beta_{l}$ we compute an inner loop additional to the explicit integration shown in (41).

Definition 9**.**

For a given quadrature, we define the discrete pseudospectrum indicator as

[TABLE]

7.6.1 Computational costs

In this subsection, we estimate the computational cost for performing the operations involved in (42). Computationally, the errors in (41) require special treatment as discussed in 7.1.3. However, for simplicity of the estimations, we disregard additional costs from integration of weakly singular kernels in higher dimensions.

We estimate the costs in terms of the evaluation of $u^{\gamma}(x_{j})$ and $\Phi(x_{i},x_{j})$ in complex double precision, which combined account for the heaviest work load in each individual term of (42). The evaluation of these two operations account for a rough computational time $t_{q}\approx 10^{-6}s$ bench-marked on the processor Intel Core i7-3770, CPU: 3.40GHz.

For the estimation, assume that we have $N_{K}$ cells $K_{i}\subset\Omega_{a}$ and choose a finite element space of degree $p$. Then, we have $N_{K}\times(p+1)^{d}$ terms in the outer loop and the inner loop requires $N_{r}\times(p+1)^{d}$ terms, where $N_{r}$ denotes the number of cells in $\Omega_{r}$. The estimate for the computational cost for $\|T^{\gamma}(\omega^{\gamma})u^{\gamma}\|^{2}_{\Omega_{a}}$ is about $N_{K}\times(p+1)^{d}\times N_{r}\times(p+1)^{d}=N_{K}\times N_{r}\times(p+1)^{2d}$ evaluations of the kernel. Assume $m$ as the size of a one dimensional partition, from where it is reasonable to have $N_{r}=m^{d}$, $N_{K}=(cm)^{d}$, for $c>1$. Then, the cost is given by $c^{d}\times(m\times(p+1))^{2d}$.

Aiming at getting a better intuition of the requirements of the computation, we can check estimations for the cost, inner time $t_{i}$, total time $t$ running on a single processor. We set $c=2$ and for a small, moderate and large problem we use $m=10,10^{2},10^{3}$ respectively. In Table 2 we show the estimations for the computational costs and times required by (41).

The presented sorting scheme is fully parallelizable and the total time of execution can be reduced by a factor of ten by using additional cores. However, the conclusion of Table 2 is that the computational cost is extremely high for realistic computations. Basically, evaluating $\delta^{\gamma}_{h}(\Omega_{a})$ as in (42) results in sorting schemes that are far more expensive than the solution of the NEP (29).

7.6.2 Sorting strategy

A pseudospectrum strategy based on Definition 9, consist of the sorting of computed pairs $(\omega_{m}^{\gamma},\xi_{m})$, solution to $\eqref{eq:matrix_nep}$, according to their respective indicator $\delta^{\gamma}_{h}(\Omega_{a})$. As discussed in Section 7.1.3, the evaluation of the $\beta_{l}$ in (41) requires special treatment such as non-standard quadrature rules similar to the ones introduced in [44, 45, 46]. Additionally, as estimated in Section 7.6.1 and Table 2, the evaluation of $\delta^{\gamma}_{m}(\Omega_{a}):=\|u_{m}-Ku_{m}\|_{L^{2}(\Omega_{a})}$ in higher dimensions is prohibitively expensive since $Ku_{m}$ is a volume integral operator. With these issues in mind, our aim is to propose an approximated version of $\delta^{\gamma}_{h}(\Omega_{a})$ from Definition 9 such that we improve in performance compared to the results in Table 2, and we avoid the use of specialized quadrature schemes for the evaluation of singular kernels. Our goal is to reduce the complexity of the sorting scheme, such that the cost of the new strategy scales linearly with the cost of the computation of the inner loop.

In the remaining of the section, we present an alternative sorting alternative based on Definition 9 with computational cost that scales with the cost of evaluating the inner loop.

7.6.3 Sorting estimations

The orderings given in Definition 8 are used to group quadrature pairs over $\Omega_{0}$, $\Omega_{r}$ and $\Omega_{a}$ separately. The resulting pairs are written as $\{\mathrm{x}_{k},\mathrm{w}_{k}\}$ for index $k\in\mathcal{I}_{q}$ with $q=0$ corresponding to quadrature pairs from $\Omega_{0}$, $q=r$ from $\Omega_{r}$, and similarly for $q=a$ from $\Omega_{a}$. Then, from Definition 9 we have

[TABLE]

In the estimate (43), we used the properties that the quadrature weights ${\rm w}_{j}$ are positive, and that $|\Omega^{h}_{q}|=\sum_{j\in\mathcal{I}_{q}}\mathrm{w}_{j}$.

We base our sorting strategy in the following definition. A function $u\in C(\Omega_{a})$ satisfies

[TABLE]

Since both the FE approximations and the Nyström approximations converge as we increase the dimension of the discretization, then we expect that the estimate (44) wraps around $\delta_{h}^{\gamma}(\Omega_{a})$ tighter for finer discretizations.

Definition 10**.**

Sorting indicator:* For a given eigenpair $(\omega_{m}^{\gamma},v_{m}^{\gamma})$ we define*

[TABLE]

Then, (44) suggests an alternative strategy as an approximation of the pseudospectrum indicator in definition 5. Basically, the goal is to test points hoping that when the number of evaluations is large we have a good approximation to the upper bound (45). In order to make the most of each evaluation, we consider the following heuristic ordering of the evaluation points:

We assume that $|T^{\gamma}(\omega^{\gamma})u^{\gamma}(x)|$ is large at points $x$ where $|u^{\gamma}(x)|$ is large.

• 2.

Due to the rapid growth of eigenfunctions, we assume that non-convergent/spurious pairs exhibit $\max_{y\in\Omega_{r}}|u(y)-K^{\gamma}u(y)|<\max_{x\in\Omega_{0}}|u(x)-K^{\gamma}u(x)|$.

• 3.

We want to avoid large errors in the quadrature rules used. Then, it is convenient to evaluate residuals away from the singularities of the kernel $\Phi(x,y)$. These are located at points $x=y$.

The last reasons suggest that for identification of spurious pairs, it is convenient to compute first the residuals at points outside the resonator. This is, for the computation of (45), we should start by evaluating points in $\Omega_{0}$.

Definition 11**.**

Sampling set:* For given constants $l_{m},\,c>0$ we define the sampling set*

[TABLE]

Remark** 12****.**

The proposed strategy consists in evaluating the residual $|(T^{\gamma}(\omega^{\gamma})u^{\gamma})(x)|$ in points $x\in X_{m}$ as suggested by (45), where the evaluation points are not clustered together. Additionally, we want to exclude points ${x}_{l}$, such that $|u^{\gamma}({x}_{l})|\approx 0$, and we want to avoid integrating over cells with a singular kernel. For this, we select points ${x}_{j}\not\in\Omega_{r}$ such that $\inf_{y\in\Omega_{r}}|x_{j}-y|\geq l_{m}$, in order to avoid the singularity peak. Consecutively, we select $l_{m}$ such that $|\Phi(\omega_{m}l_{m})|\leq C_{1}$, which guarantees that the integrand is bounded and quadrature rules remain accurate. The benchmarks in ${\mathbb{R}}^{d}$, $d=1,2,3$ presented in this paper indicate that for practical computations using $N_{s}\approx 10$ evaluations is a good choice. Ultimately, we use the normalization $v_{m}:=u_{m}/\|u_{m}\|_{\Omega}$, with $\Omega:=\Omega_{a}\cup\Omega_{P\!M\!L}$ for the PML formulations. Then, we filter out added PML eigenvalues described in Definition 6 by requiring the condition $\|v_{m}(x)\|^{2}_{\Omega_{a}}/|\Omega_{a}|>\|v_{m}\|^{2}_{\Omega_{P\!M\!L}}/|\Omega_{P\!M\!L}|$.

The result of applying Remark 12 is an effective and inexpensive way of testing the computed pair $(\omega^{\gamma},\xi)$, where the cost scales linearly with the inner loop. If the resulting sorting indicator is larger than a user pre-defined threshold, then we disregard the eigenpair and continue testing the next eigenpair in the collection.

8 Applications to metal-dielectric nanostructures

In this section we study four interesting metal-dielectric configurations, from where numerical approximations to resonances and resonant modes are computed. Consecutively, eigenpairs are tested and solutions are sorted according to their corresponding pseudospectrum indicator (45). The sorting strategy is tested on problems where exact pairs are known. Additionally, we consider as reference problem a test case used in [51].

The first three configurations serve as benchmarking strategies for non-dispersive and piecewise constant material properties. Then, we apply the sorting algorithm on a configuration introduced in [1] where a metal coating is motivated from realistic applications in nano-photonics. Here, three different relative permittivity models are used: $\epsilon_{v}:=1$ (Vacuum), $\epsilon_{s}:=2$ (Silica), and $\epsilon_{metal}$ (Gold), modeled by a sum of Drude-Lorentz terms (7). For $\epsilon_{metal}$ we use the data given in table 3 gathered in [52]. This model of Gold has been extensively tested and has validity for $\omega\in[0.5,6.5]\,eV$, where $eV$ denotes electron volt.

We introduce a demanding configuration where the refractive index is a continuous function of space motivated from the so-called graded materials. Finally, we consider an acoustic benchmark problem in $\displaystyle\mathbb{R}^{3}$.

8.1 Modeling details

In finite precision arithmetic we prefer to work with dimensionless quantities, where we transform from dimensionless variables to physical variables (denoted with ~). We use common physical constants in SI units: $\hbar$ is the scaled Planck’s constant, $c$ is the speed of light in vacuum, and $e$ is the electron charge. In the numerical computations, we use the scaling factors $W=eV/\hbar$ in Hertz and $L=2\pi c/W$ in meters. Then, we define the dimensionless quantities

[TABLE]

The resulting length factor is $L=1239.842\,nm$, from where our spectral window becomes numerically equivalent to $eV$ scaling.

8.2 Benchmarks in 2D

The next two problems have radial symmetry centered at the origin, and the solutions expressed in polar coordinates $(r,\theta)$, will be written in terms of Bessel and Hankel functions of integer order $m$. In this simple case outgoing solutions of (9) satisfy

[TABLE]

where  supp$\,(n-1)\subset B(0,a)$. In subsections 8.2.1 and 8.2.3, we present solutions satisfying (9) and (48) for specific permittivity profiles.

8.2.1 Standard benchmark: Single disk problem (SD)

Denote by $u=u_{1},\,n=n_{1}$ the restrictions of $u,n$ to $\Omega_{1}:=B(0,R)$, and set $n=n_{2}=1$ elsewhere. The corresponding exact eigenfunctions to (9) and (48) read:

[TABLE]

The eigenvalues $\omega$ corresponding to $m=0$ are simple and those corresponding to $m>0$ are semi-simple and have multiplicity $\alpha=2$. The exact eigenvalue relationship for TM and TE can be written as

[TABLE]

where $g=n_{1},\,g=1/n_{1}$ corresponds to the TM polarization and TE polarization, respectively. In the case of a semi-simple eigenvalue $\omega_{j}$, the exact solution consist of two eigenfunctions $u_{j,1}$ and $u_{j,2}$. Then, we define the linear combination $v_{j}=c_{1}u_{j,1}+c_{2}u_{j,2}$ for the constants $c_{1},c_{2}\in\displaystyle\mathbb{C}$. In this way, after we have computed an eigenpair $(\omega^{\gamma},u^{\gamma})$, we compare $u^{\gamma}$ with the closest exact solution $v_{j}$. That is, we find $c_{1},c_{2}$ by solving the minimization problem

[TABLE]

This can be done in the discrete setting with a standard least-squares optimization routine.

For numerical computations we use $R=1$, and $a=2,3$. Additionally, we place the disk center a distance $s=0.2$ from the origin such that many terms are needed in the DtN for approximation of resonances.

Comparison with respect to exact solutions and LS-projections: We denote with $u^{\gamma}$ the computed FE eigenfunction, with $u$ the best fitted exact solution and with $K^{\gamma}u^{\gamma}$ the Lippmann-Schwinger projected eigenfunction. The best fitted exact eigenfunction $u$ is obtained by solving the minimization problem (51).

For convenient visualization and comparison of the eigenfunctions we define the straight line $\mathcal{L}\subset\Omega_{0}$ as

[TABLE]

where $s=0.2$. The vertices defining $\mathcal{L}$ are located a small distance away from the resonator and the boundary $\Gamma_{a}$ by setting $\epsilon=10^{-3}$. This line has been chosen in order to avoid computations over an edge of an element such that we avoid super-convergence of residuals.

In Figure 5, we gather representative results for the single disk problem. The depicted solutions in the figure correspond to the TM polarization with $p=2$, ref$\,=1$, and $a=3$. Vertically we show different representative eigenfunctions that are drawn for $x\in\mathcal{L}$. In the left panel of the figures we present real parts of $u$ (black), $u^{\gamma}$ (red), $Ku^{\gamma}$ (blue), and imaginary parts in the middle panel. In the right most panel we present the residuals $|u(x)-u^{\gamma}(x)|$ in red dashed and $|T^{\gamma}u^{\gamma}(x)|=|u^{\gamma}(x)-K^{\gamma}u^{\gamma}(x)|$ in blue dotted lines.

Results: In Figure 5 we observe that for eigenpairs that are good approximations then $|T^{\gamma}u^{\gamma}(x)|$ follows very closely $|u(x)-u^{\gamma}(x)|$, showing an important resemblance. This presents a measure convergence of FE computations. Additionally, some $u^{\gamma}$ are seen to grow exponentially away from the resonator. However, the corresponding $K^{\gamma}u^{\gamma}$ grows at a slower rate and seems to agree more with $u$ than $u^{\gamma}$ itself. Finally, in most cases, $|T^{\gamma}u^{\gamma}(x)|$ is smoother and more flat than $|u(x)-u^{\gamma}(x)|$.

The application of the sorting scheme described in Remark 12 to this problem gives the results presented in Figures 6 and 7. Exact eigenvalues are marked with dots and computed eigenvalues are marked with colored circles, where the color is given by approximating the sorting indicator (45) as described in Remark 12. In Figure 6 we present plots for the TM polarization and in Figure 7 results for the TE polarization. The left panels are computed by placing the DtN at $a=2$, while $a=3$ in the panels on the right, and both discretizations have the same FE up to $a=2$.

From the figures we observe that increasing $a$ results in an increase in the number of eigenvalues in the given spectral window, which is expected from the discussion in Section 6.2. Moreover, the added eigenvalues from $a=3$ pollute larger regions in the spectral window compared to $a=2$. This conclusion can also be obtain by noticing that the minimum in $\tilde{\delta}$ increases for larger $a$, which is related to the fact that exponential growth of eigenfunctions becomes a challenge for FE discretizations.

Additionally, we see that computed eigenvalues with small sorting indicator resemble the pattern drawn by the exact eigenvalues. As the value of the indicator increases the corresponding eigenvalues loose this property and become erratic. In this way we can choose a threshold and for example keep eigenvalues $\omega^{\gamma}_{m}$ with $\tilde{\delta}_{m}\leq 10^{-2}$ to obtain a set of eigenvalues resembling the behavior of the exact eigenvalues. Moreover, in Figure 6 we illustrate results for computations on a sequence of finer meshes. We observe that computation on finer meshes improve the accuracy of computed eigenvalues and exhibit less spurious eigenvalues than computation on coarser meshes. Additionally, we obtain lower values for $\tilde{\delta}_{m}$ for finer discretizations compared with coarser discretizations.

Finally, it should be noted that the indicator $\tilde{\delta}$ resembles the FE discretization error. Namely, $\tilde{\delta}_{m}$ increases with $\operatorname{Re}\omega_{m}$, which is the expected behavior for the FE discretization error of non-dispersive Helmholtz problems [15, 53, 1]. Additionally, eigenvalues with small indicator values appear very close to the exact eigenvalues of the problem.

8.2.2 Reference configuration: Single square problem (SS)

For this problem, we set $\Omega_{1}:=[-s,s]^{2}$, $n=n_{1}$ for $x\in\Omega_{1}$ and $n=1$ elsewhere. The correct parameters for obtaining the results presented in [51, Fig. 10] are $s=1.5$ and $n_{1}=0.2$ for TM polarization.

Results: For this problem, we compute all eigenvalues in the spectral window only once, and cheaply test each computed eigenpair. The results are plotted in Figure 8, by using a DtN in the left panel, and PML in the right panel, where both formulations use the same FE basis in $\Omega_{a}$. Again we can see the effectiveness of using the indicator $\tilde{\delta}$. Correct approximations to resonances are easily identified from the overwhelming rest of computed eigenvalues. Additionally, we see that the FE-DtN produces better approximations to resonances than the FE-PML. For example, we mark with dots the two reference values given in [51, Fig. 9], and we see that for the presented discretizations in Figure 8 the reference eigenvalues have a lower indicator $\tilde{\delta}_{m}$ for the FE-DtN than with the FE-PML.

Comparison from the results in Figure 8 against those in [51, Figs. 10, 11, 12] not only illustrates the reliability of the identification of true approximation to resonances by using the scheme in Definition 5, but it also shows its simplicity, flexibility and large coverage in the complex plane.

8.2.3 Benchmark with dispersion: Single coated disk problem (SCD)

In this configuration, we consider a resonator consisting of a dielectric disk with a uniform coating layer. The geometry is described by two concentric circumferences of radii $0<R_{1}<R_{2}$, with vacuum as surrounding medium. The inner disk has constant relative permittivity index, and is coated by a layer of gold. We set $n_{1}=\sqrt{\epsilon_{s}},$ and $n_{2}:=\sqrt{\epsilon_{metal}}$ is the value such that $\operatorname{Im}\{n_{2}\}$ (absorption coefficient) is positive.

The exact solutions satisfy (9) and (48) with $R\geq R_{2}$. The resonance relationship reads

[TABLE]

where for TM, $g:=(n_{1},n_{2},n_{2},n_{1},1,n_{2},n_{2},1)$, and for TE, $g:=(n_{2},n_{1},n_{1},n_{2},n_{2},1,1,n_{2})$. The parameters used for the computation are $R_{1}=0.8$, $R_{2}=1.0$ with scaling factor $L=1239.842\,nm$.

A complex Newton root finder [54] is then used to compute very accurate approximations of the resonances. For each $m$ in equation (53), we search numerically the resonances $\omega_{m,1},\omega_{m,2},\ldots$ with machine precision stopping criterion. In [1, Table 2], we list a selection of resonances computed from (53).

Results: The results for both polarizations are gathered in Figure 53. Here, we observe typical behavior of dispersive resonators $i$) there exist clustering of eigenvalues close to the poles and zeros of the Drude-Lorentz model (7), $ii$) in the TE polarization we have clustering of eigenvalues to the so-called plasmonic branch points of the model, which are those values $\omega$ such that $\epsilon_{\text{m}etal}(\omega)=-\epsilon_{j}$, with $j=0,1$. In other words, $-\epsilon_{\text{m}etal}(\omega)$ matches the value of a neighboring dielectric constant (juncture). In the figure, the bottom panels are close up windows showing such accumulations and the corresponding values for the resulting eigenpair indicator $\tilde{\delta}_{m}$ as described in Remark 12. As expected, the value of $\tilde{\delta}_{m}$ increases when approaching a critical value (pole, or zero), this behavior is expected since close to a critical point the value $|n(\omega)\omega|$ increases and the resulting eigenfunction oscillates more rapidly, and the FE error increases. This is further discussed in [1].

8.2.4 Configuration with continuous $n(x)$: Bump problem

Consider the following refractive index:

[TABLE]

subject to the compatibility conditions: $P_{3}(0)=1,\,P_{3}(R)=0,\,P^{\prime}_{3}(0)=0,\,P^{\prime}_{3}(R)=0$, and $R=1$.

8.3 Acoustic benchmark in 3D: Single ball problem (SB)

This case is analogous to the Single Disk problem 8.2.1 with $d=3$.

Denote by $u=u_{1},\,n=n_{1}$ the restrictions of $u,n$ to $\Omega_{1}:=B(0,R)$, and set $n=n_{2}=1$ elsewhere. The corresponding exact outgoing resonant modes of (9) read:

[TABLE]

The eigenvalues $\omega$ corresponding to $m=0$ are simple and those corresponding to $m>0$ are semi-simple and have multiplicity $\alpha=2m+1$. The exact eigenvalue relationship can be written as

[TABLE]

*Results for problems 8.2.4, and 8.3: * Similarly as pointed out in former discussions of this section, the results gathered in Figures 10 and 11, confirm the effectiveness from the application of the sorting scheme described in Remark 12. The results in the figures illustrate a clear method to effectively sort the computed eigenpairs by using an approximation of the sorting indicator (45).

9 Conclusions

We have presented a sorting scheme, based on the Lippmann–Schwinger equation, for marking of potentially spurious scattering resonant pairs in $\displaystyle\mathbb{R}^{d}$ Helmholtz problems. For all computations and for TM and TE polarizations, numerical experiments on a broad range of problems illustrate that the sorting scheme provides valuable information on the location of potentially spurious resonances at low computational cost. This information can, for example, be used to determine the most important modes in quasimodal expansions.

Acknowledgments

This work is founded by the Swedish Research Council under Grant No. 621-2012-3863.