## 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 result in approximations of the true resonances but in practice also a large number of solutions that are unrelated to the true resonances. Those *spurious eigenvalues* are a major problem in resonance computations. The origin of unphysical eigenvalues in scattering resonance computations is *spectral instability*, which is common for non-normal operators trefethen07 ; davies97 .
Spectral instability is known to be less problematic with volume integral equations compared with formulations based on differential operators. Therefore, we proposed in araujo+engstrom+2017
to use a Lippmann-Schwinger type of integral equation for removing spurious solutions in a one-dimensional setting.

In this paper, we show that it is possible to extended the approach in araujo+engstrom+2017 to higher dimensions. In particular, we show that the used sorting scheme can be computed cheaply and give valuable information on the location of spurious eigenvalues. This test is motivated by spectral stability properties of the Lippmann-Schwinger formulation of the problem. The idea is then to test each computed pair and obtain an pseudospectral indicator . Then we sort all computed pairs from the smallest to the largest value of . Finally, we can choose a user defined tolerance , and drop all pairs with . In araujo+engstrom+2017 we presented computations on several test cases for space dimension . The strategy was successful, since the proposed filtering scheme was able to identify physical pairs. However, the application of the ideas presented in araujo+engstrom+2017 for turn out to be more challenging. The effective application of these ideas in higher dimensions requires further considerations that are addressed in the current work. The new results confirm the efficiency of the proposed sorting scheme in , .

## 2 Electromagnetic and acoustic scattering resonances

Assume that is independent of and consider electromagnetic waves propagating in the -plane. The -independent electromagnetic field is then decomposed into transverse electric (TE) polarized waves and transverse magnetic (TM) polarized waves Cessenat1996 . This decomposition reduces Maxwell’s equations to one scalar equation for and one scalar equation for . The TM-polarized waves and the TE-polarized waves satisfy formally

(1) |

respectively. For the scattering resonance problems, and are assumed to be locally integrable functions that satisfy an outgoing condition MR1350074 ; MR1037774 .

Let the physical domain be an open ball of radius with boundary ,
and let be the bounded domain defining the resonators. Hence, we assume that the relative permittivity in is one. Furthermore, let denote the union of disjoint resonators , satisfying for all , as shown in figure 1. A *scattering resonance* was in gopal08 defined as a complex number for which the Lippmann-Schwinger equation

(2) |

has a non-zero solution . The integral operator in (2) is for TM/TE waves given by

(3) |

Here, is known as the outgoing Green function in free space for the Helmholtz equation MR1350074 ; MR1037774 . Notice that while we are interested in , the integration in (2) is only performed over , since the integration over the air region 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 in using a standard laptop; see e.g. bra13 ; araujo+engstrom+2017 . However, accurate computations of eigenpairs of (2) in higher dimensions would require huge computer resources; See MR2406839 and the discussion in Section 7.6.

### Acoustic scattering resonances in

Sound-soft materials are characterized by the speed of sound and we assume that acoustic resonators are defined by . Then, the acoustic pressure satisfies formally the Helmholtz equation

(4) |

where the outgoing condition can be expressed as satisfying an expansion in spherical harmonics outside the open ball ihl98 . Moreover, is a scattering resonance pair if (2) holds with

(5) |

Note that the acoustic problem in , is analogous to TM-polarized electromagnetic waves with .

### 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 Wang2016 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 geometry MR3622414 . The most popular method to compute resonances in , 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) MR2485441 and DtN maps araujo+engstrom+jarlebring+2017 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 jar14 ; jar15J . Discretization with FE of the PML and HIP 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 is -independent and in the general case the non-linearity in is completely determined by . 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 denote an unbounded closed linear operator in a Hilbert space with domain , spectrum , and resolvent set . Then exhibits high *spectral instability* if for a very small there exist many and such that

(6) |

even though is not close to MR1912874 . This is closely related to the pseudospectrum which is defined as the union of and all in the resolvent set for which it exists an such that (6) holds. The generalization of those results to an operator function is straightforward and we will in some of the numerical computations rely on the following alternative characterization of the pseudospectrum:

