Comparison of integral equations for the Maxwell transmission problem with general permittivities

Two recently derived integral equations for the Maxwell transmission problem are compared through numerical tests on simply connected axially symmetric domains for non-magnetic materials. The winning integral equation turns out to be entirely free from false eigenwavenumbers for any passive materials, also for purely negative permittivity ratios and in the static limit, as well as free from false essential spectrum on non-smooth surfaces. It also appears to be numerically competitive to all other available integral equation reformulations of the Maxwell transmission problem, despite using eight scalar surface densities.


Introduction
The Maxwell transmission problem is about determining the fields that result when an incident time-harmonic electromagnetic wave is scattered from and transmitted into a bounded dielectric object. Two pioneering integral equation reformulations (IERs) of this problem, which are still popular, are the Müller IER [28,Section 23] and the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) IER [27,34]. These IERs are of Fredholm's second and first kind, respectively. Over the years, much effort has been made to find more efficient IERs. Focus has been on avoiding dense-mesh and topological low-frequency breakdown, on avoiding false resonances, and on providing unique solutions for wider ranges of material parameters. Among later contributions we mention [10,12,13,19,22,25,26,33].
In this work we compare the numerical performance of two recently derived IERs of the Maxwell transmission problem. The two IERs, which are of Fredholm's second kind with singular integral operators, are referred to as "Dirac" and "HK 8-dens". "Dirac" is derived in [22] by embedding the timeharmonic Maxwell's equations into a Dirac equation and by tuning the six free parameters of this equation as to optimize numerical performance. "HK 8-dens" is derived in [19] by extending a classic IER of the Helmholtz transmission problem [23] via the use of four uniqueness parameters. Both our IERs use eight unknown scalar surface densities for modeling. This is more than other popular IERs use. Typical numbers are four [10,25,27,28,34] or six [12,26,33]. The major advantage with our new IERs, however, is that they offer unique solutions, that is they are free from false eigenwavenumbers, in particular for plasmonic scattering. Another advantage, from a programming point of view, is that both IERs require only bounded integral operators with double and single layer type kernels, and in particular do not use (compact differences of) hypersingular integral operators.
Uniqueness for a wide range of parameters is an important property of an IER of a parameter-dependent partial differential equation (PDE) for many reasons. First, one may actually be interested in solving the PDE for a wide range of parameters. Then uniqueness is obviously important. Second, even if one is not interested in a wide range of parameters, non-uniqueness outside the parameter regime of interest can seriously affect the conditioning of an IER inside the parameter regime of interest. Third, theoretical studies of the solvability of a PDE, for example the time-harmonic Maxwell's equations in Lipschitz domains, are often based on IERs. Then it is crucial that the IER has the same solvability properties as the PDE it models.
While "Dirac" and "HK 8-dens" are similar in many respects, there are also differences. "Dirac" is derived in [22] for general magnetic materials and assuming only Lipschitz regularity of the surface of the scattering object. "HK 8-dens" has the advantage that its surface densities have immediate interpretations in terms of boundary limits of physical fields. Two of these surface densities are always zero and are only needed, for uniqueness, in the main linear system that is to be solved. They are omitted in field evaluations, whose cost then are the same as for six-density IERs.
The original papers [19,22] use mutually different notation and contain numerical 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 unified and programming-friendly notation and to conduct a series of numerical experiments for their full versions in three dimensions as to see which IER is the most efficient. For this comparison, we limit the discussion to nonmagnetic materials, for which "HK 8-dens" is formulated. Moreover, the free parameters in "Dirac" are likely to need further tuning in the magnetic case, which will appear in a forthcoming publication.
We conclude the paper in Section 11 by comparing some salient properties of our IERs to those of two available competitive IERs for the Maxwell transmission problem: the "Debye" and the "DFIE" IERs [9,33]. With the aim of achieving a comparison as accurate as possible, given the available

