Comparison of integral equations for the Maxwell transmission problem with general permittivities
CComparison of integral equations for the Maxwelltransmission problem with general permittivities
Johan Helsing ∗ , Anders Karlsson † , and Andreas Ros´en ‡ Abstract
Two recently derived integral equations for the Maxwell transmis-sion problem are compared through numerical tests on simply con-nected axially symmetric domains for non-magnetic materials. Thewinning integral equation turns out to be entirely free from false eigen-wavenumbers for any passive materials, also for purely negative permit-tivity ratios and in the static limit, as well as free from false essentialspectrum on non-smooth surfaces. It also appears to be numericallycompetitive to all other available integral equation reformulations ofthe Maxwell transmission problem, despite using eight scalar surfacedensities.
We compare the numerical performance of two recently derived integralequation reformulations (IERs) of the Maxwell transmission problem – thatis, the problem where an incident time-harmonic electromagnetic wave isscattered from and transmitted into a bounded dielectric object. The twoIERs, which are of Fredholm’s second kind with singular integral operators,are referred to as “Dirac” and “HK 8-dens”. “Dirac” is derived in [15] byembedding the time-harmonic Maxwell’s equations into a Dirac equationand by tuning the free parameters of this equation as to optimize numericalperformance. “HK 8-dens” is derived in [12] by extending a classic IER ofthe Helmholtz transmission problem [16] via the use of certain uniquenessparameters. Both our IERs use eight unknown scalar surface densities formodeling. This is more than other popular IERs use. Typical numbers arefour [7, 18] or six [8, 19, 23]. The major advantage with our new IERs,however, is that they offer unique solutions, that is they are free from false ∗ Centre for Mathematical Sciences, Lund University, Box 118, 221 00 Lund, Sweden([email protected]). † Electrical and Information Technology, Lund University, Box 118, 221 00 Lund, Swe-den ([email protected]). ‡ Mathematical Sciences, Chalmers University of Technology and University of Gothen-burg, 412 96 G¨oteborg, Sweden ([email protected]). a r X i v : . [ phy s i c s . c o m p - ph ] S e p igenwavenumbers, for a wider range of material parameters than do otherIERs. Another advantage, from a programming point of view, is that bothIERs require only bounded integral operators with double and single layertype kernels, and in particular do not use (compact differences of) hyper-singular integral operators.Uniqueness for a wide range of parameters is an important property of anIER of a parameter-dependent partial differential equation (PDE) for manyreasons. First, one may actually be interested in solving the PDE for a widerange of parameters. Then uniqueness is obviously important. Second, evenif one is not interested in a wide range of parameters, non-uniqueness outsidethe parameter regime of interest can seriously affect the conditioning of anIER inside the parameter regime of interest. Third, theoretical studies of thesolvability of a PDE, for example the time-harmonic Maxwell’s equations inLipschitz domains, are often based on IERs. Then it is crucial that the IERhas the same solvability properties as the PDE it models.While “Dirac” and “HK 8-dens” are similar in many respects, there arealso differences. “Dirac” is derived assuming only Lipschitz regularity ofthe surface of the scattering object. “HK 8-dens” has the advantage thatits surface densities have immediate interpretations in terms of boundarylimits of physical fields. Two of these surface densities are always zero andare only needed, for uniqueness, in the main linear system that is to besolved. They are omitted in field evaluations, whose cost then are the sameas for six-density IERs.The original papers [12, 15] use mutually different notation and containnumerical examples for reduced and two dimensional versions of the IERs.The purpose of this work is to present “Dirac” and “HK 8-dens” in a unifiedand programming-friendly notation and to conduct a series of numericalexperiments for their full versions in three dimensions as to see which IERis the most efficient.We conclude the paper in Section 11 by comparing some salient proper-ties of our IERs to those of two available competitive IERs for the Maxwelltransmission problem: the “Debye” and the “DFIE” IERs [6, 23]. We state the Maxwell transmission problem. Let Ω + be a bounded domainin R with boundary surface Γ and an unbounded, connected, exterior Ω − .The outward unit normal at position r on Γ is ν . We consider time-harmonicfields with time dependence e − i t , where the angular frequency is scaled toone. The relation between time-dependent fields F ( r , t ) and complex fields F ( r ) is F ( r , t ) = (cid:60) e (cid:8) F ( r ) e − i t (cid:9) . (1)2 a) (b) 0 /20/2 Figure 1: (a) Geometry in R . Outside Γ the volume is Ω − and the wavenumber k − . Inside Γ the volume is Ω + and the wavenumber k + . The outward unitnormal is ν at r and ν (cid:48) at r (cid:48) . (b) The gray region and the solid black linesconstitute a set of points (arg( k − ) , arg( k + )) for which the Maxwell transmissionproblem has at most one solution. Circles are not included. The domains Ω − and Ω + are homogeneous and non-magnetic with wavenum-bers k − and k + . See Figure 1(a).Given incident fields E in and H in , generated in Ω − , we seek the totalelectric and magnetic fields E ( r ) and H ( r ), r ∈ Ω − ∪ Ω + , which, for ( k − , k + )with arguments in the set of points shown in Figure 1(b) and with ˆ (cid:15) suchthat ˆ (cid:15) = k /k − and ˆ (cid:15) (cid:54) = − , (2)solve Maxwell’s equations ∇ × E ( r ) = i k − H ( r ) , r ∈ Ω − ∪ Ω + , ∇ × H ( r ) = − i k − E ( r ) , r ∈ Ω − , ∇ × H ( r ) = − i k − ˆ (cid:15) E ( r ) , r ∈ Ω + , (3)except possibly at an isolated point in Ω − where the source of E in and H in is located, subject to the boundary conditionslim Ω − (cid:51) r (cid:48) → r ν × E ( r (cid:48) ) = lim Ω + (cid:51) r (cid:48) → r ν × E ( r (cid:48) ) , r ∈ Γ , (4)lim Ω − (cid:51) r (cid:48) → r ν × H ( r (cid:48) ) = lim Ω + (cid:51) r (cid:48) → r ν × H ( r (cid:48) ) , r ∈ Γ , (5) r / | r | × E sc − H sc ( r ) = o (cid:16) | r | − e (cid:61) m { k − }| r | (cid:17) , | r | → ∞ , (6) r / | r | × H sc + H sc ( r ) = o (cid:16) | r | − e (cid:61) m { k − }| r | (cid:17) , | r | → ∞ . (7)The scattered fields E sc and H sc are source free in Ω − and defined, along3ith the transmitted fields E tr and H tr , by E ( r ) = (cid:26) E in ( r ) + E sc ( r ) , r ∈ Ω − , E tr ( r ) , r ∈ Ω + . H ( r ) = (cid:26) H in ( r ) + H sc ( r ) , r ∈ Ω − , H tr ( r ) , r ∈ Ω + . (8)The incident fields satisfy ∇ × E in ( r ) = i k − H in ( r ) , r ∈ R , ∇ × H in ( r ) = − i k − E in ( r ) , r ∈ R , (9)except at the possible isolated source point in Ω − . In what follows, theMaxwell transmission problem (3)–(7) will be referred to as the MTP( k − , k + ) . Remark: since R is assumed to be non-magnetic, the permeability is ev-erywhere µ = 1. Furthermore, the local permittivities (cid:15) − and (cid:15) + of Ω − andΩ + do not enter the problem formulation, although ˆ (cid:15) of (2) can be expressedas ˆ (cid:15) = (cid:15) + /(cid:15) − . In other words: k − and k + are the only material dependentparameters needed in our formulation of the MTP( k − , k + ). k − , k + ) Let us first mention that there is yet no IER of the MTP( k − , k + ) knownthat is equivalent to the MTP( k − , k + ) itself for all pairs of complex-valued( k − , k + ), not even for all 0 ≤ arg( k − ) , arg( k + ) ≤ π . Neither is it known forwhich ( k − , k + ) there exist unique solutions to the MTP( k − , k + ). The mostcomprehensive result we are aware of, for uniqueness with Lipschitz regular Γand with µ = 1, says that the MTP( k − , k + ) has at most one solution when(arg( k − ) , arg( k + )) belongs to the set of points shown in Figure 1(b) [15,Proposition 8.2]. As for existence of solutions, the same set of points applywith the restriction that ˆ (cid:15) (cid:54) = − (cid:15) is outside of ageometry-dependent interval on the negative real axis, containing ˆ (cid:15) = − k − , k + )satisfies the conditions of [15, Proposition 8.4]. For “HK 8-dens”, theexistence of unique solutions is proven only for smooth Γ and Ω + sim-ply connected. All other IERs of the MTP( k − , k + ) in the electromag-netic literature seem to have unique solutions only under more restric-tive conditions. In particular, they can not guarantee uniqueness for any(arg( k − ) , arg( k + )) = (0 , π/ k − , k + ), which have re-ceived much attention recently, are dense-mesh and topological low-frequencybreakdown , see [23] and [7, Section 2]. Dense-mesh low-frequency breakdownrefers to catastrophic cancellation that destroys the numerical accuracy in4he computed fields in the static limit. Topological low-frequency break-down is a, perhaps, more elusive phenomenon that broadly seems to refer toan increased ill-conditioning of the integral equation itself in the quasi-staticlimit k ± →
0, with ˆ (cid:15) = ( k + /k − ) fixed, for Γ with non-zero genus. “Dirac”is proven to be free from dense-mesh low-frequency breakdown and, at large,also from topological low-frequency breakdown. See Section 7 below.Additional problems, which have not received as much attention as low-frequency breakdown but still cause numerical degradation of IERs, in-clude false near-eigenwavenumbers and false essential spectrum . An eigen-wavenumber is a pair of wavenumbers ( k − , k + ) for which the IER does nothave a unique solution. If the MTP( k − , k + ) has a unique solution we speakof a false eigenwavenumber, whereas we call it a true eigenwavenumber ifeven the MTP( k − , k + ) fails to have a unique solution. If we only considerwavenumber pairs from a set X and there is an eigenwavenumber z outsidebut close to X , then a pair x ∈ X near z which locally maximizes the con-dition number of the IER, is referred to a near-eigenwavenumber . This canbe a true or false such, depending on the nature of z .False essential spectrum may appear for certain ( k f − , k f+ ) with ˆ (cid:15) f =( k f+ /k f − ) real and negative and a merely Lipschitz regular Γ. More precisely:for a pair of wavenumbers ( k f − , k f+ ), we say that the IER has false essentialspectrum if it is not a Fredholm operator, even though the MTP( k f − , k f+ )defines a Fredholm map.Let k − = k f − and k + → k f+ in such a way that ˆ (cid:15) → ˆ (cid:15) f from aboveor from below in the complex plane. At a point ( k f − , k f+ ) where we havefalse essential spectrum, the typical numerical behavior of the IER is that ithas the same unique limit solution from above and below and this solutioncoincides with the solution to the MTP( k f − , k f+ ), while the IER is not solvablefor ( k − , k + ) = ( k f − , k f+ ). Γ From now on we assume that Γ is axially symmetric, since our numericalexamples cover such domains. Note, however, that the full 3D formulationsof our IERs are given in [12, 15]. We use both Cartesian coordinates x , y , z and cylindrical coordinates ρ , θ , z . A half-plane A in R is defined by θ = 0,a generating curve γ is defined by the intersection of Γ and A , and a generalpoint in A is r = ( ρ, z ) . (10)The outward unit normal and a tangent on γ , in A , are ν = ( ν ρ , ν z ) , (11) τ = ( ν z , − ν ρ ) . (12)5n R , the position with Cartesian basis vectors is r ≡ ( x, y, z ) = ( ρ cos θ, ρ sin θ, z ) . (13)We also use the cylindrical unit vectors ρ = (cos θ, sin θ, , (14) θ = ( − sin θ, cos θ, , (15) z = (0 , , . (16)On Γ, the unit normal is ν and the tangential unit vectors are τ and θν = ( ν ρ cos θ, ν ρ sin θ, ν z ) , (17) τ ≡ θ × ν = ( ν z cos θ, ν z sin θ, − ν ρ ) . (18)Vector fields, off from Γ, will often be expressed using ρ , θ , and zE = ρ E ρ + θ E θ + z E z , H = ρ H ρ + θ H θ + z H z . (19)The causal fundamental solution to the Helmholtz’ equation and its gradientare Φ k ( r , r (cid:48) ) = e i k | r − r (cid:48) | π | r − r (cid:48) | , (20) ∇ (cid:48) Φ k ( r , r (cid:48) ) = r − r (cid:48) π | r − r (cid:48) | (1 − i k | r − r (cid:48) | ) e i k | r − r (cid:48) | . (21) The integral equations of “Dirac” and “HK 8-dens” can both be written inthe general form ( I + G ) h ( r ) = 2 N f in ( r ) , r ∈ Γ , (22)with G = P E k + N (cid:48) − N E k − P (cid:48) . (23)Here h contains eight unknown scalar surface densities h = (cid:2) h h h h h h h h (cid:3) T , (24) f in contains the field components f in = (cid:2) ν · H in τ · H in θ · H in ν · E in τ · E in θ · E in (cid:3) T , (25)6 k is the Cauchy singular 8 × − K ν (cid:48) k K θ (cid:48) k − K τ (cid:48) k S k K ν × ν (cid:48) k − K ν k − K ν × θ (cid:48) k K ν × τ (cid:48) k S ν · ν (cid:48) k S ν · θ (cid:48) k − S ν · τ (cid:48) k K τ × ν (cid:48) k − K τ k − K τ × θ (cid:48) k K τ × τ (cid:48) k S τ · ν (cid:48) k S τ · θ (cid:48) k − S τ · τ (cid:48) k K θ × ν (cid:48) k − K θ k − K θ × θ (cid:48) k K θ × τ (cid:48) k S θ · ν (cid:48) k S θ · θ (cid:48) k − S θ · τ (cid:48) k S k − K ν (cid:48) k − K θ (cid:48) k K τ (cid:48) k S ν · ν (cid:48) k − S ν · θ (cid:48) k S ν · τ (cid:48) k − K ν × ν (cid:48) k − K ν k − K ν × θ (cid:48) k K ν × τ (cid:48) k S τ · ν (cid:48) k − S τ · θ (cid:48) k S τ · τ (cid:48) k − K τ × ν (cid:48) k − K τ k − K τ × θ (cid:48) k K τ × τ (cid:48) k S θ · ν (cid:48) k − S θ · θ (cid:48) k S θ · τ (cid:48) k − K θ × ν (cid:48) k − K θ k − K θ × θ (cid:48) k K θ × τ (cid:48) k (26)with operator entries detailed in Appendix A, and P , P (cid:48) , N , N (cid:48) are diagonalmatrices which are specified by “Dirac” and “HK 8-dens”, respectively.Once (22) is solved for h , the fields E and H in Ω − and Ω + , decomposedas in (19), can be evaluated from (8) and E sc ρ = 12 (cid:104) ˜ S ρ · ν (cid:48) k − − ˜ S ρ · θ (cid:48) k − ˜ S ρ · τ (cid:48) k − − ˜ K ρ × ν (cid:48) k − − ˜ K ρ k − − ˜ K ρ × θ (cid:48) k − ˜ K ρ × τ (cid:48) k − (cid:105) h − ,E sc θ = 12 (cid:104) ˜ S θ · ν (cid:48) k − − ˜ S θ · θ (cid:48) k − ˜ S θ · τ (cid:48) k − − ˜ K θ × ν (cid:48) k − − ˜ K θ k − − ˜ K θ × θ (cid:48) k − ˜ K θ × τ (cid:48) k − (cid:105) h − ,E sc z = 12 (cid:104) ˜ S z · ν (cid:48) k − S z · τ (cid:48) k − − ˜ K z × ν (cid:48) k − − ˜ K z k − − ˜ K z × θ (cid:48) k − ˜ K z × τ (cid:48) k − (cid:105) h − , (27) E tr ρ = 12 (cid:104) ˜ S ρ · ν (cid:48) k + − ˜ S ρ · θ (cid:48) k + ˜ S ρ · τ (cid:48) k + − ˜ K ρ × ν (cid:48) k + − ˜ K ρ k + − ˜ K ρ × θ (cid:48) k + ˜ K ρ × τ (cid:48) k + (cid:105) h + ,E tr θ = 12 (cid:104) ˜ S θ · ν (cid:48) k + − ˜ S θ · θ (cid:48) k + ˜ S θ · τ (cid:48) k + − ˜ K θ × ν (cid:48) k + − ˜ K θ k + − ˜ K θ × θ (cid:48) k + ˜ K θ × τ (cid:48) k + (cid:105) h + ,E tr z = 12 (cid:104) ˜ S z · ν (cid:48) k + S z · τ (cid:48) k + − ˜ K z × ν (cid:48) k + − ˜ K z k + − ˜ K z × θ (cid:48) k + ˜ K z × τ (cid:48) k + (cid:105) h + , (28) H sc ρ = 12 (cid:104) ˜ K ρ × ν (cid:48) k − − ˜ K ρ k − − ˜ K ρ × θ (cid:48) k − ˜ K ρ × τ (cid:48) k − ˜ S ρ · ν (cid:48) k − S ρ · θ (cid:48) k − − ˜ S ρ · τ (cid:48) k − (cid:105) h − ,H sc θ = 12 (cid:104) ˜ K θ × ν (cid:48) k − − ˜ K θ k − − ˜ K θ × θ (cid:48) k − ˜ K θ × τ (cid:48) k − ˜ S θ · ν (cid:48) k − S θ · θ (cid:48) k − − ˜ S θ · τ (cid:48) k − (cid:105) h − ,H sc z = 12 (cid:104) ˜ K z × ν (cid:48) k − − ˜ K z k − − ˜ K z × θ (cid:48) k − ˜ K z × τ (cid:48) k − ˜ S z · ν (cid:48) k − − ˜ S z · τ (cid:48) k − (cid:105) h − , (29) H tr ρ = ˆ k (cid:104) ˜ K ρ × ν (cid:48) k + − ˜ K ρ k + − ˜ K ρ × θ (cid:48) k + ˜ K ρ × τ (cid:48) k + ˜ S ρ · ν (cid:48) k + S ρ · θ (cid:48) k + − ˜ S ρ · τ (cid:48) k + (cid:105) h + ,H tr θ = ˆ k (cid:104) ˜ K θ × ν (cid:48) k + − ˜ K θ k + − ˜ K θ × θ (cid:48) k + ˜ K θ × τ (cid:48) k + ˜ S θ · ν (cid:48) k + S θ · θ (cid:48) k + − ˜ S θ · τ (cid:48) k + (cid:105) h + ,H tr z = ˆ k (cid:104) ˜ K z × ν (cid:48) k + − ˜ K z k + − ˜ K z × θ (cid:48) k + ˜ K z × τ (cid:48) k + ˜ S z · ν (cid:48) k + − ˜ S z · τ (cid:48) k + (cid:105) h + . (30)7ote that the additional zero entry for the z -components is due to axialsymmetry. Here the layer-potential entries are detailed in Appendix B andˆ k = k + /k − , h + = N (cid:48) h , h − = P (cid:48) h . (31)For “Dirac” we can improve the field evaluation by instead using the pro-jected densities h + = ( I + E k + ) N (cid:48) h , h − = ( I − E k − ) P (cid:48) h (32)as input to the field evaluation (27)–(30). The two choices (31) and (32)will always produce the same fields in Ω ± , and we discuss their uses in thefollowing subsections. P , P (cid:48) , N , N (cid:48) for “Dirac” Given the formalism above, “Dirac” specifies the diagonal matrices P = diag (cid:104) − i δ arg(ˆ c )ˆ c +1 − i δ arg(ˆ c ) 1 √ ˆ c + | ˆ c | √ ˆ c √ ˆ c | ˆ c | ˆ c + | ˆ c | ˆ (cid:15) ˆ (cid:15) +1 (cid:105) , (33) P (cid:48) = diag (cid:104) √ ˆ c + | ˆ c | √ ˆ c √ ˆ c c +1 1ˆ c +1 (cid:105) , (34) N = diag (cid:104) ˆ c ˆ c +1 − i δ arg(ˆ c ) ˆ c √ ˆ c + | ˆ c | √ ˆ c √ ˆ c c ˆ c + | ˆ c | (cid:15) +1 (cid:105) , (35) N (cid:48) = diag (cid:104) | ˆ c | √ ˆ c + | ˆ c | √ ˆ c √ ˆ c ˆ c ˆ c +1 ˆ c ˆ c +1 (cid:105) , (36)where ˆ c = 1 / ˆ k and ˆ (cid:15) = ˆ k as in (2). The matrices P and N differ in their firstelements from those used in [15, Theorem 2.3]. This corresponds to anotherchoice of the parameter β in [15, Section 8], using here β = 1 + i δ arg(ˆ k )instead of β = 1, to avoid false eigenwavenumbers when (arg( k − ) , arg( k + )) =( π/ , δ = 0 . /π . (37)Computations show that this is close enough to δ = 0 to not affect speed andaccuracy, but large enough to eliminate false eigenwavenumbers and near-eigenwavenumbers. In the original formulation of “Dirac” in R from [15],false eigenwavenumbers appear when (arg( k − ) , arg( k + )) = ( π/ , G in (23)contains Cauchy singular differences of operators. However, the particularchoice of P , P (cid:48) , N , N (cid:48) in the original “Dirac”, that is for δ = 0, makes G a compact operator on smooth Γ and as a consequence the spectrum of G has zero as its only accumulation point. This should be an advantage whenusing iterative solvers for (22). 8hen evaluating E and H with “Dirac”, one can use either (31) or (32)for h ± . The reason for preferring (32), which we use for smooth Γ in thenumerical examples of Section 10, is that, like in (25), components 1 and 5of h ± from (32) are zero. This leads to at most five non-zero densities in theevaluation of each field. However, for non-smooth Γ the numerical methodused in Section 10 is less compatible with (32) and we use the simpler (31).For “Dirac”, the densities computed with (32) satisfy P h + = N ( h − + f in ) , (38)so only one of the Cauchy integrals in (32) needs to be computed. P , P (cid:48) , N , N (cid:48) for “HK 8-dens” “HK 8-dens” specifies the diagonal matrices P = diag (cid:104) ˆ k γ η γ η − ˆ kλ ˆ (cid:15)η λ ˆ (cid:15)η − ˆ kη η ˆ kη η ˆ (cid:15)γ η γ η − ˆ (cid:15)η η λ ˆ (cid:15)η λ ˆ (cid:15)η − λ ˆ (cid:15)η λ ˆ (cid:15)η (cid:105) ,P (cid:48) = diag (cid:2) − − − − (cid:3) ,N = diag (cid:104) − γ η λ ˆ (cid:15)η η − η − γ η η − λ ˆ (cid:15)η λ ˆ (cid:15)η (cid:105) ,N (cid:48) = diag (cid:2) − ˆ c ˆ c ˆ c − ˆ c − ˆ (cid:15) − ˆ (cid:15) − − (cid:3) , (39)where γ , γ , η , and λ are uniqueness parameters whose determination isdiscussed in [12, Section 11.1]. A valid choice when (arg( k − ) , arg( k + )) =(0 ,
0) is (cid:2) γ γ η λ (cid:3) = (cid:2) ˆ (cid:15) − (cid:3) . (40)This is also a valid choice in the part of the uniqueness region of Figure 1(b)that is in the vicinity of (arg( k − ) , arg( k + )) = ( π/ , k − ) , arg( k + )) = ( π/ , k − ) , arg( k + )) = (0 , π/
2) is (cid:2) γ γ η λ (cid:3) = (cid:2) iˆ (cid:15) − − i i (cid:3) . (41)The surface densities h of (24) have, with “HK 8-dens”, the physicalinterpretations (cid:2) h h h h h h h h (cid:3) = (cid:2) σ E (cid:37) M J θ J τ σ M (cid:37) E M θ M τ (cid:3) , (42)where − i k − σ E and − i k − σ M are exterior limits of the electric and magneticvolume charge densities on Γ, (cid:37) E and (cid:37) M are the equivalent electric andmagnetic surface charge densities on the exterior side of Γ, and J θ , J τ , M θ , M τ are components of the equivalent electric and magnetic surface currentdensities on the exterior side of Γ. See [12, Remark 10.2], where it also is9hown that h = σ E and h = σ M must be zero on theoretical grounds.Therefore, the preferred choice for field evaluations with “HK 8-dens” isalways via (31), which is the Stratton–Chu representation, and (32) is neverneeded since it does not lead to any reduced numerical costs.We remark that “HK 8-dens” was discovered independently and priorto “Dirac”, and that it was only later realized that it can be written in theform (22)–(23), just like “Dirac”. However, it is not the case that “HK 8-dens” is a special case of “Dirac”, corresponding to a certain choice of Diracparameters as in [15, Section 8]. Indeed, “HK 8-dens” is not derivable fromjump matrices M and M (cid:48) as in [15, Theorem 2.3]. A key operator for the MTP( k − , k + ) is the static, k →
0, limit K d g ( r ) = p.v. (cid:90) Γ ν ( r (cid:48) ) · ∇ (cid:48) Φ ( r , r (cid:48) ) g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (43)of the acoustic double layer operator K ν (cid:48) k , appearing in the (1 ,
1) and (5 , E k . That is, K d equals the Neumann–Poincar´e operator K NP , possibly modulo a sign depending on convention [3]. Its essentialspectrum σ ess ( K d ) in the fractional Sobolev space H / (Γ), that is the setof λ for which λI − K d fails to be a Fredholm operator, is a compact subsetof the interval ( − , R , σ ess ( K d ) isnot necessarily symmetric with respect to 0 for Γ ⊂ R : see examples for thespectrum of the adjoint of K NP , σ ess ( K ∗ NP ) in H − / (Γ), and Γ with axiallysymmetric conical points in [14, Section 7.3].We map σ ess ( K d ) onto the negative real axis and defineΣ(Γ) = { ˆ (cid:15) ∈ C ; (1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) ∈ σ ess ( K d ) } . (44)For E k and our IERs, the relevant function space is H = H / (Γ) ⊕ H − / (Γ) ⊕ H − / (curl , Γ) ⊕ H / (Γ) ⊕ H − / (Γ) ⊕ H − / (curl , Γ) , (45)as discussed in [15, Section 5]. Here H − / (curl , Γ) denotes the tangentialvector fields in H − / with tangential curl in H − / , with suitable modifica-tion for non-smooth Γ. By inspection of the proof of [15, Proposition 8.5] itis immediate that the “Dirac” IER is a Fredholm operator in H if and onlyif the MTP( k − , k + ) defines a Fredholm map. By the latter we mean thatin a bounded neighborhood of Γ in R , the L -norms of the scattered andtransmitted fields are bounded by that of the incoming fields, modulo com-pact operators. Furthermore, as we shall prove in a forthcoming publication,the MTP( k − , k + ) defines a Fredholm map if and only if ˆ (cid:15) / ∈ Σ(Γ).10nother key operator for the MTP( k − , k + ) is the magnetic dipole oper-ator, K m g ( r ) = ν ( r ) × p.v. (cid:90) Γ ∇ (cid:48) Φ ( r , r (cid:48) ) × g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (46)acting on tangential vector fields g . We note that the static limit of theoperator appearing in the (3:4,3:4) and (7:8,7:8) size 2 × E k , is − K ∗ m . Moreover, the static limit of the normal derivative of theacoustic single layer potential, K ν k , appearing in the (2 ,
2) and (6 ,
6) diagonalelements, equals − K ∗ d .By Hodge decomposition of H − / (curl , Γ), it can be shown that theessential spectrum of K m is σ ess ( K d ) ∪ ( − σ ess ( K d )). The corresponding resultfor eigenvalues is in [2, Proposition 4.7] and in [17]. However, these resultsare proved for smooth Γ, on which K m is compact and σ ess ( K d ) = { } .In the diagonal blocks of the matrix E k , we therefore find operators withessential spectra σ ess ( K d ) as well as − σ ess ( K d ). To avoid the false essentialspectrum ( − σ ess ( K d )) \ σ ess ( K d ), the matrices P, P (cid:48) , N, N (cid:48) need to be chosencarefully. Consider the diagonal operator blocks of (26) for “Dirac” and“HK 8-dens” when ˆ (cid:15) is negative real and ˆ k/ i is positive, referred to as theplasmonic case in Section 10. The only operator block which can fail to bea Fredholm operator for “Dirac” is (6,6). This happens when ˆ (cid:15) ∈ Σ(Γ), thatis, if and only if the MTP( k − , k + ) itself fails to define a Fredholm map.The diagonal operator blocks of “HK 8-dens” that can fail to be Fredholmoperators are the (1,1), (2,2) and (7:8,7:8) blocks. The (1,1) and (2,2) blockseach fails to be a Fredholm operator when ˆ (cid:15) − ∈ Σ(Γ), whereas the (7:8,7:8)size 2 × (cid:15) ∈ Σ(Γ) orˆ (cid:15) − ∈ Σ(Γ). However, in the analysis of the full IERs also the non-diagonalblocks need to be taken into account, as the plots of densities in Section 10clearly show.The spectral properties of the diagonal blocks show that “Dirac” has nofalse essential spectrum, but also indicate that “HK 8-dens” may have falseessential spectrum. Section 10.3 contains numerical results that supportthis. However, since the space H is a space of mixed ± / In the quasi-static limit k ± →
0, with ˆ k = k + /k + fixed, the operator G from(23) simplifies considerably. The diagonal matrices P, P (cid:48) , N, N (cid:48) are all fixedwhereas our basic Cauchy integral operator becomes a 4 / E = − K d K , K ν × ν (cid:48) K ∗ d K , K , K , − K ∗ m − K d K , − K ν × ν (cid:48) K ∗ d K , K , K , − K ∗ m , (47)where we have written rows/columns 3:4 and 7:8 in vector/block notation.We keep indexing 1:8, and all operators are as in (26) but with k = 0. Nowall single layer operator entries vanish, as they contain a factor k ± , and theequations for the magnetic and electric fields, which correspond to the twodiagonal blocks, decouple. For spectral properties of the operators K d and K m defined in Section 6, used in this section, we refer to [5, Sections 5.1–5.2].These results generalize to Lipschitz surfaces, if we use the spaces H / (Γ)and H − / (curl , Γ) for K d and K m respectively.Let MTP(0 , , ˆ k ) denote the Maxwell transmission problem in the quasi-static limit with ˆ k fixed. The MTP(0 , , ˆ k ) amounts to two decoupled di-vergence and curl free vector fields E and H in Ω ± , decaying at infinity,having continuous tangential parts (4)–(5). We now also explicitly need torequire the Gauss jump conditionslim Ω − (cid:51) r (cid:48) → r ν · E ( r (cid:48) ) = ˆ (cid:15) lim Ω + (cid:51) r (cid:48) → r ν · E ( r (cid:48) ) , r ∈ Γ , (48)lim Ω − (cid:51) r (cid:48) → r ν · H ( r (cid:48) ) = lim Ω + (cid:51) r (cid:48) → r ν · H ( r (cid:48) ) , r ∈ Γ . (49)We note that since we consider non-magnetic materials, the jump conditionfor all components of H are the same. It is only the (5:8,5:8) diagonalblock in (47) which needs to be inverted in the quasi-static limit. Moreover,one can show that the MTP(0 , , ˆ k ) defines an invertible map if and only if(1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) / ∈ σ ( K d ).In Section 10.5 we show numerical results for condition numbers of“Dirac” as well as of “HK 8-dens” in this quasi-static limit, but restrictour analysis here to “Dirac”. For simplicity, assume δ = 0 in the definitionof P, P (cid:48) , N, N (cid:48) in Section 5.1. This makes the blocks in G corresponding to K , , K , , K , and K , vanish, and invertibility of I + G is deter-mined by the diagonal blocks. ( δ ≈ K , ≈ (cid:15) ) / (1 − ˆ (cid:15) ) ∈ σ ( K d ). Since the spectra of K d , K ∗ d and K ∗ m aresubsets of [ − , k (cid:54) = 0 and ˆ k (cid:54) = ∞ .With the source in Ω − , the end-point case ˆ k → k → ∞ as an exterior Dirichlet problem. Suitably12nterpreted, ˆ k = ∞ is a true eigenwavenumber, while ˆ k = 0 is not (assumingthat the exterior domain is connected). To analyze the behavior of “Dirac”,with δ = 0 as above, in these end-point cases, we examine which spectralpoints are used for the diagonal blocks in I + G . From [15, Equation (134)],we see that for the diagonal (2,2), (3:4,3:4) and (5:5) blocks, the spectralpoints are uniformly bounded away from [ − ,
1] as ˆ k → , ∞ . The (6,6)block is I − − ˆ (cid:15) (cid:15) K ∗ d as discussed above. The (1,1) block is I − ˆ k − k +1 K d (50)and the (7:8,7:8) block is I + ˆ k − k +1 K ∗ m . (51)Here σ ( K d ) ∩ {− , } = {− } . Furthermore σ ( K ∗ m ) ∩ {− , } = ∅ , assumingthat Ω + is simply connected. This shows that “Dirac” exhibits a false eigen-wavenumber in the quasi-static end-point case, k ± = 0, ˆ k →
0, due to the(1,1) block. If one is only interested in the exact quasi-static case, this prob-lem can be avoided by only using the (5:8,5:8) block of I + G for solving forthe electric fields as discussed above. (Or even simpler IERs available in thiscase.) If one is also interested in near quasi-statics and ˆ k ≈
0, assuming thatΩ + is simply connected, we can tune the Dirac parameters r, β, γ, α (cid:48) , β (cid:48) , γ (cid:48) from [15, Section 8] as follows. Rather than choosing α (cid:48) = β (cid:48) = 1 / ˆ k as donefor “Dirac”, we choose α (cid:48) = | ˆ k | / ˆ k and β (cid:48) = 1 / ( | ˆ k | ˆ k ). With these choices,the spectral point for the (1,1) block stays uniformly bounded away from[ − , k → + , where +1 / ∈ σ ( K ∗ m ). Pre-conditioning similarly to [15, Theorem 2.3], we arrive at a version of “Dirac”given by P = (cid:20) β √ | ˆ c | ˆ c + | ˆ c | β √ ˆ c + | ˆ c | √ ˆ c √ ˆ c | ˆ c | ˆ c + | ˆ c | ˆ (cid:15) | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c (cid:21) , (52) P (cid:48) = (cid:104) √ | ˆ c | √ ˆ c + | ˆ c | √ ˆ c √ ˆ c (cid:15) +1 1 | ˆ c | ˆ c | ˆ c | ˆ c (cid:105) , (53) N = (cid:20) √ | ˆ c | β | ˆ c | / ˆ c ˆ c √ ˆ c + | ˆ c | √ ˆ c √ ˆ c c ˆ c + | ˆ c | | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c | ˆ c (cid:21) , (54) N (cid:48) = (cid:104)(cid:112) | ˆ c | | ˆ c | √ ˆ c + | ˆ c | √ ˆ c √ ˆ c (cid:15) +1 (cid:105) , (55)where β = 1 − i δ arg(ˆ c ) is chosen as in Section 5.1. For simply connected Ω + ,with k ± ≈ k ≈
0, this IER will be better conditioned than “Dirac”.See Section 10.5 for numerical results. Further tuning of the six complexDirac parameters from [15, Section 8] may be needed when Ω + is not simplyconnected. 13 Surface plasmon standing waves and plasmons
Let, for the moment, R be divided into two, not necessarily bounded, do-mains Ω and Ω separated by a surface Γ and with real-valued permittivities (cid:15) and (cid:15) such that (cid:15) > (cid:15) < Surface plasmon waves , surface plas-mon standing waves , and quasi-static plasmons are then particular typesof electromagnetic fields that may appear as solutions to the MTP( k − , k + ).These fields are briefly described here from a classical electrodynamics pointof view and they appear in the numerical examples of Section 10.Surface plasmon waves (SPWs) are surface waves that can propagatealong Γ. A necessary extra condition for their existence on planar Γ isˆ (cid:15) ≡ (cid:15) /(cid:15) < − λ spw = λ (cid:112) − / | ˆ (cid:15) | (56)and decay exponentially in the directions normal to Γ. Here λ = 2 π/k isthe free space wavelength in Ω . The SPWs can, further, propagate withoutattenuation along surfaces that are invariant only in the propagation direc-tion. They can also propagate along surfaces with a finite radius of curvaturein the propagation direction, but are then attenuated due to radiation. Theradiation increases with the ratio of λ spw to the radius of curvature.Consider now, with notation as in Section 2, a bounded object Ω + andassume that the wavenumber k − is positive real and that ˆ (cid:15) < −
1, to enableSPWs along Γ. At certain k − the SPWs form standing waves. Such SPWsare here called surface plasmon standing waves (SPSWs). In some literaturethe name surface plasmon resonances (SPRs) is used to emphasize theirresonant nature. However, SPR is not a direct synonym for SPSW since itis also used for other plasmonic phenomena. Due to radiation, the SPSW isa damped resonant field and the ( k − , k + ), for which it appears, is close to atrue eigenwavenumber outside the set of points shown in Figure 1(b). Whenan SPSW is excited by an incident field, its amplitude becomes large, whichresults in a large scattering cross section. This is one reason why SPSWsare of interest in optics.A quasi-static plasmon is an electromagnetic field that approaches aneigenfield, referred to as a static plasmon , as k ± →
0. It appears aroundobjects that are much smaller than the wavelength λ − . The quasi-staticplasmon is also denoted surface plasmon (SP) in optics, a name that maylead to confusion since it is sometimes used as a synonym for SPSW. Aquasi-static plasmon is a resonant field, but in contrast to SPSWs it is nota standing wave field. For smooth Γ there is an infinite discrete set of ˆ (cid:15) corresponding to static plasmons. The electric field of a static plasmonis distributed so that the electric energies in Ω − and Ω + exactly canceleach other, while the magnetic field is zero. Indeed, Green’s identity showsthat (cid:82) Ω + | E | dx = − ˆ (cid:15) (cid:82) Ω − | E | dx so that, in accordance with Poincar´e’s14ariational principle [3, Theorem 2.4], the static plasmon corresponds to theeigenvalue (1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) for K d .Quasi-static and static plasmons can be classified as being bright or dark depending on whether they can be excited by a uniform incident field ornot. A bright plasmon radiates as an electric dipole and its far-field isvery large considering the object is small compared to λ − . A dark plasmonradiates as a higher-order multipole and its far-field is weak. See [2, Lemma5.3], where the full multipole expansion of far-fields is given. Furthermore,there is a close connection between bright static plasmons and the imaginarypart of an object’s limit polarizability [13, 21]. For an object with sharpcorners, edges, or points there is a continuous, possibly punctured, intervalof ˆ (cid:15) where a special type of static plasmons, with infinite (but cancelling)electric energies in Ω − and Ω + , can occur together with a partially embeddeddiscrete set of regular static plasmons [10, 14, 21]. This section reviews shapes of scattering objects and discretization schemesthat are used for the numerical tests in Section 10.
Examples with smooth Γ come from a generating curve γ parameterized as r ( s ) = (1 + α sin(5 s ))(cos( s ) , sin( s )) , s ∈ [ − π/ , π/ , (57)where α is a shape parameter. The choice α = 0 corresponds to Γ being theunit sphere. With α = 0 .
25 the shape of Γ resembles a “rotated starfish”.Examples with non-smooth Γ come from a γ parameterized as r ( s ) = sin( πs ) (sin((0 . − s ) α ) , cos((0 . − s ) α )) , s ∈ [0 , . , (58)where α is the opening angle of the conical point at the origin. With α > π ,the shape of Γ resembles a “tomato”. See Figure 1(a) for an illustration with α = 31 π/
18. With α < π the shape of Γ resembles a “drop with a sharptip”. The generating curve γ of (58) has previously been used for numericalexamples in [12, 14]. The integral equation (22) is solved on axially symmetric Γ using high-order Fourier–Nystr¨om discretization. An azimuthal Fourier transform isfirst applied to (22), yielding a sequence of modal problems for the Fouriercoefficients h ( n ) of h ( I + G ( n ) ) h ( n ) ( r ) = 2 N f in( n ) ( r ) , r ∈ γ , n = 0 , ± , ± , . . . , (59)15hich are solved independently using high-order panel-based Nystr¨om dis-cretization. The modal solutions h ( n ) are then synthesized to give the fullsolution h .High-order panel-based Fourier–Nystr¨om discretization for solving inte-gral equations modeling scattering problems on axisymmetric surfaces wasmade popular by Young, Hao, and Martinsson in 2012 [24]. Since then,several authors have worked on extensions of such schemes and related is-sues. The aim has been to include a broader range of integral operators,to reach higher achievable accuracy, to evaluate near fields more efficiently,and to cope with problems related to wavenumbers with large imaginaryparts [7, 18, 11, 1]. Our version of the Fourier–Nystr¨om scheme is theone used in [12]. It uses a combination of 16-point and 32-point underlyingGauss–Legendre quadrature, a variety of explicit kernel-splits, semi-analyticproduct integration computed on the fly, and is compatible with the recur-sively compressed inverse preconditioning method (RCIP) [9]. The RCIPaccelerates and stabilizes Nystr¨om discretization in the presense of singularboundary points on γ , such as corners.
10 Numerical examples
The numerical efficiency of “Dirac” and “HK 8-dens” is now compared. Thatis, we solve the discretized modal systems (59) and compute fields via (27)–(30) with P , P (cid:48) , N , N (cid:48) as in Sections 5.1 and 5.2. We also compute conditionnumbers of the modal systems. Three material parameter cases are used: • the positive dielectric case , where ˆ k = 1 . k − is positive real; • the plasmonic case , where ˆ k = i √ . k − is positive real; • the reverse plasmonic case , where ˆ k = (cid:0) i √ . (cid:1) − and k + is positivereal.These parameter cases are taken from previous work on time-harmonictransmission problems [4, 12, 15].In all examples involving field images, the incident field is a linearlypolarized plane wave E in ( r ) = x e i k − z with x = (1 , , y = 0, where the field components E θ , H ρ , and H z arezero due to symmetry. To save space in the examples, we only show E ρ and H θ . Generally, these two seem to exhibit more pronounced field patternsthan the omitted component E z .When Γ is non-smooth, it may happen that ˆ k = i √ . k = (cid:0) i √ . (cid:1) − correspond to that the MTP( k − , k + ) does not define a Fred-holm map. We then compute limit solutions as ˆ (cid:15) approaches − . − / . H from (45) and are given16 downarrow superscript. For example, the limit of the field E is denoted E ↓ .Our codes are implemented in Matlab , release 2018b, and executedon a workstation equipped with an Intel Core i7-3930K CPU and 64 GBof RAM. Large linear systems are solved iteratively using GMRES witha stopping criterion threshold of machine epsilon ( (cid:15) mach ) in the estimatedrelative residual. When assessing the accuracy of computed field quantitieswe adopt a procedure where to each numerical solution we also computean overresolved reference solution, using roughly 50% more points in thediscretization of the system under study. The absolute difference betweenthese two solutions is denoted the estimated absolute error . In examplesinvolving field images, the fields are computed at 10 target points on arectangular Cartesian grid in the computational domains shown.Subsections 10.1 and 10.5, on eigenwavenumbers and the quasi-staticlimit, concern the unit sphere. An advantage with doing experiments onspheres is that semi-analytical results given by Mie theory can be used forverification. The remaining subsections concern less trivial object shapesand there we limit ourselves mainly to the plasmonic case. The positivedielectric case is not deemed challenging enough for thorough testing. Theperformances of “Dirac” and “HK 8-dens” in the reverse plasmonic case israther similar to their respective performances in the plasmonic case – apartfrom the new risk of being close to a true eigenwavenumber and, for “HK8-dens”, also close to a false eigenwavenumber or near-eigenwavenumber. We compute condition numbers of system matrices of the discretized modalintegral equations (59) and with Γ being the unit sphere. Results for theazimuthal modes n = 0 , ,
10 are shown in Figure 2. A number of 768 dis-cretization points are placed on a grid on γ , making the total system size6144 × ,
300 data points are used to capture the rapid varia-tions in the condition numbers as k − and k + are swept through the intervals[0 ,
10] and [0 , k − and k + . Note that the regularly recurringhigh peaks that are common to Figures 2(c,d) correspond to true near-eigenwavenumbers, while the additional peaks in Figure 2(d) correspond tofalse near-eigenwavenumbers of “HK 8-dens”, as discussed in the last para-graph of Section 3. The true near-eigenwavenumbers are associated withSPSWs. Note also that the 14 eigenwavenumbers with k + ∈ [0 , (a) (b) (c) (d) (e) (f) Figure 2:
The unit sphere. Condition numbers of system matrices of the dis-cretized modal integral equations (59) with azimuthal indices n = 0 , , andwavenumber k − ∈ [0 , or k + ∈ [0 , : (a,b) the positive dielectric case with ˆ k = 1 . ; (c,d) the plasmonic case with ˆ k = i √ . ; (e,f) the reverse plas-monic case with ˆ k = (cid:0) i √ . (cid:1) − ; (a,c,e) “Dirac”; (b,d,f) “HK 8-dens”. in Figure 2(e), are true eigenwavenumbers in the reverse plasmonic case, asconfirmed to at least 13 digits by comparison with semi-analytical results.In contrast, “HK 8-dens” exhibits here, in addition, around 25 false eigen-wavenumbers and false near-eigenwavenumbers. See Figure 2(f). Elevenof the 14 peaks in Figure 2(e) form a periodic pattern and correspond toeigenwavenumbers with eigenfields that are SPSWs, whereas the other three18 a) -1 -0.5 0 0.5 1 x -1-0.500.51 z -6-4-20246 (b) -1 -0.5 0 0.5 1 x -1-0.500.51 z -2-1.5-1-0.500.511.52 (c) -1 -0.5 0 0.5 1 x -1-0.500.51 z -15.5-15-14.5-14-13.5-13-12.5-12-11.5 (d) -1 -0.5 0 0.5 1 x -1-0.500.51 z -15.5-15-14.5-14-13.5-13-12.5 (e) -1 -0.5 0 0.5 1 x -1-0.500.51 z -15.5-15-14.5-14-13.5-13-12.5-12 (f) -1 -0.5 0 0.5 1 x -1-0.500.51 z -15.5-15-14.5-14-13.5-13-12.5-12 Figure 3:
Field images for scattering from the “rotated starfish” in theplasmonic case and with E in ( r ) = x e i6 z : (a,c,e) the scattered/transmitted E ρ ( r , -field; (b,d,f) the scattered/transmitted H θ ( r , -field; (c,d) log ofestimated absolute field error for “HK 8-dens”; (e,f) the same for “Dirac”. peaks correspond to eigenwavenumbers with eigenfields that are not boundto the surface. 19 We compute images of scattered and transmitted fields E sc ρ ( r , H sc θ ( r , E tr ρ ( r , H tr θ ( r , α = 0 .
25. A number of 768 discretization points are usedfor each integral operator on γ in (59).We first test the positive dielectric case (ˆ k = 1 .
5) with k − = 10. Thisis a simple problem and both “Dirac” and “HK 8-dens” give fields with anestimated absolute precision of at least twelve digits – also close to Γ. Werefrain from showing images. The only noticeable difference between themethods is the number of GMRES iterations needed for full convergence.The “Dirac” system converges to an estimated relative residual of (cid:15) mach in111 iterations while the “HK 8-dens” system is a bit better and only needs68 iterations.We next test the plasmonic case (ˆ k = i √ . k − = 6. Figure 3shows that the accuracy achieved by “Dirac” and “HK 8-dens” is similaralso in this case, although the spectral properties of “Dirac” are now muchbetter than those of “HK 8-dens”. The condition number of the systemmatrix in the discretized “Dirac” system is 9 . · while the correspondingcondition number for “HK 8-dens” is 1 . · and this is reflected in thenumber of GMRES iterations needed for full convergence. The “Dirac”system converges to an estimated relative residual of (cid:15) mach in 170 iterationswhile the “HK 8-dens” system needs 623 iterations.The symmetry of E in is such that only the two azimuthal modes n = − n = 1 are present. The Fourier coefficients of the surface densitiesof these modes, h (1) and h ( − of (59), are either identical or have oppo-site signs. Note also the SPWs propagating along Γ in the images of Fig-ure 3(a,b). Their wavelength is roughly 15% shorter than the wavelengthgiven by (56). We repeat the experiments of Section 10.2, but now for the “tomato”, thatis, with Γ non-smooth and generated by γ of (58) with α = 31 π/
18 asillustrated in Figure 1(a). A number of 576 discretization points are usedfor each integral operator on γ in (59).In the positive dielectric case, and with a larger wavenumber k − = 18 asto compensate for the “tomato” being smaller than the “rotated starfish”,the results for both methods are even (marginally) better than in the pre-vious example. The estimated pointwise precision in the field images isbetween twelve and 13 digits (no images shown). The “Dirac” system needs96 GMRES iterations for full convergence while the “HK 8-dens” systemneeds 87 iterations. We conclude that, thanks to the RCIP method for deal-ing with the conical point, the non-smooth “tomato” is as simple as the20 a) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -4-3-2-101234 (b) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -6-4-20246 (c) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -15.5-15-14.5-14-13.5-13-12.5-12-11.5-11 (d) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -15.5-15-14.5-14-13.5-13-12.5-12-11.5 (e) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -15.5-15-14.5-14-13.5-13-12.5-12-11.5 (f) -0.5 0 0.5 x -0.4-0.200.20.40.60.811.2 z -15.5-15-14.5-14-13.5-13-12.5-12 Figure 4:
Field images for scattering from the “tomato” in the plasmoniccase and with E in ( r ) = x e i5 z : (a,c,e) the scattered/transmitted E ρ ( r , -field; (b,d,f) the scattered/transmitted H θ ( r , -field; (c,d) log of estimatedabsolute field error for “HK 8-dens”; (e,f) The same for “Dirac”. The colorbarrange in (a) is set to [ − . , . . smooth “rotated starfish” from a numerical point of view.In the plasmonic case and with k − = 5 we observe some very interestingfeatures. The wavenumber ratio corresponds to that (1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) is in − σ ess ( K d ), but not in σ ess ( K d ). Given the discussion about essential spec-21 -10 -5 -10 -10 -5 -5 -10 -5 -5 -10 -5 -5 -10 -5 -10 -10 -5 -10 -5 -10 -5 -10 -5 -5 -10 -5 -5 -10 -5 -5 -10 -5 -5 -10 -5 -10 -5 -10 -5 -10 -5 -5 -10 -5 -5 Figure 5:
Asymptotics for the “tomato” in the plasmonic case. The densities h i , i = 1 , . . . , , for azimuthal mode n = 1 , as functions of the arc lengthdistance along γ to the conical point at the origin. Row 1-2: “HK 8-dens”.Row 3-4: “Dirac”. tra in Section 6, it may not come as a complete surprise that “HK 8-dens”exhibits false essential spectrum in this case, even though correct limit so-lutions can be computed, see the last paragraph of Section 3. “Dirac”, onthe other hand, is free from this problem and correct fields can be computedwithout a limit process. Figure 4 shows that “Dirac” clearly achieves betterfield accuracy than “HK 8-dens” in this case. In terms of GMRES conver-gence the difference is even greater: “Dirac” needs 88 iterations while “HK8-dens” needs 326 iterations. There are SPWs propagating along Γ with awavelength that is roughly 20% shorter than the wavelength given by (56).It is also enlightening to inspect the asymptotics of h close to the con-ical tip of the “tomato” in the plasmonic case. Thanks to the symmetryof E in , as discussed in the last paragraph of Section 10.2, it is enough tostudy the eight densities (Fourier coefficients) contained in the modal solu-tion h (1) . Figure 5 shows these eight densities both for “Dirac” and “HK8-dens”. The individual densities are denoted h i , i = 1 , . . . ,
8, with the az-imuthal index omitted. Note that the densities h and h of “HK 8-dens”22 a) -0.6 -0.4 -0.2 0 0.2 0.4 x z -8-6-4-202468 (b) -0.6 -0.4 -0.2 0 0.2 0.4 x z -1-0.500.51 (c) -0.6 -0.4 -0.2 0 0.2 0.4 x z -15-14-13-12-11-10 (d) -0.6 -0.4 -0.2 0 0.2 0.4 x z -15.5-15-14.5-14-13.5-13-12.5 Figure 6:
Field images for scattering from the “drop with a sharp tip” inthe plasmonic case and with E in ( r ) = x e i5 z : (a,c) the scattered/transmitted E ↓ ρ ( r , -field; (b,d) the scattered/transmitted H ↓ θ ( r , -field; (c,d) log ofestimated absolute field error for “Dirac”. The colorbar range in (a) is set to [ − . , . and in (b) set to [ − . , . . are approximately zero, as they should be according to Section 5.2.The strongest singularity observed in Figure 5 is h i ∝ s − . , (60)where s is the arc length distance along γ to the conical point at the origin.For “HK 8-dens”, the singularity of (60) is observed for the densities h , h ,and h , corresponding to (cid:37) E , M θ , and M τ according to (42). For “Dirac”,the singularity of (60) is observed only for the density h . Note that eachof these functions is more singular than H / (Γ), as allowed by the functionspace H . We repeat some of the experiments of Section 10.3, but shrink the openingangle of the conical point to α = 5 π/
18, so that the shape of Γ now resembles23 -10 -5 -10 -10 -5 -10 -5 -10 -5 -10 -5 -10 -10 -5 -10 -5 -10 -5 -10 -5 -10 -5 -10 -5 -10 -5 -10 -5 -5 -10 -5 -10 -5 -5 -10 -5 -5 Figure 7:
Asymptotics for the “drop with a sharp tip” in the plasmonic case.Densities displayed as in Figure 5. Row 1-2: “HK 8-dens”. Row 3-4: “Dirac”. that of a “drop with a sharp tip”. A number of 576 discretization points areagain used for each integral operator on γ in (59).Field images for the plasmonic case along with error estimates for “Dirac”are shown in Figure 6. No SPWs are excited along Γ. A number of 83 itera-tions were needed for full GMRES convergence. The corresponding numberfor “HK 8-dens” is 314 iterations. Note that (1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) is now in theessential spectrum of K d and that E ↓ ρ and H ↓ θ are oscillatory and unboundedat the origin. The colorbar ranges in Figure 6(a,b) are therefore restrictedto the most extreme field values away from the origin. The estimated point-wise absolute error close to the origin is, of course, now larger than than inprevious examples with everywhere bounded fields.Figure 7 is analogous to Figure 5 and exhibits the same general features.The strongest singularity observed is now h i ∝ s − . . , (61)clearly visible for h of “Dirac” and for h , h , and h of “HK 8-dens”.As we are now in the essential spectrum, we do not expect the densities tobelong to H , and indeed the above mentioned densities just fail to belong24 (a) -2 -1 0 1 210 (b) -2 -1 0 1 210 (c) -2 -1 0 1 210 (d) Figure 8:
The unit sphere. Condition numbers of system matrices of the dis-cretized modal integral equations (59) with azimuthal index n = 0 in the quasi-static limit k + = k − = 0 . The plasmonic case is x ∈ ( − , and the positivedielectric case is x ∈ ( − , − ∪ (1 , . (a) “HK 8-dens”; (b) “Dirac”; (c) The“Dirac” (5:8,5:8) diagonal block; (d) The “Dirac” version (52) – (55) . to H − / (Γ). Moreover, the densities h and h , and for “Dirac” also h ,just fail to belong to H / (Γ). This second strongest singularity observedis s − . . . Again h and h are zero for “HK 8-dens”, modulorounding errors. Similar to Section 10.1 we plot condition numbers of our IERs for the unitsphere to detect possible false eigenwavenumbers, but also to get a notion ofhow well-conditioned the two IERs under study are. This time we considerthe quasi-static limit k ± → x = (1 + ˆ (cid:15) ) / (1 − ˆ (cid:15) ) in Figure 8. With a slightabuse of notation, we keep speaking of (true/false/near) eigenwavenumbers.Here x < − k > x = − − < x < k/ i > x = 0 corresponds to the essential spectrum,0 < x < < ˆ k/ i < x = 1 corresponds25o an exterior Neumann problem, and x > < ˆ k < K d and thepeaks correspond to the static plasmons discussed in Section 8. The peakat x = − / < x <
1. The full “Dirac” shows only a false eigen-wavenumber in the limit x = 1. As shown in Figure 8(d), this peak can beremoved by using the version of “Dirac” specified by (52)–(55). Note that“Dirac” specified by (33)–(36), is better conditioned near x = −
1, wherethe MTP(0 , , ˆ k ) has a dark static plasmon which for the unit sphere is themonopole field E = r / | r | , H = for | r | >
1, and E = H = for | r | <
11 Conclusions and discussion
Our numerical results show that “Dirac” wins over “HK 8-dens” in almostall tests. “Dirac” does not have any false near-eigenwavenumbers for anypassive materials, whereas “HK 8-dens” exhibits such in the plasmonic case,and even false eigenwavenumbers in the reverse plasmonic case. We haveseen that “HK 8-dens” exhibits false essential spectrum on domains withcorners, and false eigenwavenumbers in the quasi-static limit. “HK 8-dens”also requires more than three times as many GMRES iterations than “Dirac”in the plasmonic cases.Since “Dirac” clearly is a competitive IER of the MTP( k − , k + ), we endby comparing it to other available IERs. We limit our discussion to non-magnetic materials, although “Dirac” was formulated for general materialsin [15], which we plan to cover in a forthcoming publication. The aspects ofthe IERs that we base our discussion on are the following.(A) Can it be used for Lipschitz regular/piecewise smooth surfaces Γ?(B) What types of operators does it employ?(C) How fast is it to compute all the scattered and transmitted fields?(D) How does it behave in the quasi-static limit?(E) How stable is it? Does it have false eigenwavenumbers, false near-eigenwavenumbers, or false essential spectrum?We now make a comparison with the “Debye” formalism in [6, 7] and the“DFIE” formalism in [23] which, according to [23, Conclusions], are thecurrent leading IERs when it comes to well-posedness for a wide range ofpassive materials. Using formalism from [15, 22], “Dirac” uses Cauchy type26ntegral representations for the Dirac equation D F = i kF for the electro-magnetic multivector field F , whereas “Debye” employs a Dirac multivectorpotential F = D G + i kG , leading to a component-wise Helmholtz equation∆ G + k G = 0, for which classical layer potentials apply. The electric andmagnetic Gauss equations translate to certain conditions on these boundarylayers, which via Hodge decompositions of tangential vector fields allows oneto eliminate all but two scalar boundary densities for each domain. In thescattering situation of the present paper, this yields a “Debye” IER withfour scalar densities. It is known that (D) it does not have dense-meshlow-frequency breakdown, and no topological low-frequency breakdown ifcomplemented by extra equations. A drawback of “Debye” is that imple-menting the Hodge projections requires the numerical solution of a secondorder surface Laplace equation on Γ, which (A) hampers the applicabilityto non-smooth Γ. According to [7, Equation (2.19)], “Debye” also (B) usescompositions of two integral operators. It is unclear to us, what has beenproven for “Debye” about (E). In [7, Conclusions], it is stated that “Debye”is invertible for all passive materials. However, this is impossible for anyIER due to the true eigenwavenumbers shown in Figure 2(e). Moreover, themain uniqueness result for “Debye”, [6, Theorem 3.2], assumes a positivereal part on permittivity and relies on [20, Theorem 69], which assumespermittivities in the first quadrant [20, p. 261]. These assumptions do notcover all passive non-magnetic materials, which correspond to the square[0 , π/ in Figure 1(b).“DFIE” is formulated in [23]. Like “Debye”, it employs potentials toreduce Maxwell’s equations to vector Helmholtz equations. With “DFIE”,the fields E and H are calculated independently of each other, by solvingan IER with six scalar densities for each field. Again it is known for “DFIE”that (D) it does not suffer from low-frequency breakdown and (E) it doesnot have false eigenwavenumbers for k − > ≤ arg(ˆ (cid:15) ) < π , but it hasfalse near-eigenwavenumbers at (arg( k − ) , arg( k + )) = (0 , π/ C regular. However, by straightforward adaption to fractionalSobolev spaces and Lipschitz regular Γ as in [15], it appears that (A) “DFIE”is applicable to non-smooth scattering. Inspecting the integral operatorsthat “DFIE” employs in [23, Appendix D], we see that (B) it uses singleand double layer type boundary operators like “Dirac”, with one exception.The exception is K in [23, Equation (89)], which is a compact differenceof hypersingular operators, and somewhat more numerically challenging toimplement [18]. To compare (C) how fast “DFIE” and “Dirac” are, since theoperators employed are roughly the same and having a good solver in mindwhich is linear in the number of non-zero operator blocks, we count these.From [23, Equation (87)] we have 2(36 −
4) = 64 non-zero operator blocksfor “DFIE”, keeping in mind that we need to solve two 6 × E and H . From (26) we see that “Dirac” uses 64 −
14 = 5027on-zero operator blocks. For field evaluations, [23, Equations (36) and(38)] show that “DFIE” requires the evaluation of three double layer andthree single layer potentials, for each component of the fields. For “Dirac”,using the projected densities (32), which are fast to compute on Γ, each fieldcomponent requires in general the evaluation of three double layer and twosingle layer potentials.From this comparison, our conclusion is that “Dirac” is the most com-petitive IER, except possibly for (C), where further testing needs to be un-dertaken to accurately judge which of “Dirac” and “DFIE” is the faster IER.Nevertheless, “Dirac” has the advantage of (A) applying to any Lipschitzregular Γ, (B) only using bounded integral operators of single and doublelayer type with explicit kernels, (E) not having any false eigenwavenumbers,false near-eigenwavenumbers or false essential spectrum, (D) not even in thegeneral quasi-static limit, or in the special quasi-static x = +1 limit if wechoose the Dirac parameters appropriately. The reason why “Dirac” is sosuccessful in avoiding false spectra, is that the choices of parameters makeits invertibility depend only on the (6,6) diagonal block involving K d . Thusthe false spectrum of K m , of which only one Hodge component is needed forthe MTP( k − , k + ), is avoided.In [2, Important Remark p. 123] it is stated that σ ( K m ) \ σ ( K d ) con-tributes to “higher-order resonances” in the near quasi-static limit. How-ever, plots of the condition numbers for “Dirac” on the sphere near thequasi-static limit (not shown) are very similar to Figure 8(b). In particular,there are no visible traces of σ ( K m ) in terms of near-eigenwavenumbers for0 < x < Acknowledgement
This work was supported by the Swedish Research Council under contract2015-03780. ***
Appendix ***
A The entries of E k The entries of the matrix E k of (26) involve two families of integral operatorsdenoted S k and K k . A given member in an operator family, S αk or K αk , isdefined by its superscript α , which can be a constant, a unit vector, a scalarproduct of unit vectors, or a cross product of unit vectors. Specifically we28ave for S k acting on a general density gS k g ( r ) = i k (cid:90) Γ Φ k ( r , r (cid:48) ) g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (62)for S u · v (cid:48) k , where u and v are unit vectors, S u · v (cid:48) k g ( r ) = i k (cid:90) Γ u ( r ) · v ( r (cid:48) )Φ k ( r , r (cid:48) ) g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (63)and for K u k and K u × v (cid:48) k K u k g ( r ) = p.v. (cid:90) Γ u ( r ) · ∇ (cid:48) Φ k ( r , r (cid:48) ) g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (64) K u × v (cid:48) k g ( r ) = p.v. (cid:90) Γ ( u ( r ) × v ( r (cid:48) )) · ∇ (cid:48) Φ k ( r , r (cid:48) ) g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ . (65)When coding, particularly when Γ is axially symmetric and when az-imuthal Fourier transforms are to be implemented, it is helpful to have S αk and K αk in the form S αk g ( r ) = i k (cid:90) Γ s α ( r , r (cid:48) )2 π | r − r (cid:48) | e i k | r − r (cid:48) | g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (66) K αk g ( r ) = p.v. (cid:90) Γ d α ( r , r (cid:48) )2 π | r − r (cid:48) | (1 − i k | r − r (cid:48) | ) e i k | r − r (cid:48) | g ( r (cid:48) ) dΓ (cid:48) , r ∈ Γ , (67)where s α ( r , r (cid:48) ) and d α ( r , r (cid:48) ) are static kernel factors expressed in terms ofquantities introduced in (10-12) and the azimuthal angle θ .Here follow s α ( r , r (cid:48) ) for the ten operators S αk of E k s ( r , r (cid:48) ) = 1 , (68) s ν · ν (cid:48) ( r , r (cid:48) ) = ν ρ ν (cid:48) ρ cos( θ − θ (cid:48) ) + ν z ν (cid:48) z , (69) s ν · θ (cid:48) ( r , r (cid:48) ) = ν ρ sin( θ − θ (cid:48) ) , (70) s ν · τ (cid:48) ( r , r (cid:48) ) = ν ρ ν (cid:48) z cos( θ − θ (cid:48) ) − ν z ν (cid:48) ρ , (71) s τ · ν (cid:48) ( r , r (cid:48) ) = ν z ν (cid:48) ρ cos( θ − θ (cid:48) ) − ν ρ ν (cid:48) z , (72) s τ · θ (cid:48) ( r , r (cid:48) ) = ν z sin( θ − θ (cid:48) ) , (73) s τ · τ (cid:48) ( r , r (cid:48) ) = ν z ν (cid:48) z cos( θ − θ (cid:48) ) + ν ρ ν (cid:48) ρ , (74) s θ · ν (cid:48) ( r , r (cid:48) ) = − ν (cid:48) ρ sin( θ − θ (cid:48) ) , (75) s θ · θ (cid:48) ( r , r (cid:48) ) = cos( θ − θ (cid:48) ) , (76) s θ · τ (cid:48) ( r , r (cid:48) ) = − ν (cid:48) z sin( θ − θ (cid:48) ) . (77)Here follow d α ( r , r (cid:48) ) for the 15 operators K αk of E k d ν ( r , r (cid:48) ) = ν · ( r − r (cid:48) ) + ν ρ ρ (cid:48) (1 − cos( θ − θ (cid:48) )) , (78)29 ν (cid:48) ( r , r (cid:48) ) = ν (cid:48) · ( r − r (cid:48) ) − ν (cid:48) ρ ρ (1 − cos( θ − θ (cid:48) )) , (79) d τ ( r , r (cid:48) ) = τ · ( r − r (cid:48) ) + ν z ρ (cid:48) (1 − cos( θ − θ (cid:48) )) , (80) d τ (cid:48) ( r , r (cid:48) ) = τ (cid:48) · ( r − r (cid:48) ) − ν (cid:48) z ρ (1 − cos( θ − θ (cid:48) )) , (81) d θ ( r , r (cid:48) ) = ρ (cid:48) sin( θ − θ (cid:48) ) , (82) d θ (cid:48) ( r , r (cid:48) ) = ρ sin( θ − θ (cid:48) ) , (83) d ν × ν (cid:48) ( r , r (cid:48) ) = ( ν (cid:48) ρ τ · r − ν ρ τ (cid:48) · r (cid:48) ) sin( θ − θ (cid:48) ) , (84) d ν × θ (cid:48) ( r , r (cid:48) ) = − τ · ( r − r (cid:48) ) + ( τ · r + ν ρ z (cid:48) )(1 − cos( θ − θ (cid:48) )) , (85) d ν × τ (cid:48) ( r , r (cid:48) ) = ( ν ρ ν (cid:48) · r (cid:48) + ν (cid:48) z τ · r ) sin( θ − θ (cid:48) ) , (86) d τ × ν (cid:48) ( r , r (cid:48) ) = − ( ν (cid:48) ρ ν · r + ν z τ (cid:48) · r (cid:48) ) sin( θ − θ (cid:48) ) , (87) d τ × τ (cid:48) ( r , r (cid:48) ) = − ( ν (cid:48) z ν · r − ν z ν (cid:48) · r (cid:48) ) sin( θ − θ (cid:48) ) , (88) d τ × θ (cid:48) ( r , r (cid:48) ) = ν · ( r − r (cid:48) ) − ( ν · r − ν z z (cid:48) )(1 − cos( θ − θ (cid:48) )) , (89) d θ × ν (cid:48) ( r , r (cid:48) ) = τ (cid:48) · ( r − r (cid:48) ) + ( τ (cid:48) · r (cid:48) + ν (cid:48) ρ z )(1 − cos( θ − θ (cid:48) )) , (90) d θ × τ (cid:48) ( r , r (cid:48) ) = − ν (cid:48) · ( r − r (cid:48) ) − ( ν (cid:48) · r (cid:48) − ν (cid:48) z z )(1 − cos( θ − θ (cid:48) )) , (91) d θ × θ (cid:48) ( r , r (cid:48) ) = − ( z − z (cid:48) ) sin( θ − θ (cid:48) ) . (92) B Layer potentials for field evaluations
The expressions for field-evaluation (27-30) involve layer potentials ˜ S αk and˜ K αk which are defined analogously to the operators S αk and K αk in Ap-pendix A. The only difference is that ˜ S αk g ( r ) and ˜ K αk g ( r ) have r ∈ Ω − ∪ Ω + while S αk g ( r ) and K αk g ( r ) have r ∈ Γ. Therefore ˜ K αk does not need theprincipal value.When coding, it is helpful to have ˜ S αk and ˜ K αk in the form˜ S αk g ( r ) = i k (cid:90) Γ s α ( r , r (cid:48) )2 π | r − r (cid:48) | e i k | r − r (cid:48) | g ( r (cid:48) ) dΓ (cid:48) , r ∈ Ω − ∪ Ω + , (93)˜ K αk g ( r ) = (cid:90) Γ d α ( r , r (cid:48) )2 π | r − r (cid:48) | (1 − i k | r − r (cid:48) | ) e i k | r − r (cid:48) | g ( r (cid:48) ) dΓ (cid:48) , r ∈ Ω + ∪ Ω − , (94)where s α ( r , r (cid:48) ) and d α ( r , r (cid:48) ) are static kernel factors, some of which arealready listed Appendix A. The static kernel factors needed, not listed inAppendix A, are s ρ · ν (cid:48) ( r , r (cid:48) ) = ν (cid:48) ρ cos( θ − θ (cid:48) ) , (95) s ρ · τ (cid:48) ( r , r (cid:48) ) = ν (cid:48) z cos( θ − θ (cid:48) ) , (96) s ρ · θ (cid:48) ( r , r (cid:48) ) = sin( θ − θ (cid:48) ) , (97) s z · ν (cid:48) ( r , r (cid:48) ) = ν (cid:48) z , (98) s z · τ (cid:48) ( r , r (cid:48) ) = − ν (cid:48) ρ , (99)30nd d ρ ( r , r (cid:48) ) = ρ − ρ (cid:48) cos( θ − θ (cid:48) ) , (100) d z ( r , r (cid:48) ) = z − z (cid:48) , (101) d ρ × ν (cid:48) ( r , r (cid:48) ) = − ( τ (cid:48) · r (cid:48) + ν (cid:48) ρ z ) sin( θ − θ (cid:48) ) , (102) d ρ × τ (cid:48) ( r , r (cid:48) ) = ( ν (cid:48) · r (cid:48) − ν (cid:48) z z ) sin( θ − θ (cid:48) ) , (103) d ρ × θ (cid:48) ( r , r (cid:48) ) = ( z − z (cid:48) ) cos( θ − θ (cid:48) ) , (104) d z × ν (cid:48) ( r , r (cid:48) ) = ν (cid:48) ρ ρ sin( θ − θ (cid:48) ) , (105) d z × τ (cid:48) ( r , r (cid:48) ) = ν (cid:48) z ρ sin( θ − θ (cid:48) ) , (106) d z × θ (cid:48) ( r , r (cid:48) ) = ρ (cid:48) − ρ cos( θ − θ (cid:48) ) . (107) References [1] L. af Klinteberg, F. Fryklund, and A.-K. Tornberg. An adaptivekernel-split quadrature method for parameter-dependent layer potentials. arXiv:1906.07713 [math.NA] , 2019.[2] H. Ammari, Y. Deng, and P. Millien. Surface plasmon resonance of nanoparti-cles and applications in imaging.
Arch. Ration. Mech. Anal. , 220(1):109–153,2016.[3] K. Ando, H. Kang, Y. Miyanishi, and M. Putinar. Spectral analysis ofNeumann–Poincar´e operator. arXiv:2003.14387 [math.SP] , 2020.[4] A.-S. Bonnet-Ben Dhia, C. Carvalho, L. Chesnel, and P. Ciarlet. On the use ofperfectly matched layers at corners for scattering problems with sign-changingcoefficients.
J. Comput. Phys. , 322:224–247, 2016.[5] D. Colton and R. Kress.
Integral equation methods in scattering theory , vol-ume 72 of
Classics in Applied Mathematics . Society for Industrial and AppliedMathematics (SIAM), Philadelphia, PA, 2013. Reprint of the 1983 original[MR0700400].[6] C. L. Epstein, L. Greengard, and M. O’Neil. Debye sources and the numericalsolution of the time harmonic Maxwell equations II.
Comm. Pure Appl. Math. ,66(5):753–789, 2013.[7] C. L. Epstein, L. Greengard, and M. O’Neil. A high-order wideband directsolver for electromagnetic scattering from bodies of revolution.
J. Comput.Phys. , 387:205–229, 2019.[8] M. Ganesh, S. Hawkins, C. Jeznach, and D. Volkov. A well-posed surface cur-rents and charges system for electromagnetism in dielectric media.
J. IntegralEqu. Appl. , 32(1):1–18, 2020.[9] J. Helsing. Solving integral equations on piecewise smooth boundaries usingthe RCIP method: a tutorial. arXiv:1207.6737v9 [physics.comp-ph] , revised2018.
10] J. Helsing, H. Kang, and M. Lim. Classification of spectra of the Neumann-Poincar´e operator on planar domains with corners by resonance.
Ann. Inst.H. Poincar´e Anal. Non Lin´eaire , 34(4):991–1011, 2017.[11] J. Helsing and A. Karlsson. An explicit kernel-split panel-based Nystr¨omscheme for integral equations on axially symmetric surfaces.
J. Comput. Phys. ,272:686–703, 2014.[12] J. Helsing and A. Karlsson. An extended charge-current formulation of theelectromagnetic transmission problem.
SIAM J. Appl. Math. , 80(2):951–976,2020.[13] J. Helsing and K.-M. Perfekt. On the polarizability and capacitance of thecube.
Appl. Comput. Harmon. Anal. , 34(3):445–468, 2013.[14] J. Helsing and K.-M. Perfekt. The spectra of harmonic layer potential oper-ators on domains with rotationally symmetric conical points.
J. Math. PuresAppl. , 118:235–287, 2018.[15] J. Helsing and A. Ros´en. Dirac integral equations for dielectric and plasmonicscattering. arXiv:1911.00788v2 [math.AP] , revised 2020.[16] R. Kleinman and P. Martin. On single integral equations for the transmissionproblem of acoustics.
SIAM J. Appl. Math. , 48(2):307–325, 1988.[17] R. Kress. On the spectrum of the magnetostatic integral operator. In
Inversescattering and potential problems in mathematical physics (Oberwolfach, 1993) ,volume 40 of
Methoden Verfahren Math. Phys. , pages 95–105. Peter Lang,Frankfurt am Main, 1995.[18] J. Lai and M. O’Neil. An FFT-accelerated direct solver for electromagneticscattering from penetrable axisymmetric objects.
J. Comput. Phys. , 390:152–174, 2019.[19] J. Li, X. Fu, and B. Shanker. Decoupled potential integral equations for elec-tromagnetic scattering from dielectric objects.
IEEE Trans. Antennas Propag. ,67(3):1729–1739, 2019.[20] C. M¨uller.
Foundations of the mathematical theory of electromagnetic waves .Revised and enlarged translation from the German. Die Grundlehren der math-ematischen Wissenschaften, Band 155. Springer-Verlag, Berlin, 1969.[21] K.-M. Perfekt. Plasmonic eigenvalue problem for corners: Limiting absorptionprinciple and absolute continuity in the essential spectrum.
J. Math. PuresAppl. , (in press 2020).[22] A. Ros´en.
Geometric multivector analysis . Birkh¨auser AdvancedTexts: Basler Lehrb¨ucher. [Birkh¨auser Advanced Texts: Basel Textbooks].Birkh¨auser/Springer, Cham, [2019] c (cid:13)
Comm.Part. Differ. Equat. , 43(2):159–184, 2018.[24] P. Young, S. Hao, and P. G. Martinsson. A high-order Nystr¨om discretizationscheme for boundary integral equations defined on rotationally symmetric sur-faces.
J. Comput. Phys. , 231(11):4142–4159, 2012., 231(11):4142–4159, 2012.