It is well known that PML and HIF based methods encounter high spectral instability MR2485441 ; araujo+engstrom+2017 . Methods based on a DtN map encounter medium spectral instability araujo+engstrom+jarlebring+2017 ; araujo+engstrom+2017 and integral equation based methods encounter low spectral instability araujo+engstrom+2017 . 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 electromagnetics, the material properties of non-magnetic metals are characterized by the complex relative permittivity function , which changes rapidly at optical frequencies . The most common accurate material model is then the Drude-Lorentz model

(7) |

where and , , are non-negative Cessenat1996 . Hence, the Maxwell eigenvalue problem in the spectral parameter is nonlinear for metal-dielectric nanostructures. Assume that the domain of the resonators can be written in the form and let

denote the characteristic function of the subset

. For material properties that are piecewise constant in , we assume a permittivity function in the form(8) |

where the dependencies on in for are of Drude-Lorentz type (7). In addition, we will consider graded material properties, meaning that is a continuous function in . In linear acoustics, the speed of sound 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 . In the following, we use the notation

(9) |

where for the TM-case and for the TE-case.

We define for the forms

(10) |

where in (10), and are functions of . Let denote the set of values that are zeros or poles of and set .

### 3.1 DtN based methods

Scattering resonances and quasi-normal modes restricted to can be determined from a problem with a Dirichlet-to-Neumann (DtN) map lenoir92 ; schenk2011optimization ; araujo+engstrom+jarlebring+2017 . Below we present variational formulations for , . Formally, is a scattering resonance pair if (9) holds in and

(11) |

where is the normal derivative.

#### 3.1.1 DtN formulation in

In one space dimension the scattering resonance problem restricted to is formally: Find a non-zero and a complex such that

(12) |

where the DtN-map at is

(13) |

Define for and the forms as in (10), and

(14) |

#### 3.1.2 DtN formulation in

In this subsection, we present a DtN formulation in polar coordinates . Let denote the Hankel function of first kind, then the DtN operator (11) on the circle has the explicit form

(16) |

and is bounded lenoir92 .

The resonance problem restricted to is formally to find non-trivial solutions such that (9) and (11) with (16) holds. The theory presented in lenoir92 can with minor changes be used in the present case to derive properties of a variational formulation of the problem.

Variational formulation: Let denote the union of the set of zeros of , and denote the operator (11) truncated after . The eigenvalues of the truncated version of (9)-(11) are determined by the following variational problem: Find and such that for all

(17) |

where the forms are defined as in (10), and

(18) |

### 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 kim09 . The method consists on attaching to a buffer layer of thickness , where outgoing solutions decay rapidly. The buffer domain is referred to as and the full computational domain is enlarged as shown in the Figure 1.

#### 3.2.1 PML formulation in

Let , and , with . The action of the PML is defined through the stretch function

(19) |

with , and where the polynomial is required to be increasing in , and is sufficiently smooth . For this we introduce the fifth order polynomial satisfying: and .

The PML problem is restricted to and the PML strength function has then the profile shown in Fig. 1. In the following sections, we consider the complex change of variable and transformation rule

(20) |

where in and in .

#### 3.2.2 PML formulation for

Approximation of resonances using a radial PML was analyzed in kim09 and we will here only consider the PML problem truncated to the ball . The used complex stretching functions are as in one-dimension of the form (19). Let denote the set of values that are zeros or poles of and let denote a partition into the PML-region and the part of the domain containing the resonators.

In the sequel we need the following definitions

(23) |

with the properties and .

It is clear that the PML coefficients are designed such that for , we obtain and , which is the same as . Hence, the PML operator restricted to corresponds to the original operator in problem (9) (no PML effect).

Variational formulation: The eigenvalues of (9) are then determined by the following variational problem

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

(25) |

Even though , and are defined in the whole , their action takes place only in .

## 4 Discretization of the Lippmann-Schwinger equation

In this section we present a collocation method for 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.6.

### 4.1 A Galerkin-Nyström method

We present a Galerkin-Nyström discretization method for linear Fredholm integral equations of the second kind, with kernels satisfying for in the compact set . In ikebe72 , this method is referred to as case (A) of the Galerkin methods, and convergence for the problem with sources is also discussed. In (gopal08, , Sect. 3.2) the method was used for resonance computations for .

Let be piecewise polynomial functions with the property , . We introduce the representation , and with the use of (2) we obtain the nonlinear eigenvalue problem: Find and such that

(26) |