Problem formulation
Let Ω + be a bounded domain in R 3 with boundary surface Γ and an unbounded, connected, exterior Ω − . The outward unit normal at position r on Γ is ν. We consider time-harmonic fields with time dependence e −iωt and angular frequency ω > 0. The relation between time-dependent fields F (r, t) and complex fields F (r) is The domains Ω ± are homogeneous with material properties described by wavenumbers k ± . See Figure 1(a). The k ± are related to the total permittivities ± and permeabilities µ ± by k ± = ω √ ± µ ± . To avoid issues regarding the choice of branch of the square root, the k ± are considered our basic parameters. Passive materials have m{ ± } ≥ 0 and m{µ ± } ≥ 0 and values of k ± in the closed first quadrant. In Maxwell's equations, ∇ × E = iωµ ± H and ∇ × H = −iω ± E in Ω ± , we rescale the magnetic field H by the wave impedance in Ω − and set E = E and H = µ − / − H in all R 3 . Throughout the paper, only non-magnetic materials are considered. That is, µ ± = µ 0 , where µ 0 is the permeability of vacuum. We now state our formulation of the Maxwell transmission problem: Given incident fields E in and H in , generated in Ω − , we seek the total electric 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ˆ such thatˆ solve Maxwell's equations except possibly at an isolated point in Ω − where the source of E in and H in is located, subject to the boundary conditions lim lim The scattered fields E sc and H sc are source free in Ω − and defined, along with the transmitted fields E tr and H tr , by The incident fields satisfy except at the possible isolated source point in Ω − . In addition, conservation of charge must hold and by that where E sc− is the exterior limit of E sc on Γ. In what follows, the Maxwell transmission problem (3)-(10) will be referred to as the MTP(k − , k + ).
3 Properties of IERs of the MTP(k − , k + ) Let us first mention that there is yet no IER of the MTP(k − , k + ) known that 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 for which (k − , k + ) there exist unique solutions to the MTP(k − , k + ). The most comprehensive 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) [22,Proposition 8.2]. As for existence of solutions, the same set of points apply with the restriction thatˆ = −1 for smooth Γ and thatˆ is outside of a geometry-dependent interval on the negative real axis, containingˆ = −1, for merely Lipschitz regular Γ [22,Proposition 8.3].
For "Dirac" it is proven that there exist unique solutions when Ω + is a bounded Lipschitz domain with connected exterior Ω − and (k − , k + ) satisfies the conditions of [22,Proposition 8.3]. Existence of unique solutions for "HK 8-dens" is proven for (k − , k + ) as described in [19,Section 8.3]. It is assumed in [19] that Γ is smooth and Ω + simply connected but the proof holds also for objects that are multiply connected. All other IERs of the MTP(k − , k + ) in the electromagnetics literature, possibly with the exception of the Debye representation [9], seem to have unique solutions only under more restrictive conditions. In particular, they can not guarantee uniqueness for any (arg(k − ), arg(k + )) = (0, π/2).
Two problems shared by many IERs of the MTP(k − , k + ), which have received much attention recently, are dense-mesh and topological low-frequency breakdown, see [33] and [10,Section 2]. Dense-mesh low-frequency breakdown refers to catastrophic cancellation that destroys the numerical accuracy in the computed fields in the static limit. Topological low-frequency breakdown is a, perhaps, more elusive phenomenon that broadly seems to refer to an increased ill-conditioning of the integral equation itself in the quasi-static limit k ± → 0, withˆ = (k + /k − ) 2 fixed, for Γ with non-zero genus. For scattering by superconductors, this phenomenon is discussed in [7,8]. "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 lowfrequency breakdown but still cause numerical degradation of IERs, include false near-eigenwavenumbers and false essential spectrum. An eigenwavenumber is a pair of wavenumbers (k − , k + ) for which the IER does not have a unique solution. If the MTP(k − , k + ) has a unique solution we speak of a false eigenwavenumber, whereas we call it a true eigenwavenumber if even the MTP(k − , k + ) fails to have a unique solution. If we only consider wavenumber pairs from a set X and there is an eigenwavenumber z outside but close to X, then a pair x ∈ X near z which locally maximizes the condi-(c) Figure 2: Notation for an axially symmetric surface Γ generated by a curve γ: (a) a point r on Γ has outward unit normal ν and unit tangent vectors τ and θ; (b) cylindrical coordinates (ρ, θ, z) of r; (c) coordinate axes and vectors in the half-plane A defined by θ = 0. tion number of the IER, is referred to a near-eigenwavenumber. This can be a true or false such, depending on the nature of z. The terminology of true and false eigenwavenumbers was introduced in [17, Section 5.3], to make more precise the common terminology of true and spurious resonances. The term spurious near resonances is used in [33].
False essential spectrum may appear for certain (k f − , k f + ) withˆ f = (k f + /k f − ) 2 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 essential spectrum if it is not a Fredholm operator, even though the MTP(k f − , k f + ) defines a Fredholm map. While the term essential spectrum is standard, the term false/spurious essential spectrum does not seem to have been used.
Let k − = k f − and k + → k f + in such a way thatˆ →ˆ f from above or from below in the complex plane. At a point (k f − , k f + ) where we have false essential spectrum, the typical numerical behavior of the IER is that it has the same unique limit solution from above and below and this solution coincides with the solution to the MTP(k f − , k f + ), while the IER is not solvable for (k − , k + ) = (k f − , k f + ).

Notation for axially symmetric Γ
From now on we assume that Γ is axially symmetric, since our numerical examples cover such domains. Note, however, that the full 3D formulations of our IERs and of the representations of E and H are given in [19,22]. These IERs are in no way limited to axially symmetric domains. In the present paper we use both Cartesian coordinates x, y, z and cylindrical coordinates ρ, θ, z, as illustrated in Figure 2. A half-plane A in R 2 is defined by θ = 0, a generating curve γ is defined by the intersection of Γ and A, and a general point in A is The outward unit normal and a tangent on γ, in A, are In R 3 , the position with Cartesian basis vectors is r ≡ (x, y, z) = (ρ cos θ, ρ sin θ, z) .
We also use the cylindrical unit vectors On Γ, the unit normal is ν and the tangential unit vectors are τ and θ ν = (ν ρ cos θ, ν ρ sin θ, ν z ) , Vector fields, off from Γ, will often be expressed using ρ, θ, and z E = ρE ρ + θE θ + zE z , The causal fundamental solution to the Helmholtz' equation and its gradient, renormalized by a factor of −2 to make the expressions of the integral operators in Appendix A simpler, are

