Removal of spurious solutions encountered in Helmholtz scattering resonance computations in R^d

by   Juan Carlos Araujo Cabarcas, et al.

In this paper we consider a sorting scheme for the removal of spurious scattering resonant pairs in 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 graded material properties as well as 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 remove spurious solutions in electromagnetic and acoustic scattering resonances computations at low computational cost. Several test cases with Drude-Lorentz dielectric resonators as well as with graded material properties are considered.



There are no comments yet.


page 21

page 22


Diffraction Tomography, Fourier Reconstruction, and Full Waveform Inversion

In this paper, we study the mathematical imaging problem of diffraction ...

Weighted Delta-Tracking with Scattering

In this work, we expand the weighted delta-tracking routine to include a...

A continuation method for building invisible obstacles in waveguides

We consider the propagation of acoustic waves at a given wavenumber in a...

Computational design of locally resonant acoustic metamaterials

The so-called Locally Resonant Acoustic Metamaterials (LRAM) are conside...

Structured model order reduction for vibro-acoustic problems using interpolation and balancing methods

Vibration and dissipation in vibro-acoustic systems can be assessed usin...

Sequencing seismograms: A panoptic view of scattering in the core-mantle boundary region

Scattering of seismic waves can reveal subsurface structures but usually...

Point-based Acoustic Scattering for Interactive Sound Propagation via Surface Encoding

We present a novel geometric deep learning method to compute the acousti...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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


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


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


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


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


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


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


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


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


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

We define for the forms


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


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


where the DtN-map at is


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


The nonlinear eigenvalue problem is then as follows: Find vectors

and satisfying


for all . Note that (15) is a quadratic eigenvalue problem if is independent of and a rational eigenvalue problem for Drude-Lorentz type of materials (7).

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


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


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


Figure 1: Left) Arbitrary configuration of resonators. Right) PML stretching function.

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


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


where in and in .

For finite element computations we restrict the domain to , and define and choose similarly as in kim09 homogeneous Dirichlet boundary conditions. Formally, the finite PML problem is then: Find the eigenpairs such that


In the following, we consider a variational formulation of (21).

Find and such that for all


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

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


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

Find and such that for all


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

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


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


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


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


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


where the corresponding matrix valued function is


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


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


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 exhibit

which 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


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


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


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


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 .

The consequences of Remark 7 in a NEP solution strategy is that the matrix from (26) must be re-assembled for each new iteration , which is very expensive especially for problems of dimension .

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.