Dirac integral equations for dielectric and plasmonic scattering

A new integral equation formulation is presented for the Maxwell transmission problem in Lipschitz domains. It builds on the Cauchy integral for the Dirac equation, is free from false eigenwavenumbers for general complex-valued permittivities, can be used for magnetic materials, is applicable in both two and three dimensions, and does not suffer from any low-frequency breakdown. Numerical results for the two-dimensional version of the formulation, including examples featuring surface plasmon waves, demonstrate competitiveness relative to state-of-the-art integral formulations that are constrained to two dimensions. However, the new formulation applies to scattering also in three dimensions, where the theory suggests that it will perform equally well from a numerical point of view.


Introduction
This paper introduces a new boundary integral equation (BIE) formulation for solving the time-harmonic Maxwell transmission problem across interfaces between domains of constant and isotropic, but otherwise general complex-valued, permittivities and permeabilities. The interfaces may be merely Lipschitz regular. As the title suggests, we are primarily interested in non-magnetic materials with permittivities in the closed upper half-plane. However, our integral equation and analysis apply also to magnetic and acoustic scattering.
We refer to our new formulation as a Dirac integral equation, since the main idea behind it is to embed Maxwell's equations into a Dirac-type equation by relaxing the constraints divE = 0 and divB = 0 on the electric and magnetic fields E and B and introducing two auxiliary Dirac variables, which make the partial differential equation (PDE) elliptic also without these divergence-free constraints. This means that our integral equation uses eight unknown scalar densities in three dimensions (3D) and four unknown scalar densities in two dimensions (2D), where eight and four, respectively, are the dimensions of the Clifford algebra which is the algebraic structure behind the Dirac equation and which we make use of in the analysis. There are also independent evidence, see [17], that eight seems to be the number of scalar densities needed to construct a Maxwell integral equation in 3D that is free from false eigenwavenumbers for all complex-valued choices of material constants.
The theory of Dirac integral equations for solving Maxwell's equations was developed by the second author in his PhD work [2,3,1,6,4] with Alan McIntosh, building on earlier works by McIntosh and Mitrea [20] and Axelsson, Grognard, Hogan and McIntosh [5] The focus of these works is on the harmonic analysis and Fredholm operator theory of the integral equations. In a series of more recent works [24,22], the second author has further developed these integral equations with focus on numerical applications. In the present work, the main objective is to design the integral equations so that not only the operator has good Fredholm properties, but also so that no false eigenwavenumbers appear. This problem may be best explained by the relation PDE • Ansatz = BIE . The PDE problem, Helmholtz or Maxwell, is to invert a map (F + , F − ) → g, where g is boundary datum, and F ± are the solutions, the sought fields in the interior domain Ω + and exterior domain Ω − , respectively. To obtain a BIE we need to choose an appropriate ansatz, that is an integral representation h → (F + , F − ) of the fields in terms of a density function h on ∂Ω. From this we obtain a BIE as the composed map h → g, via (F + , F − ). A main objective is to use an invertible ansatz with as good condition number as possible, so that the BIE is invertible whenever the PDE problem is well posed.
Our main results, the Dirac integral equations for solving Maxwell's equations in R 2 as well as in R 3 , are formulated in Section 2 and make use of ansatzes which are invertible for all choices of material constants. At a more technical level we make the following comments, where k ± are wavenumbers in Ω ± , Im k ± ≥ 0, andk = k + /k − .
• Propositions 5.1 and 7.2 show sufficient conditions for the Helmholtz and Maxwell transmission problems to have unique solutions. For non-magnetic materials, our result is the hexagonal region shown in Figure 1, which strictly contains the regions earlier obtained for the Helmholtz transmission problem in [19,18]. See, however, [17,Remark 3.1] on a minor flaw in [19]. When also 0 /2 0 /2 Figure 1. The hexagonal region of (arg(k − ), arg(k + )) for whichk ∈ WP (k − , k + ), with WP (k − , k + ) from Definition 2.1.
the PDE problem defines a Fredholm map, something that can fail in the hexagon only when | arg(k + ) − arg(k − )| = π/2, then we also have existence of solutions. Outside this region of existence and uniqueness, and for a given k − , the uniqueness of solutions fails only for a discrete set of k + . • The Dirac BIE, presented in Section 2, is only one possible choice from a family of Dirac BIEs derived in Sections 6 and 7. In R 2 , these BIEs depend on four Dirac parameters r, β, α , β . The choice for r is critical for numerical performance, whereas the precise choices for the other parameters seem to be of lesser importance as long as they are chosen so that no false eigenwavenumbers appear. In R 3 , there are two additional parameters γ, γ .
The Dirac BIE is constructed in Propositions 6.3 and 7.6 using a dual PDE problem as ansatz. Our choice of dual PDE problem giving the Dirac BIEs in Section 2, has to effect that these BIEs are invertible whenever the PDE problem is well posed. • Propositions 5. 2 and 7.4 show that the PDE problem defines a Fredholm map provided that the quotientˆ of the permittivities is not in an interval [−C, −1/C] for some C ≥ 1 depending on the Lipschitz regularity of the interface. These estimates use Hodge potentials for the fields and concern the physical energy norm, corresponding to a boundary function space, H 2 or H 3 , which is a suitable mix of the function space H 1/2 (∂Ω) ⊂ L 2 (∂Ω) of traces to Sobolev H 1 (Ω) functions and its dual space H −1/2 (∂Ω) ⊃ L 2 (∂Ω). More precisely, in R 2 one can show that the PDE problem defines a Fredholm map when ±(ˆ + 1)/(ˆ − 1) is outside the essential spectrum of the double layer potential operator where Φ is the Laplace fundamental solution and ν is the outward unit normal. In R 3 , the Fredholm-map property of the Maxwell problem is controlled by the essential spectrum of K d together with the essential spectrum of the magnetic dipole operator acting on tangential vector fields f . This follows by writing out the results in [4] in classical vector notation. Indeed, the cosine operator (E k N + N E k )/2 appearing there, is the direct sum of operators K d , and in R 3 also K m . Compact perturbations of the operators K d and K m , and their adjoints, appear along the diagonals of our basic Cauchy integral operators (38) and (40). In particular, K m appears in the (3:4, 3:4) and (7:8,7:8) size 2 × 2 diagonal blocks in (40). It is important to note that in the energy norm, the essential spectra of K d and K m are subsets of (−1, 1). If considering instead the L 2 (∂Ω) norm, then on domains in R 2 with one corner, the essential spectrum of K d is a lying "figure eight", centered around 0. See Fabes, Jodeit and Lewis [11] and the final comments and (121) in the present paper.
• Although our main concern is wavenumbers k ± = 0, it is important to have a good behaviour of the BIEs as frequency ω → 0. Typical problems which can occur are dense-mesh low-frequency breakdown and topological lowfrequency breakdown. We refer to [25,10] for a more detailed discussion and for BIEs that are immune to such low-frequency breakdown, and remark that this immunity is shared by the Dirac BIEs presented in Section 2.
One reason for this is that our four densities include the gradient of the field in R 2 and our eight densities include the electric and the magnetic fields in R 3 , so that no numerical differentiation is needed. Moreover, if k ± → 0 withk constant in (10) or (20), then P, P , N, N are all constant whereas the Cauchy singular integral operators E k ± converge to E 0 in operator norm. The proof that the limit operator is invertible is outlined in Remark 7.3.
The paper ends with Section 8, where the numerical performance of the Dirac integral equations is studied in R 2 . Among other things, the performance is compared to that of a state-of-the-art system of integral equations of direct (Green's theorem method) type [18,17] which is only applicable in R 2 , which involves only two unknown scalar densities, and which for certain k ± coincides with a 2D version of the classic Müller system [21, p. 319]. Not surprisingly this special-purpose system performs best, but the new 4 × 4 Dirac system is not far behind. In particular it performs well under computationally challenging plasmonic conditions, that is when is negative and real.
In future work, we plan to test the new 8 × 8 Dirac system for the Maxwell transmission problem in R 3 using the numerical techniques in [17]. All BIE systems that we are aware of in the literature, with size less than 8 × 8, exhibit false eigenvawenumbers -primarily at the left middle corner point in Figure 1. In contrast, the theory in this paper suggests that the Dirac integral equations are entirely free from false eigenwavenumbers and will perform just as well in R 3 as in R 2 .

Matrix Dirac integrals
This section states, in classical vector and matrix notation, the integral equations which we propose for solving the Maxwell transmission problems. The derivation of these equations in later sections makes use of multivector algebra. We consider a bounded interior domain Ω + , separated by the interface ∂Ω from the exterior domain Ω − := R n \ Ω. Fix a large R < ∞ and define Ω − R := {x ∈ Ω − ; |x| < R}. For simplicity we assume throughout this paper that Ω − is connected, although many results do not need this assumption.
The parameters that we use for problem description in the transmission problems (4) and (14), that we aim to solve, are wavenumbers k + and k − in Ω + and Ω − respectively, and a jump parameterˆ . We assume that Im k ± ≥ 0, but unless otherwise stated we do not assume any relation between k + , k − andˆ . We denote the ratio between the wavenumbers byk = k + /k − . Our main interest is in scattering for non-magnetic materials whereˆ =k 2 . In applications, the parameterˆ appears as the ratioˆ = + / − of the permittivities ± in Ω ± . Similarly for magnetic materials, we have a ratioμ = µ + /µ − of the permeabilities µ ± in Ω ± . In this general case, we have the relationˆ =k 2 /μ.
We remark that in what follows, E ± denotes the standard electric fields, but we have rescaled the magnetic fields B ± so that what we here call B ± , in standard notation reads B ± / √ ± µ ± . To formulate our results, we need the following conical subsets of C. We use the argument −π < arg(z) ≤ π. Definition 2.1. Let k − , k + ∈ C\{0} be such that Im k − ≥ 0 and Im k + ≥ 0. Define Our notation is to denote fields in the domains by F, U, . . ., with suitable superscript ±, and to write f, u, . . . for respective boundary traces. On ∂Ω we denote by ν the unit normal vector pointing into the exterior domain Ω − . Furthermore {ν, τ } and {ν, τ, θ} denote positive ON-frames on ∂Ω depending on dimension, so that τ is the counter-clockwise tangent on curves ∂Ω ⊂ R 2 . On surfaces ∂Ω ⊂ R 3 , the theory we develop works for any choice of tangent ON-frame {τ, θ}. We denote the directional derivative in direction v by ∂ v , and with slight abuse of notation , ∂ ν u denotes normal derivative on ∂Ω although this use values of U in a neighbourhood of ∂Ω.
Our main result for 2D scattering concerns the Helmholtz transmission problem where Ω + ⊂ R 2 , n = 2, and with u 0 ∈ H 1/2 (∂Ω) being the trace of an incoming wave U 0 , and we want to solve for U + ∈ H 1 (Ω + ) and U − ∈ H 1 (Ω − R ). Define Theorem 2.2. Let Ω + ⊂ R 2 be a bounded Lipschitz domain, with Ω − being connected and notation as above. Then there exists a constant 1 ≤ C(∂Ω) < ∞ depending on the Lipschitz constants of the parametrizations of Ω + by smooth domains, so the following holds.
The transmission problem (4) is well posed if For its solution, consider the Dirac integral equation where E k is the singular integral operator (38) which we introduce in Section 4, P, P , N, N are the constant diagonal matrices (5)- (8), and The operator in (10) is invertible on the energy trace space H 2 from (56), introduced in Section 6, whenever (9) holds andk ∈ C \ (−∞, 0]. Moreover, the solution to (4) in Ω ± is obtained from h + = N h and h − = P h as Proof. This result is derived in Sections 5 and 6, where we arrive at equation (82) with α =ˆ . We precondition (82) by multiplying from the left by P and writing h = P h to obtain (10) with N = P M , N =kM P and P (kM + M )P = I. We remark that the parameter ratiosk < 0 andˆ = −1 must be excluded for P and P to be well defined. However, the only problem whenk < 0 is that we cannot precondition (82) with P and P to achieve a non-integral term I in (10).
In applications of the Dirac problem (55) to the Helmholtz problem (4), as in Example 6.1, we have f 0 = ik − u 0 + ∇u 0 , that is (11). Equations (12) and (13) for the fields are seen from the first and last two rows in (38) respectively, and the expressions for h ± are seen from the right factor in (75). Note that h − happens to coincide with the densityh, which is introduced in (78) in a different context. We next formulate the analogous result for the Maxwell transmission problem in R 3 . Our transmission problem is where E 0 and B 0 are the incoming electric and magnetic fields, and we want to solve for E ± and B ± . All fields are assumed to be L 2 integrable in a neighbourhood of ∂Ω.
Define matrices (where we writeĉ = 1/k andμ =k 2 /ˆ ) Let Ω + ⊂ R 3 be a bounded Lipschitz domain, with Ω − being connected and notation as above. Then there exists a constant 1 ≤ C(∂Ω) < ∞ depending on the Lipschitz constants of the parametrizations of Ω + by smooth domains, so the following holds.
The transmission problem (14) is well posed if For its solution, consider the Dirac integral equation where E k is the singular integral operator (40) which we introduce in Section 4, P, P , N, N are the constant diagonal matrices (15)- (18), and with field components in the frame {ν, τ, θ}. The operator in (20) is invertible on the energy trace space H 3 from (85), introduced in Section 7, whenever (19), holds. Moreover, the solution to (14) in Ω ± is obtained from h + = N h and h − = P h as whereK andS are the layer potentials (41) and (42) from Section 4, with v × · denoting the map x → v × x.
Proof. This result is derived in Section 7, where we arrive at equation (117) with α =ˆ . We precondition (117) by multiplying from the left by P and writingh = P h to obtain (20) with N = P M , N = M P and P (k −1 M + M )P = I. We remark that the parameter ratiosk < 0,μ =k 2 /ˆ = −1 andˆ = −1 must be excluded for P and P to be well defined. However, the only problem whenk < 0 is that we cannot precondition (117) with P and P to achieve a non-integral term I in (20). In applications of the Dirac problem (55) to the Maxwell problem (14), or in multivector notation (89), we have F 0 = E 0 + B 0 , that is (21). Equations (22) and (23) for the fields are seen from the first and last two rows in (40) respectively.

Multivector algebra
This section contains a very brief explanation of the multivector algebra which we use. For a complete account of the theory of multivector algebra and Dirac equations which we use, we refer to [23], in particular Chapter 9 therein. Roughly speaking multivectors are the same as Cartan's alternating forms, and multivector fields are same as differential forms, and amounts to an algebra of not only the one-dimensional vectors but also k-dimensional algebraic objects for 0 ≤ k ≤ n in R n . Concretely, multivectors in R n are the following type of 2 n dimensional objects, where our interest is in n = 2, 3. Denote by {e 1 , . . . , e n } the standard vector basis for R n . We write ∧R n for the complex 2 n dimensional space of multivectors in R n , which is spanned by basis multivectors e s , where s ⊂ {1, 2, . . . , n}. We write ∧ j R n for the subspace spanned by e s , where j is the number of elements in the index set s. In R 2 , a multivector is of the form and so amounts to two scalars a, b and a vector v = v 1 e 1 + v 2 e 2 . In R 3 , a multivector is of the form (25) w = a + v 1 e 1 + v 2 e 2 + v 3 e 3 + u 1 e 23 + u 2 e 31 + u 3 e 12 + be 123 , and so amounts to two scalars a, b and two vectors v = v 1 e 1 + v 2 e 2 + v 3 e 3 and u = u 1 e 1 + u 2 e 2 + u 3 e 3 via Hodge duality.
We use three products on ∧R n : the exterior product u ∧ w, the (left) interior product u w and the Clifford product uw, corresponding to signed unions, set differences and symmetric differences of index sets respectively. We refer to [23,Chapters 2,3] for details, and remark that our convention is that of positive Clifford squares so that (26) e i e j + e j e i = 2 e i , e j .
When s = {s 1 , . . . , s k } with s 1 < . . . < s k , then the basis multivector e s equals the Clifford (as well as exterior) product of the basis vectors e s 1 , . . . , e s k . It is mainly the case when the first factor u is a vector, which we use. Replacing w by a multivector field F : R n → ∧R n and u by the nabla symbol ∇, we obtain differential operators ∇ ∧ F (x), ∇ F (x) and DF (x) = ∇F (x): the exterior, interior and Dirac derivatives. We refer to [23,Chapters 7,8] for details.
The time-harmonic wave Dirac equation which we consider in this paper is for multivector fields F . Note that in (27), the product between e j and the multivector ∂ e j F (x) is the Clifford product, which for example in 3D is a combination of vector and scalar products. The main application are Maxwell's equations. Define the total electromagnetic (EM) field to be the multivector field where E is the electric field and B is the magnetic field, and the ∧ 0 R 3 and ∧ 3 R 3 parts are zero. In this formalism, the Dirac equation (27) for F = F em given by (28) coincides with the time-harmonic Maxwell's equations. Indeed, the ∧ 0 R 3 ,

The R n Cauchy integral
A main reason for using the Dirac framework is that it provides us with a Cauchytype reproducing formula, which allows for a generalization of complex function theory to n ≥ 2 and k = 0. See [23,Chapters 8,9] for further details. More precisely, if F satisfies (27) in a domain Ω with boundary ∂Ω, then a Cauchy-type reproducing formula holds. We write dy and ν ∈ ∧ 1 R n for the standard measure and outward pointing unit normal on ∂Ω respectively, and the integrand uses two Clifford products. The first factor Here Φ k (x) ∈ ∧ 0 R n = C is the Helmholtz fundamental solution, and we use the normalization from [8] so that (∆ + k 2 )Φ k = −2δ(x). Hence the factor −1/2 in (30). In dimension n = 3 we have and in dimension n = 2 we have in terms of the Hankel function H 0 (k|x|). The classical theory of complex Hardy spaces generalizes from complex function theory to our Dirac setting. Our basic operator, acting on a suitable space H of functions h : ∂Ω → ∧R n on ∂Ω, is the Cauchy principal value integral which reduces to the classical Cauchy integral when n = 2, k = 0. The basic operator algebra is that E 2 k = I, and (34) The operator E ± k projects onto its range, the subspace of H which we denote E ± k H and consists of all traces F | ∂Ω of fields satisfying DF = ikF in Ω ± . For the exterior domain Ω − these fields also satisfy a Dirac radiation condition; See (55) below for its formulation.
For computations, we express the Cauchy singular operator E k as a matrix with entries being double and single layer potential-type operators, using the notation respectively. Here v(x, y) is a vector field and a(x, y) is a scalar function, depending on x, y ∈ ∂Ω. The single layer S a k is a weakly singular integral operator, and so is also K v k if ∂Ω is smooth and v = ν is the normal direction at x or y. Otherwise K v k is a principal value singular integral, but is bounded on many natural function spaces H, also for general Lipschitz interfaces ∂Ω.
Consider first dimension n = 2, with a curve ∂Ω and Φ k given by (32). Here we use a positively oriented frame {ν, τ } at x ∈ ∂Ω, with ν = ν(x) being the normal vector into Ω − and τ = τ (x) the tangential vector counter clock-wise from ν. The corresponding frame at y ∈ ∂Ω we write as ν = ν(y) and τ = τ (y). In the plane, the Clifford algebra ∧R 2 is four dimensional and spanned by {1, e 1 , e 2 , e 12 }. We prefer the ordering {1, e 12 , e 1 , e 2 }, since Clifford multiplication by vectors then will be represented by block off-diagonal matrices. At x ∈ ∂Ω, we write using instead the vector frame {ν, τ }. Here ντ = e 12 ∈ ∧ 2 R 2 does not depend on x, even though ν and τ do so. By writing out ∇Φ k and the Clifford products in (33), we obtain in the multivector frame {1, ντ, ν, τ }. Next consider dimension n = 3, with a surface ∂Ω and Φ k given by (31). Here we use a positively oriented ON-frame {ν, τ, θ} at x ∈ ∂Ω, with ν = ν(x) being the normal vector into Ω − and τ = τ (x) and θ = θ(x) being tangential vector fields. The frame vectors at y ∈ ∂Ω, we write as ν = ν(y), τ = τ (y) and θ = θ(y). A multivector field at x ∈ ∂Ω we write as Here ντ θ = e 123 ∈ ∧ 3 R 3 does not depend on x, although in general each of the three vectors do so. Again, we prefer this ordering of the frame multivectors since Clifford multiplication by vectors then is block off-diagonal. In the multivector frame {1, τ θ, θν, ντ, ντ θ, ν, τ, θ}, we have Finally we note that we similarly can write the Cauchy integral (29) for the fields in Ω ± , in matrix form. The only difference is the normalization factor 2 and the fact that we choose the frame {e 1 , . . . , e n } at the field point x ∈ Ω ± . We also allow for a general function h : ∂Ω → ∧R n and not only a trace f of a solution to the Dirac equation in (29). By the associativity of the Clifford product, this still yields a field F solving the Dirac equation, but h = F | ∂Ω in general. In this case, when x / ∈ ∂Ω, we denote the layer potentials byK v k andS a k , with v and a now only depending on y ∈ ∂Ω. We also need generalizations of these potentials, and writẽ where A(y) : R n → R n is a matrix function and v(y) is a vector field.

Helmholtz existence and uniqueness
This section surveys basic solvability results for the Helmholtz transmission problem (4), which are valid in any dimension n ≥ 2. The proofs follow [4], with a translation from the Dirac to the Helmholtz framework. Consider first uniqueness of solutions U ± .
Proposition 5.1. Assume that Ω − is connected and that the incoming wave vanishes so that u 0 = 0. If Proof. From the jump relations we have Apply Green's first identity for Ω + and Ω − R , use the Helmholtz equation for U ± and the radiation condition for U − , and multiply the equation byˆ k − /i to obtain as R → ∞. Denote the three integrals appearing in (45) by I + , I − R and I R respectively, including i −1 in the first two, and set φ ± := | arg(k ± /i)| as in Definition 2.1, let z :=ˆ k − /k + , and define the sector We note that (47) and similarly for I − R . We verify that the condition (43) implies U ± = 0, by examining (45) in the nine cases φ ± = 0, 0 < φ ± < π/2 and φ ± = π/2 as follows.
If also φ + = π/2 then Re I + = 0, and we conclude that lim R→∞ I R = 0 also when z < 0, and can in the same way conclude that U − = 0 = U + .
Next consider the existence of solutions U ± .
Proposition 5.2. Let Ω be a bounded Lipschitz domain. Then there exists 1 ≤ C(∂Ω) < ∞ such that if Proof. It suffices to estimate (49) by the norm of u 0 and compact terms. Indeed, perturbation theory for Fredholm operators applies to show that the index of the map (U + , U − ) → u 0 is zero, by perturbing the parameters (k − , k + ,ˆ ) to (i, i, 1), where the map clearly is invertible.
To establish the estimate of (49), we construct certain auxiliary potentials V ± to the gradient vector fields ∇U ± . In Ω − R , we simply use the given scalar potential V − = U − . In Ω + , we find a bivector field V + : Ω + → ∧ 2 R n and a vector field (Modulo the termṼ + , this means that V + is a conjugate function when n = 2 and a vector potential when n = 3.) The existence and compactness of the map H 1 (Ω + ) → L 2 (Ω + ) 2 : U + → (V + ,Ṽ + ) follows from Hodge decompositions as in [6], after translation from the spacetime framework using [23, Sec. 9.1]. To complete the construction of V ± , we extend the potentials V ± to compactly supported functions on R n , with ∇V − and ∇ V + belonging to L 2 (R n ).
Pairing the jump relations with v ± , we have and (54) Therefore, adding the equations (51) yields an estimate of (49) whenever we do not haveˆ < 0. Whenˆ << 0, we instead subtract the equations to conclude. When ≈ 0 we can also obtain an estimate of (49) by instead starting with a bivector potential V + in Ω + , and a scalar potential V − in Ω − .

The 2D Dirac integral
Let Ω + ⊂ R n be a bounded Lipschitz domain and consider the Dirac transmission problem (55) This is our master transmission problem, into which we embed the Helmholtz and Maxwell transmission problems.
Throughout this section n = 2. In the frame (37), our boundary function space for f 0 and f ± is (56) In the domains Ω + and Ω − R this corresponds to where F ± 0 and F ± 2 denote the scalar and bivector parts of F ± respectively, that is the first two components in the basis {1, e 12 , e 1 , e 2 }. The matrix M , to be chosen to be constant and diagonal in the frame (37), specifies jump relations across ∂Ω. Example 6.1 (Transverse magnetic (TM) scattering). We consider applications of the Dirac equation to the scattering of EM fields as in (28), but independent of the e 3 -coordinate and polarized so that To write a Helmholtz equation for this EM field, we normalize by a left Clifford multiplication and define the field where ± is suppressed. Since the Clifford product is associative, we have Writing F = ikU + V 1 e 1 + V 2 e 2 , the Dirac equation DF = ikF amounts to V = ∇U and divV = −k 2 U , that is the Helmholtz equation ∆U + k 2 U = 0 for the scalar function U . We have With this setup for both domains Ω ± , jump relations for the electromagnetic field specify the jump matrix (65) M = diag k aˆ 1 for F ± in the frame {1, ντ, ν, τ }. The parameter a can be chosen freely since F 2 = 0 for the field F .
The following Dirac well-posedness result exploits that the Dirac equation in the plane consists of two coupled Helmholtz equations.  Then the operator B 2 : Proof. Consider the equation with a given source f 0 , and write f + = f + 0 + f + 1 + f + 2 e 12 with scalar functions f + 0 and f + 2 and a vector field f + 1 , and similarly for f − and f 0 . We first prove uniqueness in two steps as follows. To this end we assume that f 0 = 0.
(i) The scalar functions F ± 2 solve the Helmholtz equation as a consequence of F ± solving the Dirac equation, with wavenumbers k ± respectively. Moreover, the vector part of DF ± = ik ± F ± shows that (70) ∇F ± 2 = (∇F ± 0 − ik ± F ± 1 )e 12 . From the assumed jump relations for f ± it therefore follows that f + 2 =kf − 2 /β and ∂ ν f + 2 = β∂ ν (kf − 2 /β), and hence Proposition 5.1 withˆ = β shows that f ± 2 = 0. (ii) Next consider the scalar functions F ± 0 , which also solve the Helmholtz equation. Since F ± 2 = 0 by (i), we have ∇F ± 0 = ik ± F ± 1 and obtain jump relations . Again Proposition 5.1 applies, now withˆ = α, and shows that f ± 0 = 0. From (70) we conclude f ± 1 = 0 and in total f ± = 0. To show existence, it suffices by perturbation theory for Fredholm operators to prove an estimate This follows as in steps (i) and (ii) but instead appealing to Proposition 5.2. In this case, we obtain in step (i) that f + 2 The next result is central to this paper, where we derive Dirac integral equations by using an ansatz obtained from an auxiliary Dirac transmission problem via duality. Then the operator By Proposition 6.2, the left factor is invertible, so it suffices to show that the right factor also is invertible. To this end, we use the (non-Hermitean) complex bi-linear duality on H 2 , where w = j (−1) j w j denotes the involution of a multivector w = j w j , w j ∈ ∧ j R n . It is readily verified that where M = diag 1/k (1/k)/β α 1 . This shows that in the natural way E ∓ k ± H 2 is the dual space of E ± k ± H 2 . Hilbert space duality theory shows that in (75), invertibility of the right factor is equivalent to invertibility of the left factor, with M replaced by M , and k − and k + swapped.

The 3D Dirac integral
In Sections 5 and 6, we derived an integral equation for solving Dirac transmission problems in R 2 , which applies to Helmholtz/TM Maxwell scattering. We here derive the completely analogous integral equation in R 3 , with applications to scattering for the full Maxwell equations and not only the Helmholtz equation.
We set out to formulate an integral equation for solving the Dirac transmission problem (55) when n = 3. In the frame (39), our boundary function space for f 0 and f ± is (85) where curl denotes surface curl on ∂Ω. In the domains Ω + and Ω − R this topology corresponds to where F j denotes the ∧ j R 3 part of the field. Example 7.1. Maxwell's equations correspond to an electromagnetic field F with F 0 = 0 = F 3 as in (28). The energy norm that we consider is simply the L 2 norm of F ± in Ω + and Ω − R , respectively, and the corresponding function space H 3 on ∂Ω is (H −1/2 (∂Ω)) 6 , with tangential curls of E and B belonging to H −1/2 (∂Ω).
In the 3D Dirac transmission problem (55), Maxwell scattering for the field F = F em from (28) specifies the jump matrix Consider the Maxwell transmission problem (14), which we write in multivector notation, with F ± 1 = E ± and F ± 2 being the Hodge dual of B ± , as follows. (89) We require the following Maxwell versions of the results in Section 5.
Proposition 7.2. Assume that Ω − is connected and that the incoming wave vanishes so that f 0 = 0. If (90)ˆ /k ∈ WP (k − , k + ), Proof. Similar to the proof of Proposition 5.1, we use the jump relations to obtain We then apply a Stokes theorem for Ω + and Ω − R , to obtain Using (92), the result follows similarly to the proof of Proposition 5.1.
Remark 7.3. Consider the quasi-static limit k ± = 0, with fixed jumpk, of Maxwell's equations (89). Two modifications are needed. The jumps (100) and (102) below specified by the Gauss law do not follow from (89) as when k ± = 0, but need to be assumed. We also need to modify the radiation condition, since in the quasi-static limit the exterior field has decay O(|x| −2 ).
The estimates of the essential spectrum in Proposition 7.4 below go through with minor changes. The proof of uniqueness of solutions in Proposition 7.2, however, presents some novelties. Note that in the quasi-static limit, the electric and magnetic fields F 1 and F 2 decouple. For the electric field, we have F ± 1 being divergence-and curl-free vector fields, with jumps To show that this forces F + = 0 = F − , we note that is curl-free in all R 3 in a distributional sense. This shows that F = ∇U in R 3 for some function U ∈ H 1 loc (R 3 ), and from the decay condition on F − we conclude that the integration constant can be chosen so that U = O(|x| −1 ) as x → ∞. Since ν F ± = ∂ ν U ± , with U ± := U | Ω ± , we can proceed as in Proposition 5.1 to obtain (95) by the norm of f 0 and compact terms. The proof is similar to that of Proposition 5.2, replacing the scalar function ik ± u ± by the vector field F ± 1 , and the vector field ∇u ± by the bivector field F ± 2 .
In Ω − R we construct potentials such that and in Ω + we construct potentials such that where subscript j refers to a ∧ j R 3 valued field. The existence and estimates of such potentials follow from Hodge decompositions as in [6], after translation from the spacetime framework using [23, Sec. 9.1].
We consider first F ± 1 , where we pair the jump relation for ν ∧ f ± 1 with u + 2 . From the jump relation for ν f ± 2 and Maxwell's equations, we obtain (100) Next consider F ± 2 , where we pair the jump relation for ν f ± 2 with u − 1 . From the jump relation for ν ∧ f ± 1 and Maxwell's equations, we obtain (102) , which we pair with u + 3 . We get Applying the general Stokes theorem, see [23,Sec. 7.3], to (101) and (103) and summing the so obtained estimates yield an estimate of (97), similar to the proof of Proposition 5.2. A main idea is that the potentials U ± j depend compactly on the fields F ± j , and for details of the estimates we refer to [4,Lem. 4.9,4.17] where all the details are fleshed out.
With this solvability result for Maxwell's equations, we next derive solvability results for the 3D Dirac equation similarly to what was done in 2D in Section 6.  Then the operator B 3 : Proof. The proof is similar to that of Proposition 6.2, but now using that in 3D the Dirac equation contains two Helmholtz equations along with Maxwell's equations. We write and proceed in three steps. We first prove uniqueness.
(ii) Writing F ± 3 = U ± 3 e 123 , as in (i) the scalar functions U ± 3 solve the Helmholtz equation. From the Dirac equation we have From this and (104), we conclude that ∂ ν u + 3 = γ∂ ν (u − 3 /γ) and u + 3 = u − 3 /γ and hence Proposition 5.1 withˆ = γ shows that F ± 3 = 0 = U ± 3 . (iii) From (i) and (ii) we conclude that F ± solve Maxwell's equations (89) witĥ = α, and Proposition 7.2 applies to show that F ± = 0. To show existence, it suffices by perturbation theory for Fredholm operators to prove an estimate This follows similarly to steps (i)-(iii) by instead appealing to Propositions 5.2 and 7.4. Then the operator Proof. This follows by duality from Proposition 7.5, entirely analogously to the proof of Proposition 6.3.
Consider the 3D Dirac integral equation involving the operator from (113), with right hand side specified by the jump condition in (55) and auxiliary densityh which we precondition ash = P h in the proof of Theorem 2.3. Similar to what we did for the 2D twin (78), we optimize (114) by choosing the parameters β, γ, α , β , γ , r. Recall that for EM fields α =ˆ . We therefore consider α as having a prescribed value. Writing (113) as in (79), we have in R 3 that (115) rM +M = diag r α +k αβ r γ + 1 k r +k α r +k α rk + 1 γ r kα β + 1 α r α k + 1 r α k + 1 , and, modulo lower order operators of the form K v k 1 − K v k 2 and S a k 1 − S a k 2 , the integral operator T = rE k + M − M E k − from (79) is the entry-wise product of E k and (116)  It is the diagonal 4 × 4 blocks in (116) which are our main concern, within which we note that the diagonal 2 × 2 blocks are weakly singular operators on smooth domains.
To summarize, for solving the Maxwell transmission problem (14), we have obtained the 3D Dirac integral equation

Numerical results for the 2D Dirac integral equation
This section shows how the 2D Dirac integral equation (10), along with the field representation formula (12), performs numerically when applied to the planar Helmholtz/TM Maxwell transmission problem. For comparison, we also investigate the performance of the 4 × 4 system of integral equations [17,Eq. (85)], along with its field representation formulas [17,]. For simplicity, we refer to the system (10) as "Dirac" and to the system [17,Eq. (85)] as "HK 4-dens". The reason for comparing with "HK 4-dens" is that this system, just like "Dirac", has a 8 × 8 counterpart in 3D which applies to Maxwell's equations.
We also compare to a state-of-the-art 2 × 2 system of integral equations [17,Eq. (88)] of Kleinman-Martin type [18] and to the 2D version of the classic Müller system [21, p. 319], which also is a 2 × 2 system, and refer to these systems as "best KM-type" and "2D Müller". While "best KM-type" is limited to planar problems it yields, in general, the best results. For interior wavenumbers arg(k + ) = 0, "best KM-type" coincides with "2D Müller".
All systems of integral equations are discretized using Nyström discretization and underlying composite 16-point Gauss-Legendre quadrature. On smooth ∂Ω we use a variant of the "Nyström scheme B" in [15]. In the presense of singular boundary points on ∂Ω such as corner vertices, which would cause performance degradation in a naive implementation, the Nyström scheme is stabilized and accelerated using recursively compressed inverse preconditioning (RCIP) [14]. Accurate evaluation of singular operators on ∂Ω and accurate field evaluation of layer potentials at field points close to ∂Ω are accomplished using panel-wise product integration. See [16,Sec. 4], and references therein, for details. Large discretized linear systems are solved iteratively using GMRES.
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. Fields are computed at 10 6 points on a rectangular Cartesian grid in the computational domains shown. When assessing the accuracy of computed fields we compare to a reference solution. The reference solution is either obtained from a system deemed to give more accurate solutions, or by overresolution using roughly 50% more points in the discretization of the system under study. The absolute difference between the original solution and the reference solution is called the estimated absolute error.
8.1. The operators. We compute condition numbers of the discretized system matrices in the systems under study. One purpose is to detect false eigenwavenumbers. Another purpose is to compare the conditions number of the system matrices with each other. The interface ∂Ω is the smooth starfish-like curve [16, eq. (92)], shown in Figure 2(a) and originally suggested for scattering problems in [12]. A number of 976 discretization points are placed on ∂Ω. We study three cases of jumpsk where we vary k − or k + : • The positive dielectric case. The exterior wavenumber is positive real, 0 < k − ≤ 20, andk = k + /k − = 1.5 so that 0 < k + ≤ 30 andˆ = 2.25. This corresponds to the lower left corner point in Figure 1, where Theorem 2.2 guarantees that no true eigenwavenumbers exist, as well as no false eigenwavenumbers for "Dirac".  shows that "Dirac" and "HK 4-dens" perform equally well, except at low frequencies where the condition number of "Dirac" fares better and is comparable to that of "2D Müller".
• The plasmonic case. Again the exterior wavenumber is positive real, 0 < k − ≤ 20, butk = k + /k − = i √ 1.1838 so that 0 < k + /i ≤ 20 √ 1.1838, k + is imaginary, andˆ = −1.1838 is rather close to the essential spectrum {−1}. This corresponds to the left middle corner point in Figure 1, where Theorem 2.2 guarantees that no true eigenwavenumbers exist, as well as no false eigenwavenumbers for "Dirac". Figure 2(c) shows that the condition number of "Dirac" is closer to that of "best KM-type" than to "HK 4-dens", in particular at high frequencies. "2D Müller" exhibits 12 false eigenwavenumbers.
• The reverse plasmonic case. Now the interior wavenumber is positive real, 0 < k + ≤ 20, andk = k + /k − = (i √ 1.1838) −1 so that 0 < k − /i ≤ 20 √ 1.1838, k − is imaginary, andˆ = −1.1838. This corresponds to the lower middle corner point in Figure 1, which does not belong to the hexagon, and indeed Figure 2(d) shows 37 true eigenwavenumbers. The different systems here have a relative performance similar to that of the plasmonic case, with "Dirac" performing rather close to "2D Müller".  8.2. Field computations. We solve the Dirac system (10) for h, compute interior and exterior fields U ± via (12), and compare with results from the other systems. Gradient fields ∇U ± are not computed, but we remark that the representation formula (13) for ∇U ± uses layer potentials with the same type of (near-logarithmic and near-Cauchy) singular kernels as the representation formula (12). The 2 × 2 systems, on the other hand, have accompanying representation formulas for ∇U ± with layer potentials that contain near-hypersingular kernels. Evaluation of these layer potentials at field points near ∂Ω may cause additional loss of accuracy. The curve ∂Ω is chosen as the boundary of the one-corner drop-like object [16,Eq. (92)], discernible in Figure 3(a). We let k − = 18 in the positive dielectric and plasmonic cases, and k + = 18 in the reverse plasmonic case. The incoming wave is a plane wave from south-west, u 0 (x, y) = e ik − (x+y)/ √ 2 , and 800 discretization points are placed on ∂Ω. Whenˆ = −1.1838, then ±(ˆ + 1)/(ˆ − 1) is in the essential H 1/2 (∂Ω)-spectrum of (1) and a homotopy-based numerical procedure is adopted whereˆ = −1.1838 is approached from above in the complex plane [13, Sec. 6.3].
• Figure 3 covers the positive dielectric case. The real parts of the scattered fields U ± are shown in Figure 3 . Field computation in the plasmonic case; (a) the scattered fields Re U ± ; (b) log 10 of the estimated absolute error using "HK 4-dens"; (c) same using "Dirac"; (d) same using "best KM-type".
systems. Most likely, this is because (12) does not exploit null-fields in the near-field evaluation. See [17,Sec. 7]. The number of GMRES iterations needed to meet a stopping criterion threshold of machine epsilon in the relative residual are 61, 65, and 44 for "HK 4-dens", "Dirac", and "2D Müller", respectively. In this case the operators in "HK 4-dens" and "Dirac" seem to have similar spectral properties, while the spectral properties of the operator in "2D Müller" are better.
• Figure 4 covers the plasmonic case and is organized as Figure 3. There are propagating waves in Ω − , exponentially decaying waves into Ω + , and a surface plasmon wave along ∂Ω. "Dirac" here performs almost on par with "best KM-type" and gives 2 more accurate digits than "HK 4-dens". The number of GMRES iterations needed are 266, 173, and 143 for "HK 4dens", "Dirac", and "best KM-type", respectively. In this case the operator in "Dirac" seems to have considerably better spectral properties than the operator in "HK 4-dens".
• Figure 5 covers the reverse plasmonic case. The results are similar to those of the plasmonic case, although a digit is lost with all systems and we have propagating waves in Ω + and exponentially decaying waves into Ω − . "Dirac" performs almost on par with "2D Müller" and gives 2 more accurate digits than than "HK 4-dens".  . Same as Figure 6, but for the plasmonic case.  Figure 8. Same as Figure 6, but for the reverse plasmonic case.
meaning that h 1 , h 2 ∈ H 1/2 (∂Ω) and h 3 , h 4 ∈ H −1/2 (∂Ω). The result is shown in Figure 6, where indeed h 1 and h 2 are seen to be continuous at the corner vertex (note that h 2 ≈ 0). The only singular density is h 3 , which is related to that it is only the (3, 3) diagonal element in the 4×4 block-operator of (10) which we do not control by the choices of parameters r, β, α , β in Section 6. Using the automated eigenvalue method of [14,Sec. 14], the asymptotic behaviour of h 3 near the corner is determined to be where η = −0.12319432456634 and t is the arc length distance to the corner vertex. So h 3 is in fact in L 2 (∂Ω). • The plasmonic case: Here the hypothesis onˆ andk in Theorem 2.2 is not satisfied sinceˆ < 0 makes ±(ˆ +1)/(ˆ −1) hit the essential H −1/2 (∂Ω)-spectrum of (1). Nevertheless, the RCIP-accelerated Nyström scheme manages to produce the limit solution h shown in Figure 7. As in the positive dielectric case, the densities h 1 , h 2 , h 4 are good, although h 4 exhibits an oscillatory behaviour. However h 3 / ∈ H −1/2 (∂Ω). More precisely, its asymptotics near the corner are as in (120) with η = −1.00000000000000 − i1.57105873276994. So h 3 ∈ H −s (∂Ω) for any s > 1/2.
• The reverse plasmonic case: the results, shown in Figure 8, are very similar to those of the plasmonic case. The asymptotics of h 3 are as in (120) with η = −1.00000000000000 + i1.57105873276994. We end with a remark on the densities h obtained in the plasmonic and reverse plasmonic cases, which fall outside the energy trace space H 2 from (56). More generally, this energy trace space belongs to a family of function spaces, where Sobolev regularity s = 1/2 and s − 1 = −1/2 is replaced by a more general regularity index s. On Lipschitz domains, the possible range is 0 ≤ s ≤ 1. In the plasmonic cases, our computed densities h belong to the larger spaces s < 1/2. For 0 ≤ s < 1, the corresponding norms of the fields are weighted Sobolev norms using |∇U ± (x)| 2 dist(x) 1−2s dx, where dist(x) denotes distance from x to ∂Ω, whereas for the endpoint s = 1, this must be replaced by a norm involving a non-tangential maximal function. A reference for the elementary results for 0 < s < 1 is Costabel [9]. The bounds on double layer potential operators for the endpoints s = 0, 1 require harmonic analysis and the Coifman-McIntosh-Meyer theorem [7].
The essential spectrum of the double layer potential operator (1) depends on the choice of function space, that is on s. For the energy trace space s = 1/2 this spectrum is a subset of the real interval (−1, 1), but in the endpoint spaces s = 0, 1, as alluded to in the introduction, this spectrum can be computed to be a lying "figure eight", parametrized as (121) ± sin(δπ(1 + iξ)/2)/ sin(π(1 + iξ)/2), ξ ∈ R, where δ = θ/π − 1 on a domain with a corner of opening angle θ. The point of departure for the investigations reported on in this paper was the spin integral equation proposed in [22,Sec. 5] for solving the Maxwell transmission problem (14). However, it was soon realized that this was not suitable for the plasmonic and reverse plasmonic cases. Indeed, the theory developed for these spin integral equations makes use of non-diagonal matrices P, P , N, N in (10) which mix H ±1/2 (∂Ω) and is limited to the function space L 2 (∂Ω), that is the end point space s = 1. As we have seen, surface plasmon waves appear in the larger function spaces s < 1/2. The numerical algorithm used here becomes unstable for the spin integral equation when the "figure eight" of (121) is approached. Whenˆ = −1.1838 is approached from above in the complex plane, this happens nearˆ = −1.1838 + i0.2168.