A unified formalism
The integral equations of "Dirac" and "HK 8-dens" can both be written in the general form with Here h contains eight unknown scalar surface densities f in contains the field components E k is the Cauchy singular 8 × 8 block operator matrix with operator entries detailed in Appendix A, and P , P , N , N are diagonal matrices which are specified by "Dirac" and "HK 8-dens", respectively.
Once (23) is solved for h, the fields E and H in Ω ± , decomposed as in (20), can be evaluated from (8) and Here the layer-potential entries are detailed in Appendix B and Note that the additional zero entry of E z and H z is due to axial symmetry. For "Dirac", the field evaluation can be improved by instead using the projected densities as input to the field evaluation (28)- (31). The two choices (33) and (34) will always produce the same fields in Ω ± , and we discuss their uses in the following subsections. We remark that the framework (23)-(27) applies to both "Dirac" and to "HK 8-dens", also for general object shapes Γ, whenever τ and θ are two arbitrary orthogonal tangential unit vectors. The corresponding representations of E and H are then given by [22,Eqs. (22) and (23)].
5.1 Choice of P , P , N , N for "Dirac" Given the formalism above, "Dirac" specifies the diagonal matrices whereĉ = 1/k andˆ =k 2 as in (2). The matrices P and N differ in their first elements from those used in [22,Theorem 2.3], and we have set µ = 1 as we consider non-magnetic materials. The modified first elements correspond to another choice of the parameter β in [22, Section 8], using here β = 1 + iδ arg(k) instead of β = 1, to avoid false eigenwavenumbers when (arg(k − ), arg(k + )) = (π/2, 0). We use Computations show that this is close enough to δ = 0 not to affect speed and accuracy, but large enough to eliminate false eigenwavenumbers and neareigenwavenumbers. In the original formulation of "Dirac" in R 3 from [22], false eigenwavenumbers appear when (arg(k − ), arg(k + )) = (π/2, 0), that is, when the wavenumber in the exterior region is imaginary. This corresponds to the circled lower corner point in Figure 1(b) and is confirmed by numerical experiments on the unit sphere and comparison with semi-analytic results given by Mie theory. For the parameters (35), with δ = 0.2/π, it follows from [22,Proposition 8.5] that there are no false eigenwavenumbers for "Dirac" in the shaded region in Figure 1(b). Not even at (arg(k − ), arg(k + )) = (π/2, 0), where there are true eigenwavenumbers. Note that "Dirac" is not a Fredholm second kind integral equation with compact operators on smooth Γ. For one thing, the block operator G in (24) contains Cauchy singular differences of operators. However, the particular choice of P , P , N , N in the original "Dirac", that is for δ = 0, makes G 4 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 when using iterative solvers for (23).
When evaluating E and H with "Dirac", one can use either (33) or (34) for h ± . The reason for preferring (34), which we use for smooth Γ in the numerical examples of Section 10, is that, like in (26), components 1 and 5 of h ± from (34) are zero. This leads to at most five non-zero densities in the evaluation of each field. However, for non-smooth Γ the numerical method used in Section 10 is less compatible with (34) and we use the simpler (33). For "Dirac", the densities computed with (34) satisfy so only one of the Cauchy integrals in (34) needs to be computed.
5.2 Choice of P , P , N , N for "HK 8-dens" "HK 8-dens" specifies the diagonal matrices where γ 1 , γ 2 , η, and λ are uniqueness parameters whose determination is discussed in [19, Section 11.1 and Appendix D]. Note that in [19], these parameters are called γ E , γ M , c, and λ. A valid choice when (arg(k − ), arg(k + )) = (0, 0) is 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 + )) = (π/2, 0). For this reason (39) is used also at (arg(k − ), arg(k + )) = (π/2, 0). A valid choice when (arg(k − ), arg(k + )) = (0, π/2) is The surface densities h of (25) have, with "HK 8-dens", the physical interpretations where −ik − σ E and −ik − σ M are exterior limits of the electric and magnetic volume charge densities on Γ, E and M are the equivalent electric and magnetic surface charge densities on the exterior side of Γ, and J θ , J τ , M θ , M τ are components of the equivalent electric and magnetic surface current densities on the exterior side of Γ. See [19,Remark 10.2], where it also is shown that h 1 = σ E and h 5 = σ M must be zero on theoretical grounds. Therefore, the preferred choice for field evaluations with "HK 8-dens" is always via (33), which is the Stratton-Chu representation, and (34) is never needed since it does not lead to any reduced numerical costs.
We remark that "HK 8-dens" was discovered independently and prior to "Dirac", and that it was only later realized that it can be written in the form (23)-(24), just like "Dirac". However, it is not the case that "HK 8dens" is a special case of "Dirac", corresponding to a certain choice of Dirac parameters as in [22,Section 8]. Indeed, "HK 8-dens" is not derivable from jump matrices M and M as in [

False essential spectra
A key operator for the MTP(k − , k + ) is the static, k → 0, limit of the acoustic double layer operator K ν k , appearing in the (1, 1) and (5, 5) diagonal blocks of E k . That is, K d equals the Neumann-Poincaré operator K NP , possibly modulo a sign depending on convention [3]. Its essential spectrum σ ess (K d ) in the fractional Sobolev space H 1/2 (Γ), that is the set of λ for which λI − K d fails to be a Fredholm operator, is a compact subset of the interval (−1, 1), for any Lipschitz surface Γ. Unlike in R 2 , σ ess (K d ) is not necessarily symmetric with respect to 0 for Γ ⊂ R 3 : see examples for the spectrum of the adjoint of K NP , σ ess (K * NP ) in H −1/2 (Γ), and Γ with axially symmetric conical points in [21,Section 7.3].
We map σ ess (K d ) onto the negative real axis and define For E k and our IERs, the relevant function space is  (26). Furthermore, as we shall prove in a forthcoming publication, the MTP(k − , k + ) defines a Fredholm map if and only ifˆ / ∈ Σ(Γ), for any wavenumbers and not only in the quasistatic limit.
Another key operator for the MTP(k − , k + ) is the magnetic dipole operator, acting on tangential vector fields g. We note that the static limit of the operator appearing in the (3:4,3:4) and (7:8,7:8) size 2 × 2 diagonal blocks in E k , is −K * m . Moreover, the static limit of the normal derivative of the acoustic single layer potential, K ν k , appearing in the (2, 2) and (6, 6) diagonal elements, equals −K * d . By Hodge decomposition of H −1/2 (curl, Γ), it can be shown that the essential spectrum of K m is σ ess (K d )∪(−σ ess (K d )). The corresponding result for eigenvalues is in [2,Proposition 4.7] and in [24]. However, these results are proved for smooth Γ, on which K m is compact and σ ess ( In the diagonal blocks of the matrix E k , we therefore find operators with essential spectra σ ess (K d ) as well as −σ ess (K d ). To avoid the false essential spectrum (−σ ess (K d )) \ σ ess (K d ), the matrices P, P , N, N need to be chosen carefully. Consider the diagonal operator blocks of (27) for "Dirac" and "HK 8-dens" whenˆ is negative real andk/i is positive, referred to as the plasmonic case in Section 10. The only operator block which can fail to be a Fredholm operator for "Dirac" is (6,6). This happens whenˆ ∈ Σ(Γ), that is, if and only if the MTP(k − , k + ) itself fails to define a Fredholm map.
The spectral properties of the diagonal blocks show that "Dirac" has no false essential spectrum, but also indicate that "HK 8-dens" may have false essential spectrum. Section 10.4 contains numerical results that support this. However, since the space H 3 is a space of mixed ±1/2 regularity, a full proof must include an analysis of the off-diagonal blocks. We plan a forthcoming publication devoted to a more careful theoretical and numerical study of the essential spectrum and the issues discussed in this section.

The quasi-static limit
In the quasi-static limit k ± → 0, withk = k + /k − fixed, the operator G from (24) simplifies considerably. The diagonal matrices P, P , N, N are all fixed whereas our basic Cauchy integral operator becomes a 4/4 block diagonal operator 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 (27) but with k = 0. Now all single layer operator entries vanish, as they contain a factor k ± , and the equations for the magnetic and electric fields, which correspond to the two diagonal blocks, decouple. For spectral properties of the operators K d and K m defined in Section 6, used in this section, we refer to [6, Sections 5.1-5.2]. These results generalize to Lipschitz surfaces, if we use the spaces H 1/2 (Γ) and H −1/2 (curl, Γ) for K d and K m respectively. Let MTP(0, 0,k) denote the Maxwell transmission problem in the quasistatic limit withk fixed. The MTP(0, 0,k) amounts to two decoupled divergence and curl free vector fields E and H in Ω ± , having continuous tangential parts (4)- (5) and with E sc and H sc decaying at infinity. We now also explicitly need to require the Gauss jump conditions lim We note that since we consider non-magnetic materials, the jump condition for all components of H are the same. It is only the (5:8,5:8) diagonal block in (46) which needs to be inverted in the quasi-static limit. Moreover, one can show that the MTP(0, 0,k) defines an invertible map if and only if (1 +ˆ )/(1 −ˆ ) / ∈ σ(K d ). In Section 10.6 we show numerical results for condition numbers of "Dirac" as well as of "HK 8-dens" in this quasi-static limit, but restrict our analysis here to "Dirac". For simplicity, assume δ = 0 in the definition of P, P , N, N in (35). This makes the blocks in G corresponding to K 1,3:4 , K 2,3:4 , K 7:8,5 and K 7:8,6 vanish, and invertibility of I + G is determined by the diagonal blocks. (δ ≈ 0 makes K 1,3:4 ≈ 0.) The main operator is the (6,6) diagonal block, which fails to be invertible for "Dirac" precisely when (1 +ˆ )/(1 −ˆ ) ∈ σ(K d ). Since the spectra of K d , K * d and K * m are subsets of [−1, 1], for any topology of Γ, the remaining diagonal operators are all seen to be invertible for anyk = 0 andk = ∞.
Preconditioning similarly to [22,Theorem 2.3], we arrive at where β = 1 − iδ arg(ĉ) is chosen as in Section 5.1. When k ± ≈ 0 andk ≈ 0, the version (51) of "Dirac" is better conditioned than that of (35). Indeed, one verifies for these choices of parameters and P, P that I + G converges in operator norm as k ± → 0 andk → 0, to a limit operator which is invertible due to the block diagonal and triangular structures and the fact that the (1,1) block is now invertible as above.
For the MTP(0, 0, ∞), "Dirac" again has a false eigenwavenumber, this time due the (6,6) block. It is not clear how to avoid this by tuning the free Dirac parameters. We plan to address this issue, as well as Ω + that are not simply connected, in a forthcoming paper.

Surface plasmon standing waves and plasmons
Let, for the moment, R 3 be divided into two, not necessarily bounded, domains Ω 1 and Ω 2 separated by a surface Γ and with real-valued permittivities 1 and 2 such that 1 > 0 and 2 < 0. Surface plasmon waves, surface plasmon standing waves, and quasi-static plasmons are then particular types of electromagnetic fields that may appear as solutions to the MTP(k − , k + ). These fields are briefly described here from a classical electrodynamics point of view and they are referenced in the discussion of the numerical examples in Section 10.
Surface plasmon waves (SPWs) are surface waves that can propagate along Γ. A necessary extra condition for their existence on planar Γ iŝ Here λ 1 = 2π/k 1 is the free space wavelength in Ω 1 . Furthermore, the following has been verified by us for a circular cylinder, using semi-analytic methods, and by numerical examples in [17,19]: The SPWs can propagate without attenuation along surfaces that are invariant only in the propagation direction. They can also propagate along surfaces which are curved in the direction of propagation, but are then attenuated due to radiation. The radiation increases with the ratio of λ spw to the radius of curvature. A natural conjecture is that SPWs appear, excited at corners or edges on Γ, when (1+ˆ )/(1−ˆ ) hits the essential spectrum of K d . This was discussed, for edges, in [17,Section 7.3], but that discussion is not exhaustive. Indeed, as we shall see in Section 10, an SPW appears in Figure 6 even though we are outside the essential spectrum, and in Figure 8 we are in the essential spectrum, but no SPW appears. The description of the precise mechanisms for the excitation of SPWs on non-smooth Γ remains an open problem. Consider now, with notation as in Section 2, a bounded object Ω + and assume that the wavenumber k − is positive real and thatˆ < −1, to enable SPWs along Γ. At certain k − the SPWs form standing waves. Such SPWs are here called surface plasmon standing waves (SPSWs). In some literature the name surface plasmon resonances (SPRs) is used to emphasize their resonant nature. However, SPR is not a direct synonym for SPSW since it is also used for other plasmonic phenomena. Due to radiation, the SPSW is a damped resonant field and the (k − , k + ), for which it appears, is close to a true eigenwavenumber outside the set of points shown in Figure 1(b). That is, SPSWs appear close to true near-eigenwavenumbers, in the terminology from Section 3. When an SPSW is excited by an incident field, its amplitude becomes large, which results in a large scattering cross section. This is one reason why SPSWs are of interest in optics.
A quasi-static plasmon is an electromagnetic field that approaches an eigenfield, referred to as a static plasmon, as k ± → 0. It appears around objects that are much smaller than the wavelength λ − . The quasi-static plasmon is also denoted surface plasmon (SP) in optics, a name that may lead to confusion since it is sometimes used as a synonym for SPSW. A quasi-static plasmon is a resonant field, but in contrast to SPSWs it is not a standing wave field. For smooth Γ there is an infinite discrete set ofˆ corresponding to static plasmons. The electric field of a static plasmon is distributed so that the electric energies in Ω − and Ω + exactly cancel each other, while the magnetic field is zero. Indeed, Green's identity shows thatˆ Ω + |E| 2 dx = − Ω − |E| 2 dx so that, in accordance with Poincaré's variational principle [3, Theorem 2.4], the static plasmon corresponds to the eigenvalue (1 +ˆ )/(1 −ˆ ) 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 or not. A bright plasmon radiates as an electric dipole and its far-field is very large considering the object is small compared to λ − . A dark plasmon radiates as a higher-order multipole and its far-field is weak. See [2, Lemma 5.3], where the full multipole expansion of far-fields is given. Furthermore, there is a close connection between bright static plasmons and the imaginary part of an object's limit polarizability [20,29]. For an object with sharp corners, edges, or points there is a continuous, possibly punctured, interval ofˆ where a special type of static plasmons, with infinite (but cancelling) electric energies in Ω − and Ω + , can occur together with a partially embedded discrete set of regular static plasmons [15,21,29].

Scattering objects and discretization
This section reviews shapes of scattering objects and discretization schemes that are used for the numerical tests in Section 10.

Two surface families
Examples with smooth Γ come from a generating curve γ parameterized as r(s) = (1 + α sin(5s))(cos(s), sin(s)) , s ∈ [−π/2, π/2] , where α is a shape parameter. The choice α = 0 corresponds to Γ being the unit sphere. With α = 0.25 the shape of Γ resembles a "rotated starfish". See Figure 2 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 sharp tip". The generating curve γ of (54) has previously been used for numerical examples in [19,21].

Fourier-Nyström discretization
The integral equation (23) is solved on axially symmetric Γ using highorder Fourier-Nyström discretization. An azimuthal Fourier transform is first applied to (23), yielding a sequence of modal problems for the Fourier coefficients h (n) of h which are solved independently using high-order panel-based Nyström discretization. The modal solutions h (n) are then synthesized to give the full solution h. High-order panel-based Fourier-Nyström discretization for solving integral equations modeling scattering problems on axisymmetric surfaces was made popular by Young, Hao, and Martinsson in 2012 [35]. Since then, several authors have worked on extensions of such schemes and related issues. 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 imaginary parts [10,25,16,1]. Our version of the Fourier-Nyström scheme is the one used in [19]. This version is of 16th order. It relies on a combination of 16-point and 32-point underlying Gauss-Legendre quadrature, a variety of explicit kernel-splits, semi-analytic product integration computed on the fly, and is compatible with the recursively compressed inverse preconditioning method (RCIP) [14]. The RCIP accelerates and stabilizes Nyström discretization in the presense of singular boundary points on γ, such as corners.

Numerical examples
The numerical efficiency of "Dirac" and "HK 8-dens" is now compared. That is, we solve the discretized modal systems (55) and compute fields via (28)-(31) with P , P , N , N as in Sections 5.1 and 5.2. We also compute condition numbers of the modal systems. Three material parameter cases are used: • the positive dielectric case, wherek = 1.5 and k − is positive real; • the plasmonic case, wherek = i √ 1.1838 and k − is positive real; • the reverse plasmonic case, wherek = i √ 1.1838 −1 and k + is positive real.
These parameter cases are taken from previous work on time-harmonic transmission problems [4,19,22].
In all examples involving field images, unless otherwise stated, the incident field is a linearly polarized plane wave E in (r) = xe ik − z with x = (1, 0, 0). The fields are plotted in the plane y = 0, where the field components E θ , H ρ , and H z are zero due to symmetry. To save space in the examples, we only show E ρ and H θ . Generally, these two seem to exhibit more pronounced field patterns than the omitted component E z .
When Γ is non-smooth, it may happen thatk = i √ 1.1838 ork = i √ 1.1838 −1 correspond to that the MTP(k − , k + ) does not define a Fredholm map. We then compute limit solutions asˆ approaches −1.1838 or −1/1.1838 from above in the complex plane. Such limit solutions have boundary traces lying outside the function space H 3 from (44) and are given a downarrow superscript. For example, the limit of the field E is denoted E ↓ . Our codes are implemented in Matlab, release 2018b, and executed on a workstation equipped with an Intel Core i7-3930K CPU and 64 GB of RAM. Large linear systems are solved iteratively using GMRES with a stopping criterion threshold of machine epsilon ( mach ) in the estimated relative residual. When assessing the accuracy of computed field quantities we adopt a procedure where to each numerical solution we also compute an overresolved reference solution, using roughly 50% more points in the discretization of the system under study. The absolute difference between these two solutions is denoted the estimated absolute error. In examples involving field images, the fields are computed at 10 6 target points on a rectangular Cartesian grid in the computational domains shown.
Subsections 10.1, 10.3, and 10.6, on eigenwavenumbers, convergence, and the quasi-static limit, involve the unit sphere. An advantage with doing experiments on spheres is that semi-analytic results given by Mie theory can be used for verification. The remaining subsections concern less trivial object shapes and there we limit ourselves mainly to the plasmonic case. The positive dielectric case is not deemed challenging enough for thorough testing. The performances of "Dirac" and "HK 8-dens" in the reverse plasmonic case is rather similar to their respective performances in the plasmonic case -apart from the new risk of being close to a true eigenwavenumber and, for "HK 8-dens", also close to a false eigenwavenumber or near-eigenwavenumber.

Condition numbers on the unit sphere
We compute condition numbers of system matrices of the discretized modal integral equations (55) and with Γ being the unit sphere. Results for the azimuthal modes n = 0, 5, 10 are shown in Figure 3. A number of 768 discretization points are placed on a grid on γ, making the total system size 6144×6144, and up to 5,300 data points are used to capture the rapid variations in the condition numbers as k − and k + are swept through the intervals [0, 10] and [0, 6]. The purpose of this study is primarily to detect possible false eigenwavenumbers, but also to get a notion of how well-conditioned the two IERs under study are. Figure 3 shows that "Dirac" and "HK 8-dens" have similar condition numbers in the positive dielectric case, while "Dirac" is better conditioned than "HK 8-dens" in the plasmonic-and reverse plasmonic cases -particularly at higher wavenumbers k − and k + . Note that the regularly recurring high peaks that are common to Figures 3(c,d) correspond to true near-eigenwavenumbers, each with multiplicity one for a given n, while the additional peaks in Figure 3(d) correspond to false near-eigenwavenumbers of "HK 8-dens", as discussed in the last paragraph of Section 3. The true near-eigenwavenumbers are associated with SPSWs. The first peak in Figures 3(c,d) corresponds to an SPSW with = 6, where is the degree of the spherical harmonic Y n , and that peak is common to the 13 n-values −6 ≤ n ≤ 6. Note also that the 14 eigenwavenumbers with k + ∈ [0, 6], visi- ble in Figure 3(e), are true eigenwavenumbers in the reverse plasmonic case, as confirmed to at least 13 digits by comparison with semi-analytic results. In contrast, "HK 8-dens" exhibits here, in addition, around 25 false eigenwavenumbers and false near-eigenwavenumbers. See Figure 3(f). Eleven of the 14 peaks in Figure 3(e) form a periodic pattern and correspond to eigenwavenumbers with eigenfields that are SPSWs, whereas the other three peaks correspond to eigenwavenumbers with eigenfields that are not bound to the surface.

Field images for the "rotated starfish"
We compute images of scattered and transmitted fields E sc ρ (r, 0), H sc θ (r, 0) and E tr ρ (r, 0), H tr θ (r, 0). The surface Γ is that of the "rotated starfish", given by (53) with α = 0.25. A number of 768 discretization points are again used for each integral operator on γ in (55). We first test the positive dielectric case (k = 1.5) with k − = 10. This is a simple problem and both "Dirac" and "HK 8-dens" give fields with an estimated absolute precision of at least twelve digits -also close to Γ. We refrain from showing images. The only noticeable difference between the methods is the number of GMRES iterations needed for full convergence. The "Dirac" system converges to an estimated relative residual of mach in 111 iterations while the "HK 8-dens" system is a bit better and only needs 68 iterations.
We next test the plasmonic case (k = i √ 1.1838) with k − = 6. Figure 4 shows that the accuracy achieved by "Dirac" and "HK 8-dens" is similar also in this case, although the spectral properties of "Dirac" are now much better than those of "HK 8-dens". The condition number of the system matrix in the discretized "Dirac" system is 9.7 · 10 3 while the corresponding condition number for "HK 8-dens" is 1.2 · 10 5 and this is reflected in the number of GMRES iterations needed for full convergence. The "Dirac" system converges to an estimated relative residual of mach in 170 iterations while the "HK 8-dens" system needs 623 iterations.
The symmetry of E in is such that only the two azimuthal modes n = −1 and n = 1 are present. The Fourier coefficients of the surface densities of these modes, h (1) and h (−1) of (55), are either identical or have opposite signs. Animations based on (1) for a sequence of t reveal that SPWs are propagating along Γ in the images of Figure 4(a,b). Their wavelength is roughly 15% shorter than the wavelength given by (52).

Convergence of field images
As pointed out by one of the referees, it is interesting to see how the numerical error in field images, produced by our 16th order Fourier-Nyström scheme of Section 9.2 in conjunction with "Dirac", evolves as the size of the discrete linear systems (55) increases. To this end, Figure 5(a) shows the average estimated absolute field error in Figure 4(e,f) under mesh refinement. The error is estimated using 90,000 field points on a Cartesian grid in the box B = {−1.3 ≤ x ≤ 1.3; −1.2 ≤ z ≤ 1.4} and behaves very stably.
It is also of interest to make convergence studies on setups which have semi-analytic solutions, and compare with these, even though such solutions themselves may not be entirely free from numerical error. To this end we let an incident field be scattered from and transmitted into the unit ball. Here j 1 is the first order spherical Bessel function of the first kind and we are still in the plasmonic  Figure 4(e,f); (b) the scattered/transmitted E(r)-field and H(r)-field resulting from the incident field (56) and Γ being the unit sphere.
case withk = i √ 1.1838 and k − = 6. Figure 5(b) shows the average absolute error, obtained by comparison with the Mie solution, in the fields E(r) and H(r) under mesh refinement. The error is computed using 90,000 field points on a Cartesian grid in the box B = {−2 ≤ x ≤ 2; −2 ≤ z ≤ 2} and normalized with the largest field amplitudes in B. The convergence in Figure 5(b) is very similar to that in Figure 5(a).

In a false essential spectrum for the "tomato"
We repeat the experiments of Section 10.2, but now for the "tomato", that is, with Γ non-smooth and generated by γ of (54) with α = 31π/18 as illustrated in Figure 1(a). A number of 576 discretization points are used for each integral operator on γ in (55), making the total system size 4608×4608.
In the positive dielectric case, and with a larger wavenumber k − = 18 as to compensate for the "tomato" being smaller than the "rotated starfish", the results for both methods are even (marginally) better than in the previous example. The estimated pointwise precision in the field images is between twelve and 13 digits (no images shown). The "Dirac" system needs 96 GMRES iterations for full convergence while the "HK 8-dens" system needs 87 iterations. We conclude that, thanks to the RCIP method for dealing with the conical point, the non-smooth "tomato" is as simple as the smooth "rotated starfish" from a numerical point of view.
In the plasmonic case and with k − = 5 we observe some very interesting features. The wavenumber ratio corresponds to that (1 +ˆ )/(1 −ˆ ) is in −σ ess (K d ), but not in σ ess (K d ). Given the discussion about essential spectra 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", on the other hand, is free from this problem and correct fields can be computed without a limit process. Figure 6 shows that "Dirac" clearly achieves better field accuracy than "HK 8-dens" in this case. In terms of GMRES convergence the difference is even greater: "Dirac" needs 88 iterations while "HK 8-dens" needs 326 iterations. There are SPWs propagating along Γ with a wavelength that is roughly 20% shorter than the wavelength given by (52). It is also enlightening to inspect the asymptotics of h close to the conical tip of the "tomato" in the plasmonic case. Thanks to the symmetry of E in , as discussed in the last paragraph of Section 10.2, it is enough to study the eight densities (Fourier coefficients) contained in the modal solution h (1) . Figure 7 shows these eight densities both for "Dirac" and "HK 8-dens". The individual densities are denoted h i , i = 1, . . . , 8, with the azimuthal index omitted. Note that the densities h 1 and h 5 of "HK 8-dens" are approximately zero, as they should be according to Section 5.2.
The strongest singularity observed in Figure 7 is where s is the arc length distance along γ to the conical point at the origin. The exponents in (57) and (58) are estimated using the automated eigenvalue method of [14,Section 14] and all displayed digits are believed to be  correct. For "HK 8-dens", the singularity of (57) is observed for the densities h 6 , h 7 , and h 8 , corresponding to E , M θ , and M τ according to (41). For "Dirac", the singularity of (57) is observed only for the density h 6 . Note that each of these functions is more singular than H 1/2 (Γ), as allowed by the function space H 3 .

In the essential spectrum for a "drop with a sharp tip"
We repeat some of the experiments of Section 10.4, but shrink the opening angle of the conical point to α = 5π/18, so that the shape of Γ now resembles that of a "drop with a sharp tip". A number of 576 discretization points are again used for each integral operator on γ in (55).
Field images for the plasmonic case along with error estimates for "Dirac" are shown in Figure 8. No SPWs are excited along Γ. A number of 83 iterations were needed for full GMRES convergence. The corresponding number Densities displayed as in Figure 7. Column 1-2: "HK 8-dens". Column 3-4: "Dirac".
for "HK 8-dens" is 314 iterations. Note that (1 +ˆ )/(1 −ˆ ) is now in the essential spectrum of K d and that E ↓ ρ and H ↓ θ are oscillatory and unbounded at the origin. The colorbar ranges in Figure 8(a,b) are therefore restricted to the most extreme field values away from the origin. The estimated pointwise absolute error close to the origin is, of course, now larger than than in previous examples with everywhere bounded fields. Figure 9 is analogous to Figure 7 and exhibits the same general features. The strongest singularity observed is now clearly visible for h 6 of "Dirac" and for h 6 , h 7 , and h 8 of "HK 8-dens". As we are now in the essential spectrum, we do not expect the densities to belong to H 3 , and indeed the above mentioned densities just fail to belong to H −1/2 (Γ). Moreover, the densities h 3 and h 4 , and for "Dirac" also h 1 , just fail to belong to H 1/2 (Γ). This second strongest singularity observed is s −0.50000+1.25455i . Again h 1 and h 5 are zero for "HK 8-dens", modulo rounding errors.

Static plasmons
Similar to Section 10.1 we plot condition numbers of our IERs for the unit sphere to detect possible false eigenwavenumbers, but also to get a notion of how well-conditioned the two IERs under study are. This time we consider the quasi-static limit k ± → 0 of the IERs as in Section 7, and plot the condition numbers as functions of x = (1 +ˆ )/(1 −ˆ ) in Figure 10. With a slight abuse of notation, we keep speaking of (true/false) eigenwavenumbers.
As discussed in Section 7, the condition number of the (5:8,5:8) diagonal block for "Dirac", shown in Figure 10(c), follows closely that of K d and the peaks correspond to the static plasmons discussed in Section 8. The peak at x = −1/3 has a bright plasmon whereas all others have dark plasmons. "HK 8-dens" shows a symmetric spectrum with an infinite number of false eigenwavenumbers for 0 < x < 1. The full "Dirac" shows only two false eigenwavenumbers in the limits x = ±1. As shown in Figure 10(d), the peak at x = 1 can be removed by using the version of "Dirac" specified by (51). Note that "Dirac" specified by (35), is better conditioned near x = −1, where false eigenwavenumbers occur both for "Dirac" and "HK 8dens". For "Dirac" this false eigenwavenumber corresponds to the monopole field E = r/|r| 3 , H = 0 for |r| > 1, and E = H = 0 for |r| < 1.

Conclusions and discussion
Our numerical results show that "Dirac" wins over "HK 8-dens" in almost all tests. "Dirac" does not have any false near-eigenwavenumbers for any passive materials, whereas "HK 8-dens" exhibits such in the plasmonic case, and even false eigenwavenumbers in the reverse plasmonic case. We have seen that "HK 8-dens" exhibits false essential spectrum on domains with corners, 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 end by comparing it to other available IERs. We limit our discussion to nonmagnetic materials, although "Dirac" was formulated for general materials in [22], which we plan to cover in a forthcoming publication. The aspects of the IERs that we base our discussion on are the following. We now make a comparison with the "Debye" formalism in [9,10] and the "DFIE" formalism in [33] which, according to [33,Conclusions], are the current leading IERs when it comes to well-posedness for a wide range of passive materials. Using formalism from [22,31], "Dirac" uses Cauchy type integral representations for the Dirac equation DF = ikF for the electromagnetic multivector field F , whereas "Debye" employs a Dirac multivector potential F = DG + ikG, leading to a component-wise Helmholtz equation ∆G + k 2 G = 0, for which classical layer potentials apply. The electric and magnetic Gauss equations translate to certain conditions on these boundary layers, which via Hodge decompositions of tangential vector fields allows one to eliminate all but two scalar boundary densities for each domain. In the scattering situation of the present paper, this yields a "Debye" IER with four scalar densities. It is known that (D) it does not have dense-mesh low-frequency breakdown, and no topological low-frequency breakdown if complemented by extra equations. A drawback of "Debye" is that implementing the Hodge projections requires the numerical solution of a second order surface Laplace equation on Γ. An implementation of "Debye" for scattering by perfect conductors and piecewise smooth Γ is in [5], but (A) it has not been used for the MTP(k − , k + ) on piecewise smooth Γ. According to [10,Equation (2.19)], "Debye" also (B) uses compositions of two integral operators. Concerning (E), it is stated in [10,Conclusions] that "Debye" is invertible for all passive materials. However, this is impossible for any IER, which is equivalent to MTP(k − , k + ), due to the true eigenwavenumbers shown in Figure 3(e). Moreover, the main uniqueness result for "Debye", [9, Theorem 3.2], assumes that e { + } > 0 and relies on [28,Theorem 69], which assumes permittivities in the first quadrant [28, p. 261]. These assumptions do not cover all passive non-magnetic materials, which correspond to the square [0, π/2] 2 in Figure 1(b). It was recently shown in [11] though, that these uniqueness results can be sharpened to cover (arg(k − ), arg(k + )) = (0, π/2) in particular. However, numerical experiments revealing possible false near-eigenwavenumbers as in Figure 3(d) are not available. "DFIE" is formulated in [33]. Like "Debye", it employs potentials to reduce Maxwell's equations to vector Helmholtz equations. With "DFIE", the fields E and H are calculated independently of each other, by solving an 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 does not have false eigenwavenumbers for k − > 0 and 0 ≤ arg(ˆ ) < π, but it has false eigenwavenumbers at (arg(k − ), arg(k + )) = (0, π/2) as shown in [18,Section 11.3.2] (at least for a 2D version of "DFIE"). The theory in [33] is presented for Hölder boundary function spaces, and it is assumed that Γ is C 2 regular. However, by straightforward adaption to fractional Sobolev spaces and Lipschitz regular Γ as in [22], it appears that (A) "DFIE" is applicable to non-smooth scattering. Inspecting the integral operators that "DFIE" employs in [33,Appendix D], we see that (B) it uses single and double layer type boundary operators like "Dirac", with one exception. The exception is K 31 in [33, Equation (89)], which is a compact difference of hypersingular operators, and somewhat more numerically challenging to implement [25]. To compare (C) how fast "DFIE" and "Dirac" are, since the operators employed are roughly the same and having a good solver in mind which is linear in the number of non-zero operator blocks, we count these. From [33,Equation (87)] we have 2(36 − 4) = 64 non-zero operator blocks for "DFIE", keeping in mind that we need to solve two 6 × 6 systems for obtaining both E and H. From (27) we see that "Dirac" uses 64 − 14 = 50 non-zero operator blocks. For field evaluations, [33,Equations (36) and (38)] show that "DFIE" requires the evaluation of three double layer and three single layer potentials, for each component of the fields. For "Dirac", using the projected densities (34), which are fast to compute on Γ, each field component requires in general the evaluation of three double layer and two single layer potentials.
Comparing with "Dirac", this has the advantage of (A) applying to any Lipschitz regular Γ and (B) only using bounded integral operators of single and double layer type with explicit kernels. We remark that although (23) seems to be a Fredholm equation of the second kind, the operator G is not compact, even on smooth Γ. However G 4 is compact, which explains why "Dirac" works well in an iterative solver. Moreover, "Dirac" does not (E) have any false eigenwavenumbers, false near-eigenwavenumbers or false essential spectrum, (D) not even in the quasi-static limit, except at x = ±1. These quasi-static endpoints are related to eddy current scattering, and is the topic of a forthcoming paper by the authors. A reason why "Dirac" is so successful in avoiding false spectra, is that the choices of parameters make its invertibility depend only on the (6,6) diagonal block involving K d . Thus the false spectrum of K m , of which only one Hodge component is needed for the MTP(k − , k + ), is avoided.
In [2, Important Remark p. 123] it is stated that σ(K m ) \ σ(K d ) contributes to "higher-order resonances" in the near quasi-static limit. However, plots of the condition numbers for "Dirac" on the sphere near the quasi-static limit (not shown) are very similar to Figure 10 The entries of the matrix E k of (27) involve two families of integral operators denoted S k and K k . A given member in an operator family, S α k or K α k , is defined by its superscript α, which can be a constant, a unit vector, a scalar product of unit vectors, or a cross product of unit vectors. Specifically we have for S 1 k acting on a general density g