The Nyström method consists in choosing the collocation points as the nodes of a high order quadrature rule. By doing this, the convergence of the scheme is considerably improved. In Section 7.6 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 Sakurai09 ; Beyn12 ; MR3423605 .

###### Remark 1.

The formulation in (2) uses information of the exact solution of the problem at every discretization node , through its fundamental solution . From where the numerical scheme is flexible in the way that it can be posed in the smallest domain , as well as in larger domains , 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 be covered with a regular and quasi uniform finite element mesh consisting of elements . The mesh is designed such that the permittivity function is continuous in each . Let be the length of the largest diagonal of the non-curved primitive and denote by the maximum mesh size .

Let denote the space of polynomials on of degree and define the dimensional finite element space

(27) |

The computations of discrete resonance pairs are for performed in the approximated domain using curvilinear elements Babuska92 . The meshes used are *shape regular* in the sense of (Schwab1998, , 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 . Let be a basis of . Then , and the entries in the finite element matrices are of the form

(28) |

The matrix eigenvalue problem is then: Find the eigenpars such that

(29) |

where the corresponding matrix valued function is

(30) |

In the case where is given as piecewise smooth function of space, we write (29) as

(31) |

with matrices

###### Remark 2.

Truncation of the DtN: Let be the smallest integer greater than or equal to . We use the rule according to the results in araujo+engstrom+jarlebring+2017 , where from the considered spectral window, 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 kim09 ; araujo+engstrom+2017 , which accounts for large enough and such that the search region is feasible. Additionally, we use the space for computations with the PML formulation.

Finally, we mention that all discretization methods presented in this work use the same FE space over .

###### Remark 4.

All formulations (LS, DtN, PML) use the FE triangulation , which is the restriction of to .

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 resonances

In this section we derive a discrete form of (2) that allow us to identify resonances from spurious solutions once we have computed FE solutions to (29) or to (31). The resulting expression for the sorting scheme is a discrete form of the condition , where is a FE solution restricted to . Let be a basis for and let be the -projection on . Then, the discrete Lippmann-Schwinger equation (26) is written in the form

###### Definition 5.

Pseudospectrum indicator: The computed eigenvalue belongs, for given , to the -*psudospectrum* if the pair satisfies . Then, for a given domain , we define the pseudospectrum indicator as

(32) |

Particularly, we want to be able to measure whether the computed eigenpair is converging to a physical pair. Naturally, non convergent pairs exhibit large values.

Additionally, we have that certain spurious eigenpairs introduced by using the PML method are easily identified and can be removed by using 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 araujo+engstrom+2017

, and eigenfunctions

exhibit oscillations in , but decay in the physical region . Then, by using the normalization , the PML critical eigenvalues exhibitwhich can be succesfully used as a filtering criterion for removing PML added eigenpairs.

## 7 Computational details

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.7. For convenience of the reader we provide a summary of the used standard techniques.

### 7.1 Master element and transfinite interpolation

Consider a physical element and let denote the master element. Numerical quadrature is used to integrate a function over and when high order polynomial spaces are used, it is convenient to compute information from the shape functions in the master element and then store it. In this way we gain in performance as computations from high polynomials are expensive. Consequently, functions defined over a physical element are mapped to functions over , where we perform integration. Then, the mapping transforms coordinates as . The action of the mapping is enforced by the Jacobian’s determinant (Strang2008, , Sec. 3.3), (pavel-solin04, , Sec. 3.4). Then, we have

(33) |

For the case where is a line, quadrilateral or a brick element, the explicit expression for is a known bilinear transformation. When has curved edges,

can be described by the so-called theory of Transfinite Interpolation

Gordon1973 , and the implementation and computational details can be found in Gordon73 , (pavel-solin04, , 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. For further details and explicit error estimates on curved elements the reader is referred to

(Strang2008, , Sec. 3.3), (Oden76, , Sec. 6.7). For the description on how transform from to , and other related details, the reader is refereed to (pavel-solin04, , Sec. 3.3).### 7.2 Evaluation of integrals

In this subsection, we revise briefly numerical integration by *Gaussian-Legendre* quadratures. In the one dimensional case, integration over the master element is approximated by formulas of the form , where are the quadrature weights, the quadrature nodes, and is the quadrature error or remainder. The coefficients are all positive (hildebrand87, , Sec. 8.4). These type of quadrature rules are derived under the assumption that , then, the Weierstrass approximation theorem (hildebrand87, , Sec 1.2) guarantees the existence of a polynomial such that , for an specified .
In this way can be set to minimize , and is integrated exactly.
The effective way of reducing is by increasing the polynomial degree until the residual is below . In the quadrature formula, increasing is equivalent to increasing the number of evaluation points . The remainder for the -point Gaussian quadrature satisfies

(34) |

where we see that if is a polynomial of order , the remainder vanishes and the quadrature gives the exact integral value. Further details on Gaussian quadratures can be revised in e.g. (hildebrand87, , Ch. 8), and implementation details are provided in (solin04, , Ch 4).

In , , we describe quadrature formulas for integration of , when a physical element is allowed to be curved. Then, we resource to the formula

(35) |

where the new weights and nodes

are the piled up tensor product version of vector of the corresponding one dimensional values. For example, with

, and , the new index is , and . Then, we obtain the transformed nodes and the corresponding weights as .Finally, we discuss the composite of a quadrature rule, when integration is performed over a domain defined by the union of several elements . The integrand is now required to be piecewise smooth , for . Then, we obtain

(36) |

where for each element , we have the quadrature pairs , similarly as in (35).

The polynomial spaces that we use for are based on the tensor product of one dimensional finite element spaces (solin04, , Sec. 2.2), dealII90 . As we work with piecewise smooth coefficients, we make sure that the jumps of coincide with the possibly curved element edges , such that for we have . This allow us to use quadrature rules in each individual element and guarantee convergence of the error of the numerical integration.

#### 7.2.1 Integrating weakly singular kernels

In the discretization of the Lippmann-Schwinger formulation (26), we encounter the situation where the integrand will contain both evaluation points in an element . For one-dimensional problems (), the kernel is continuous, but has a jump in the derivative at points . In the troublesome element , we can always split the integration interval and perform two separate quadrature integrations. Then, by using Gauss-type of quadratures, it is possible to avoid the evaluation of araujo+engstrom+2017 .

In higher dimensions () the kernel is weakly singular (Colton+Kress1983, , Sec. 2.3), what makes the integration in (26) more demanding. This difficulty can be overcome by specializing the quadratures as done for example in Kaneko1994 , Duan2009 , and Anand2016 . There, extra effort was spent in refining adaptively on elements , where the integrand is unbounded. Then a Nyström type of high order quadratures, combined with interpolation in polar coordinates along with other techniques were used in order to keep small to desired order. As expected, the challenge becomes more pronounced in higher dimension as can be seen in Duan2009 and in Anand2016 . In those papers the aim was to solve a scattering problem through the Lippmann-Schwinger formulation for a given incoming wave. However, our case is very different as we look for scattering resonances, where the corresponding eigensolver is computationally much more demanding than a linear solve.

### 7.3 Solution of the nonlinear eigenvalue problems

The approximation of resonances based on the DtN formulation for leads to the matrix problem in (29). The solution of this NEP is based on the solution strategy presented in araujo+engstrom+jarlebring+2017 , where we use a specialization of the Infinite Arnoldi method Jarlebring:2012:INFARNOLDI ; Jarlebring:2015:WTIARTR 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 in araujo+campos+engstrom+roman , which is a specialization of the solver in guttel17 implemented in the SLEPc library slepc05+roman .

### 7.4 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 8.2 we will numerically show the desired stability of the spectrum of the resulting discrete operator for perturbations of .

###### Remark 7.

The resulting system matrices has dimensions comparable to FE matrices: , with . However, the matrices are dense ( storage), non-symmetric, and the elements of are transcendental functions of .

### 7.5 Memory requirements

Let be the number of cells in a triangulation in space dimension , the polynomial degree of the basis functions in use, and bytes is the memory required to store a complex number in double precision. Given the number of non zeros elements in a matrix, the memory required to store it is .

In the collocation method given in Sec. 4.1, matrices are dense and we get . In turn, the FE matrices from Sec. 5.2 are sparse, and each cell in the triangulation contributes with a block of size support points. Additionally, we have scattered connections of order with neighboring cells that we omit for simplicity. A simple estimation gives . Furthermore, by assuming the division in a one dimensional partition, it is reasonable to have , and for a small, moderate and large problem respectively.

Comments

There are no comments yet.