Boosting the Maxwell double layer potential using a right spin factor

We construct new spin singular integral equations for solving scattering problems for Maxwell's equations, both against perfect conductors and in media with piecewise constant permittivity, permeability and conductivity, improving and extending earlier formulations by the author. These differ in a fundamental way from classical integral equations, which use double layer potential operators, and have the advantage of having a better condition number, in particular in Fredholm sense and on Lipschitz regular interfaces, and do not suffer from spurious resonances. The construction of the integral equations builds on the observation that the double layer potential factorises into a boundary value problem and an ansatz. We modify the ansatz, inspired by a non-selfadjoint local elliptic boundary condition for Dirac equations.


Introduction
The classical principal value double layer potential is the operator Kf (x) = 2p.v. ∂Ω (∇Φ)(y − x) · ν(y)f (y)dσ(y), x ∈ ∂Ω, acting on functions f the boundary of a bounded Lipschitz domain Ω ⊂ R n , and using the Laplace fundamental solution Φ and the outward pointing unit normal vector field ν for its kernel. Here σ is standard surface measure and we choose to normalize by a factor 2. The method of boundary integral equations for solving the classical Dirichlet and Neumann boundary value problems for the Laplace equation on Ω is as follows. To solve the Dirichlet problem with datum g 1 , we solve from which the harmonic function u is obtained as the double layer potential with density f 1 . Similarly, to solve the Neumann problem with datum g 2 , we solve from which the harmonic function u is obtained as the single layer potential with density f 2 . For smooth domains, K is a weakly singular integral, which gives a compact operator on many function space and invertibility can be deduced by classical Fredholm theory. For Lipschitz domains, K is a singular integral operator (modulo the factor ν), and its L p boundedness, 1 < p < ∞, follows from the seminal work by Coifman, McIntosh and Meyer [8]. For the rest of this paper, we will restrict attention to the most fundamental function space for singular integrals: L 2 . On a strongly Lipschitz domain, that is when ∂Ω is locally a Lipschitz graph, Rellich identities replaces the Fredholm arguments to show that I ± K is a Fredholm operator on L 2 (∂Ω).
In this paper we derive a new integral equation for solving the Maxwell scattering problem again perfect conductors. This makes use of Clifford algebra and an embedding of Maxwell's equations into a Dirac equation. To explain our ideas and results, we discuss in this introduction the double layer potential in the complex plane where all the main ideas are present but the algebra is simpler: Clifford algebra simplifies to complex algebra, Maxwell's equations simplify to Cauchy-Riemann's equations, and Dirac solutions simplify to analytic and anti-analytic functions. To write K in complex notation when dimension is n = 2, we replace points x and y by complex numbers z and w. In the integrand, ∇Φ(y −x) becomes (2π(w−z)) −1 and n(y)dσ(y) becomes −idw. With the expression z · w = Re(zw) for the real inner product, we obtain (1) Kf Here we recognize the Cauchy integral from complex analysis, and we define the principal value Cauchy integral Superficially, E looks not much different from K. Clearly, boundedness of E on a given curve implies boundedness of K on that curve. However, we are more concerned with invertibility of I ± K, and in this case E is a much simpler object than K. In fact E ± = 1 2 (I ± E) are projections, although in general not orthogonal. Here E + projects onto the interior Hardy space: the subspace consisting of traces of analytic functions in Ω. The null space of E + , the subspace along which it projects, is the exterior Hardy space consisting of traces of analytic functions in Ω − which vanishes at ∞. This explains the structure of E itself. It is a reflection operator since E = E + − E − , which reflects the exterior Hardy space E − L 2 = R(E − ) across the interior Hardy space E + L 2 = R(E + ). In particular we see that the spectrum is σ(E) = {+1, −1}.
Hiding behind K there is also a second reflection operator, namely pointwise complex conjugation Nf (z) = f (z), z ∈ ∂Ω, which comes along with its two spectral projections N + f = Re f and N − f = Im f . Note that although we use complex numbers, we will still regard them mainly as vectors. In particular, we consider the two basic reflecion operators E and N as real linear operators. In terms of these projections, the operator 1 2 (I + K) used for solving the Dirichlet problem, is the composition of two (restrictions of) projections, namely the interior Cauchy projection restricted to the interior Hardy subspace. This is readily seen from (1). Equivalently, 1 2 (1 + K) is the compression N + E + N + of the Cauchy projection E + to the subspace N + L 2 of real functions.
The factorization into (3) and (4) explains the relation between the boundary value problem and K. However, we need to re-interpret the Dirichlet problem as a boundary value problem for analytic functions: We regard harmonic functions as real parts of analytic function, neglecting some possible minor topological obstructions. Having switched in this way from the Laplace equation to the Cauchy-Riemann system, the Dirichlet problem now amounts to the Hilbert problem of finding the analytic function in Ω which has a prescibed real part at the boundary. Thus, in terms of operators, solving the boundary value problem means inverting the map (4).
As for the right factor (3), this should be viewed as an ansatz for the solution (or more precisely the trace of the analytic function). The reason d'être for (3) is that it replaces the non-locally defined Hardy space E + L 2 by the locally defined space N + L 2 , which of course is crucial for numerical computations. Ideally, we would like (3) to map N + L 2 bijectively onto E + L 2 , with a not to large condition number. Recall that the condition number κ(T ) = T T −1 of a linear operator is a measure of how easy T is to invert numerically. Unfortunately, the ansatz (3) may not be injective. Moreover, we show in Section 6 that even the Fredholm condition number (3) is comparable to 1/θ as θ → 0, when ∂Ω contains a corner with angle θ. This demontrates that numerically, the integral equation (I + K)h = g may be much worse than the boundary value problem N + h = g, h ∈ E + L 2 , that it is used for solving.
In this paper, we propose a new integral equation for solving the same boundary value problem, which arise upon modifying the maps (3) and (4) to respectively. Our motivation is a non-selfadjoint local elliptic boundary value problem for Dirac equations Details concerning the Dirac operator D is found in Section 2. At the boundary S − : f → 1 2 (1 − ν)f acts from the left by Clifford multiplication and yields a projection. We follow conventions from physics that Clifford squares of real vectors are positive. Stokes' theorem shows that the eigenvalue problem Df = ikf , νf = f only give spectrum in Im k < 0. This gives an indication of that the boundary value problem (7) has good properties. Indeed (7) is as well posed as one can possibly hope for when studying time-harmonic waves with wave number Im k ≥ 0. See Section 4 for more details. Replacing N by the reflection operator S : f → νf above, yields an ansatz (5) with better properties than (3). Due to its origin, we will refer to this new ansatz as the spin ansatz. However, to be able to use the spin ansatz, we need to embed our differential equation into a Dirac equation. This means that we extend the operators N and E to act in a larger spaceL 2 of multivector valued functions, in which the map (4) encoding the original boundary value problem is the restriction to an invariant subspace of a larger map (6) encoding a Dirac boundary value problem. In two dimensions, the Dirac equation amounts to a pair of functions in Ω, one which is analytic and one which is anti-analytic. The resulting equation is as follows. Given Dirichlet datum g : ∂Ω → R, we solve the real linear singular integral equation for h : ∂Ω → C, where |dw| and t(z) denote the scalar measure and unit tangent vector along ∂Ω respectively. The solution u + iv in Ω, u solving the Dirichlet problem and v its harmonic conjugate function, is then the Cauchy integral of h.
The derivation of this spin integral formulation for the Dirichlet problem is found in Example 4.5. For the idea above to work, it is important that the Dirac boundary value problem (6) is well posed just like the original boundary value problem (4) which we embed. In the static case k = 0 in two dimensions above, both the original boundary value problem and the Dirac boundary value problem, are in general well posed only in Fredholm sense. For time-harmonic Maxwell's equations however, the wave number is k = 0, and in this case we will have a well posed boundary value problem with a unique solution in a connected exterior domain. We embed Maxwell's equations into a larger Dirac equation, and also the Dirac boundary value problem will be well posed. Proceeding similar to above in three dimensions, and with k = 0, we show in Example 4.6 how the Maxwell scattering problem in R 3 \ Ω against a perfect conductor Ω, can be solved by a singular integral equation of the form where M is a multiplier involving ν and where Ψ k = ∇Φ modulo weakly singular terms. The auxiliary function h on ∂Ω in this case has four scalar components, but there are no algebraic, differential or integral constraints imposed on h. It follows from Section 4 that this equation is uniquely solvable with a Fredholm condition number (that is using the Calkin algebra norm for operators) comparable to that for the Maxwell boundary value problem. In Section 5 we generalise (9) beyond the case of a perfect conductor, and formulate an integral equation with a spin ansatz for solving the Maxwell scattering problem against a finite number of objects with different scalar och constant permittivity, permeability and conductivity.
As is often the case in research, both before and after Newton, one is sitting on the shoulders of giants. In my case, many of the ideas underlying this paper is a heritage from my PhD supervisor and colleague Alan McIntosh. As already mentioned, the L 2 boundedness of singular integrals of the kind employed in this paper on Lipschitz surfaces follows from the work of Coifman, McIntosh and Meyer [8], by Calderón's method of rotation from the one dimensional case. A direct proof in R n using Clifford algebra was given by McIntosh [14]. In my thesis [2,3,1,6,4], I elaborated on the ideas of McIntosh to solve Maxwell's equations by embedding into the elliptic Dirac equation, and this build on the earlier works of McIntosh and Mitrea [15] and Grognard, Hogan and McIntosh [5]. Although the idea is older than so, the understanding of a boundary value problem in terms of subspaces like in Figure 1 is for me a heritage from McIntosh. In the study of smooth boundary value problems for Dirac operators, see for example the book [7] by Booß-Bavnbek and Wojciechowski, it is standard to formulate boundary conditions in terms of subspaces. But less so in the study of non-smooth boundary value problems for Maxwell's equations or second order elliptic equations.
It is surprising that it is still a somewhat open problem to find a numerically well behaved boundary integral equation for solving scattering problems for Maxwell's equation. See Epstein and Greengard [11] and Epstein, Greengard and O'Neil [12] for recent new Debye formulations, and Colton and Kress [9,10] for the classical formulations. The spin integral equations that we propose in the present paper are based on the McIntosh singular integral approach with Clifford algebra from [15,5,2]. However, the integral formulations there are not suitable for numerical computations, since they suffer from spurious resonances and the same problems as the classical double layer potential equation. The work in the present paper began in [16], where an equation in the spirit of (9) was formulated for solving the Maxwell scattering problem against a perfect conductor, using the formalism from [3]. Both (9) and the spin integral equation from [16] have the advantages of not suffering from spurious resonances and having an improved condition number, at least in the Fredholm sense for Lipschitz boundaries, compared to classical formulations. However, the equation in [16] is for eight unknown scalar functions, as compared to the four unknown scalar functions for the equation in this paper. In both cases, the main novelty lies is the use of an auxiliary spin boundary condition 7, to obtain a singular integral operator with improved condition number. To our knowledge, this local non-selfadjoint Dirac boundary condition has not been exploited in this way before, with numerical computations in mind.

Higher dimension algebra
In this section, we fix notation and survey the higher dimensional algebra which we need for Dirac equations. See [2,16] for more details.
In particular we recall in Example 2.2 how Maxwell's equations fit into this framework. We denote by Ω + a bounded domain in R n with a strongly Lipschitz boundary Σ = ∂Ω + . This means that locally around each point on the boundary, Σ coincides with the graph of a Lipschitz regular function, suitably rotated. The unbounded exterior domain we denote by Ω − = R n \ Ω + . Sometimes we abbreviate Ω = Ω + . The unit normal vector field on Σ pointing into Ω − we denote by ν. The R n standard basis is denoted {e j } n j=1 . Functions to be used are defined on subsets of R n , which take values in the complex exterior algebra for R n , which we denote by ∧C n . This is the 2 n -dimensional complex linear space spanned by basis multivectors The scalars C = ∧ 0 C n ⊂ ∧C n is the one-dimensional subspace corresponding to j = 0, and the vectors C n = ∧ 1 C n ⊂ ∧C n is the n-dimensional subspace corresponding to j = 1. General objects in the exterior algebra we refer to as multivectors. For a given 0 ≤ j ≤ n, we denote the subspace of j-vectors by ∧ j C n , and ∧ ev C n := ⊕ j ∧ 2j C n and ∧ od C n := ⊕ j ∧ 2j+1 C n denotes the even subalgebra and odd subspace of the exterior algebra. We use the hermitean inner product (·, ·) on ∧C n for which the above basis multivectors is an ON-basis. The function space on Σ where we consider our integral equations, is the space On ∧C n we use three complex bilinear products: The exterior product u ∧ w, the (left) interior product u w, and the Clifford product uw. In the special case when the left factor is a vector u ∈ ∧ 1 C n , these products are given by the following basis formulas. For s 1 < s 2 < . . . < s j , let s = {s 1 , . . . , s j } and e s = e s 1 ∧ e s 2 ∧ . . . ∧ e s j , and denote by ǫ(i, s) the number of indices s k ∈ s which are < i. Then The exterior and Clifford products are the associative complex algebra products with 1 = e ∅ ∈ ∧ 0 C n as identity, which are uniquely determined by these basis formulas. The definition of the interior product, which is not associative, for two general multivectors is by duality: We require whenever u has real coordinates. Two useful unitary operations on ∧C n , which are automorphisms with respect to all three products above, are the involution w and reversion w, given by for w = e s = e s 1 ∧ e s 2 ∧ . . . ∧ e s j . We also make use of coordinate-wise complex conjugation which we denote by w c . For a ∈ ∧ 1 C n and w ∈ ∧C n , the Clifford, interior and exterior products are related by aw = a w + a ∧ w and conversely by the Riesz formulas Similar to the discussion in the introduction, we now aim to define two pairs of complementary L 2 projections, using the above algebra. We define Theese yield projection operators onto the subspaces of multivector fields f which are tangential and normal pointwise all over Σ respectively. To explain this algebra, construct at a given point x ∈ Σ an ON-basis with ν(x) as first basis vector. Then write a given multivector w in the induced basis for ∧C n . It is normal if all non-zero terms contain a factor ν. It is tangential if no non-zero terms contain a factor ν. When computing ν ∧ w, a factor ν is added to all tangential terms, and all normal terms are nulled. When computing ν w, a factor ν is removed from all normal terms, and all tangential terms are nulled. The reflection operator N = N + − N − , which reflects normal multivectors across tangential multivectors, can be written The expression (14) is readily seen to be correct by writing f (x) in an ON-basis adapted to ν(x), and (anti-)commute one of the factors ν through f . The boundary value problems we aim at, use N for the description of the boundary conditions. We now turn to the differential equation in the domains, which generalizes the Cauchy-Riemann system: where Df = n j=1 e j (∂ j f (x)) is a Dirac operator. The wave number k ∈ C is always assumed to satisfy Im k ≥ 0, where our main interest k ∈ R \ {0} corresponds to undampened time-harmonic waves. From the factorization (D+ik)(D−ik) = ∆+k 2 of the Helmholtz operator we see that when f solves Df = ikf , then each coordinate function solves the Helmholtz equation. Similar to the Cauchy integral formula for analytic functions, we have a reproducing formula for solutions to Df = ikf . For this we need fundamental solutions: For the Helmholtz operator ∆ + k 2 we use the fundamental solution In R n , the Cauchy singular integral operator on Σ is Note that the kernel uses two Clifford products. Analogous to the static classical two dimensional case, this operator E k is a reflection operator. The spectral projection E + k = 1 2 (I + E k ) projects onto the subspace consisting of traces of solutions to projects onto the subspace consisting of traces of solutions to Df = ikf in Ω − which satisfy a Dirac radiation condition. See [16].
After these definitions, we formulate boundary value problems in this Dirac framework. We can restrict the projections N + and N − to one of the two subspaces E + k L 2 and E − k L 2 . This gives the four maps These represents four different boundary value problems. In (17) we look for a solution in Ω + with a prescribed tangential part on Σ. In (18) we look for a solution in Ω + with a prescribed normal part on Σ. In (19) we look for a solution in Ω − with a prescribed tangential part on Σ. In (20) we look for a solution in Ω − with a prescribed normal part on Σ.
Conversely, we can restrict the projections E + k and E − k to one of the two subspaces N + L 2 and N − L 2 . This gives the four maps These do not represent boundary value problems, but rather are ansatzes for (traces of) solutions to the Dirac equation in the domains. They can be combined with the above boundary value maps respectively to yield integral equations similar to the classical double layer potential equation.
Example 2.1. We describe in two dimensions n = 2, how the Dirac framework in this section is related to the complex analysis used in the Introduction. For complex analysis we do not regard the imaginary unit algebraically as i = √ −1 but rather geometrically as the unit bivector j = e 1 ∧ e 2 . Note that with Clifford algebra j 2 = e 1 e 2 e 1 e 2 = −1 and that reversion (11) gives complex conjugation with respect to j. In this way, a real multivector can be viewed as a pair of complex numbers (z, w), More generally, a complex multivector corresponds to a pair of bicomplex numbers (z, w) where the components z 1 , z 2 , w 1 , w 2 are complex numbers with the imaginary unit i. This is only needed when considering differential equations with non-zero wave number k = 0, like for Maxwell below. For real multivectors and k = 0 as in the introduction, we note that the Dirac operator acts as on a pair of complex valued functions f and g, representing a multivector field f + e 1 g : Ω → ∧R 2 . Here ∂ = ∂ 1 − j∂ 2 and ∂ = ∂ 1 + j∂ 2 . Therefore a solution to D(f + e 1 g) = 0 is a complex analytic function f and an anti-analytic function g, and the Cauchy singular integral operator E = E 0 from (16) acts as the classical Cauchy integral (2) on f , with i replaced by j, and as its anti-analytic analogue on g. For the boundary conditions, we note that whereν = e 1 ν is the unit normal vector ν represented as a complex number. Thus we see that N + (f +e 1 g) gives the real part of f and the tangential part of g, whereas N − (f +e 1 g) gives the imaginary part of f and the normal part of g, on the boundary ∂Ω.
Example 2.2. We are now in three dimensions n = 3. Consider Maxwell's equations in Ω − consisting of a uniform, isotropic and conducting material, so that electric permittivity ǫ, magnetic permeability µ and conductivity σ are constant and scalar. We study electromagnetic wave propagation in Ω − , with material constants ǫ, µ, σ, and Ω + is assumed to be a perfect conductor. Maxwell's equations in Ω − then take the form for the (rescaled) electric and magnetic fields E and H, where the wave number k satisfies k 2 = (ǫ + iσ/ω)µω 2 . The data in the scattering problem we seek to solve are incoming electric and magnetic fields E 0 and H 0 solving (38) in Ω − . Our problem is to solve for the scattered electric and magnetic fields E and H solving (38) in Ω − , an inhomogeneous boundary condition at Σ and the Silver-Müller radiation condition at infinity. At Σ, since the fields vanish in the perfect conductor Ω + , we have boundary conditions with incoming fields E 0 and H 0 from Ω − and surface charges and currents ρ s and J s .
In R 3 , a multivector w can be viewed a collection of two scalars α and β, and two vectors a and b, where w = α + a + * b + * β ∈ ∧C 3 . Here * denotes the Hodge star defined by * 1 = e 1 ∧ e 2 ∧ e 3 , * e 1 = e 2 ∧ e 2 , * e 2 = −e 1 ∧ e 3 and * e 3 = e 1 ∧ e 2 , which identifies ∧ 0 C 3 and ∧ 3 C 3 , and ∧ 1 C 3 and ∧ 2 C 3 respectively. Following the setup in [16], we write the full electromagnetic field as the multivector field F = E + * H, and note that Maxwell's equations implies the Dirac equation 15 for this F . However, F is not a general solution to 15, but satisfies the constraint that the ∧ 0 C 3 and ∧ 3 C 3 parts of F vanishes.
The Maxwell boundary conditions 26 means that ν × E and ν · H are prescribed on Σ. In terms of F , this means that N + F is prescribed. In [16,Sec. 5], it was shown that the constraint ∇ T × E 0 T = ikH 0 N , which the incoming fields will satisfy, will imply that the Dirac solution F to the boundary value problem with g = −N + (E 0 + * H 0 )| Σ , indeed is a Maxwell field in the sense that the F = E + * H for two vector fields E and H solving the Maxwell boundary value problem.
Example 2.3. Consider k = 0 and the reflection operators E = E 0 and N which we use to encode Dirac boundary conditions. In the Introduction, we saw in two dimensions n = 2, that the double layer potential K equals the compression of E to the subspace N + L 2 . More precisely, following Example 2.1, we mean that E and N are restricted to the subspace of complex valued, that is ∧ 0 R 2 ⊕ ∧ 2 R 2 -valued, functions, where E acts by the classical Cauchy singular integral and N acts by complex conjugation. We now explain how K also fits into the Dirac framework in this section, in dimensions n ≥ 3. With the projections N + and N − , it is not possible to compress E to ∧ 0 C n when n ≥ 3. Nevertheless, both K and its adjoint appear in different invariant subspaces for the operators N ± EN ± . For K * , we note that when f is scalar, so that νf is a normal vector field, then So the subspace of normal vector fields is invariant under N − EN − , and its action there is given by −K * upon identifying scalars and normal vector fields. For K, we compute for a scalar function f ∈ H 1 (Σ), that where ∇ T denotes tangential gradient. So the subspace of tangential gradient vector fields is invariant under N + EN + , and its action there is given by K on the potential.

Known well posedness results
In this section we survey the known invertibility results from [15,5,2] for the maps (17) -(24) on the space L 2 (Σ; ∧C n ) on bounded strongly Lipschitz surfaces in R n . We include the proofs since they serve as background later in Section 4. Proof. We demonstrate this for the map (19); the proofs for the other three are similar. For the proof we need the Riesz formula valid for any f ∈ ∧C n and vector ν ∈ ∧ 1 C n . This is the reversed version of (12). The strong Lipschitz and compactness assumptions on Σ shows the existence of a smooth and compactly supported vector field θ such that inf Σ (θ, ν) > 0. Using this we calculate for The identity θν + νθ = 2(θ, ν) is a special case of (13). The last identity above is an application of Stokes' theorem. This leads to the norm estimate This proves the claim since f ∧ ν = N + f and since the Cauchy integral acts as a compact operator L 2 (Σ) → L 2 (U). Proof. We show how the lower estimate Assume g ∈ N − L 2 and apply the above lower bound to f : Since N + g = 0, this implies the claimed lower bound g E + k g modulo a compact term. Proof. By the method of continuity it suffices to consider the case k = 0, since k → E k is a continuous map. It is in fact an operator-valued analytic map. Writing E = E 0 , we note the dualities (Nf, νg) = −(f, ν(Ng)) and (Ef, νg) = −(f, ν(Eg)), for any f, g ∈ L 2 (Σ; ∧C n ). Now consider two of the restricted projections which have the same null space, for example (19) and (22). Computing for f ∈ E − L 2 and g ∈ N − L 2 that (N + f, νg) = (f, νg) = (f, ν(E + g)), it follows that the map dual to (19) is similar to (22). Since (19) and (22) have the same finite dimensional null space, their index must be zero. Proof. For a solution f to Df = ikf in Ω + , we apply Stokes' theorem and obtain For a solution f to Df = ikf in Ω − , we apply Stokes' theorem and obtain , then (f, νf ) = 0 on Σ. If Im k > 0, this forces f = 0. since in this case f decays exponentially at ∞. When k ∈ R \ {0} we instead conclude that |x|=R (f, x |x| f )dσ = 0 for all large R. Next we note the identity 2|f | 2 = |( x |x| − 1)f | 2 + 2( x |x| f, f ). Integrating this over the sphere |x| = R, we obtain lim R→∞ |r|=R |f | 2 dσ = 0 using the Dirac radiation condition satisfied by f ∈ E − k L 2 for the term ( x |x| − 1)f . Rellich's lemma shows that f = 0 since Ω − is connected and ∆f + k 2 f = 0. See [15,16] for more details.

The spin ansatz and new integral equations
In this section we construct new ansatzes for the boundary value problems (17)-(20). Instead of using the ansatzes (21)-(24), we replace the reflection operator N by the reflection operator x ∈ Σ.
To explain why E k together with boundary conditions S should yield boundary value problems with better solvability properties than E k and N, we need to explain some operator algebra. In this we follow [3] and consider the operator algebra generated by two reflection operators A and B, that is A 2 = B 2 = I, abstractly on a Hilbert space H. We think of B as encoding the differential equation through the abstract Hardy subspace projections B ± = 1 2 (I ± B), and of A as encoding boundary conditions through the two complementary projections A ± = 1 2 (I ± A). The most important operators to describe the geometry between A and B are the cosine operator C = 1 2 (AB + BA) and the rotation operators AB and BA. Note that (BA) −1 = AB. Note that (iv) is symmetric under swapping A and B, and therefore so is (i), (ii) and (iii).
Proof. We have identities  Consider now B = E k and A = S. In this case we note that the rotation operator BA is given by x ∈ Σ, since ν 2 = 1. This is really the core observation of this paper. For k = 0, we note that this yields a skew-symmetric operator (ES) * = −ES, with purely imaginary spectrum. In particular it stays well away from ±1, and therefore it is clear that all boundary value problems described by E and S are well posed. Note that this follows from Proposiition 4.1 by abstract arguments and is not using the strong Lipschitz assumption on Σ, in contrast to well posedness for the pair E and N in Section 3. For non-zero k, we have at least that (E k S) * = −E k S modulo compact operators. In this way, we see from abstract considerations only, that all restricted projections are Fredholm operators with index zero. We summarize and complement with injectivity results in the following proposition. Proof. The bounds of the (Fredholm) inverses follows from the skew-adjointness of ES by the formulas in the proof of Proposition 4.1. The fact that the index is zero follows from the method of continuity for k = 0, since k → E k is continuous.
The injectivity results that remains to prove are that The idea of proof is similar to Proposition 3.4. For a solution f to Df = ikf in Ω + , we apply Stokes' theorem and obtain If Im k ≥ 0, this forces f = 0 on Σ and therefore by the Cauchy formula also f = 0 in Ω + .
For f ∈ E − k L 2 ∩ S − L 2 we have (f, νf ) = −|f | 2 on Σ, and a similar application of Stokes' theorem yields If Im k > 0, this forces f = 0 as before, letting R → ∞. since in this case f decays exponentially at ∞. Also when k = 0, f has enough decay for us to conclude. (However, the case k = 0 is already taken care of by the skew-adjointness of ES.) When k ∈ R \ {0} we instead conclude that |x|=R (f, x |x| f )dσ ≤ 0 for all large R. Integrating the identity x |x| f, f ) over the sphere |x| = R, we obtain lim R→∞ |r|=R |f | 2 dσ = 0 since the Dirac radiation condition for f at infinity shows the the first term on the right vanishes at infinity. As in Proposition 3.4, this forces f = 0 by Rellich's lemma if Ω − is connected. If Ω − have bounded connected components, we can argue as for E + k L 2 ∩ S + L 2 = {0} to conclude that f = 0 also in these components. Proposition 4.3 is the main result that we need to obtain the announced spin integral equations for solving Dirac boundary value problems. The idea is to use the ansatz 32, which is always invertible by Proposition 4.3, for the interior boundary value problems (17) and (18). This leads to the integral equations Similarly, for the exterior boundary valur problems (19) and (20), we use the ansatz (33), which is also always invertible by Proposition 4.3. This gives integral equations  (20) respectively, are well posed, and both the domains and ranges are simple pointwise defined subspaces of L 2 (Σ). We can however improve these integral equations a little further for numerical implementation, so that both the domains and ranges are the same subspace, and not depending on ν(x) like S ± L 2 and N ± L 2 do. To this end, we apply again the abstract setup for boundary value problems described in this section. Consider the reflection operator given by pointwise involution of the multivector field. The corresponding spectral subspaces are Computing the relevant cosine operators, we have The proof of Proposition 4.1 give us explicitly invertible maps between subspaces S ± L 2 and subspaces N ± L 2 on the one hand, and explicitly invertible maps between subspaces S ± L 2 and subspaces T ± L 2 on the other hand. Indeed, we note that if A and B are two reflection operators on a Hilbert space H satisfying AB + BA = 0, then the associated eight restricted projections are pairwise inverse, up to a factor 2, as follows.
We can now formulate the main result of this paper, namely spin integral equations for solving the boundary value problems for the differential equation Df = ikf with prescribed tangential or normal part of the field at the boundary. Let Ω + ⊂ R n be a bounded strongly Lipschitz domain, with exterior domain Ω − , and consider a wave number Im k ≥ 0.
• The interior boundary value problem to find a solution f to Df = ikf in Ω + with prescribed tangential/normal part N ± f = g at Σ is well posed in the sense that N ± : E + k L 2 → N ± L 2 is invertible, if and only if the singular integral equation In this case the solution to the boundary value problem is f = E + k S − h at Σ. • The exterior boundary value problem to find a solution f to Df = ikf in Ω ± with prescribed tangential/normal part N ± f = g at Σ is well posed in the sense that N ± : E − k L 2 → N ± L 2 is invertible, if and only if the singular integral equation In this case the solution to the boundary value problem is f = E − k S + h at Σ.
Proof. For the interior boundary value problems, the ansatz E + k : S − L 2 → E + k L 2 is an invertible map for any Im k ≥ 0 by Proposition 4.3. For the exterior boundary value problems, the ansatz E − k : S + L 2 → E − k L 2 is an invertible map for any Im k ≥ 0 by Proposition 4.3. We have also seen that T + S ± : N ± L 2 → T + L 2 and S ± : T + L 2 → S ± L 2 are invertible maps. These invertible maps enable us to fomulate the boundary value problems as singular integral equations on the subspace L 2 (Σ; ∧ ev C n ) as stated.
Example 4.5. We saw in Example 2.1 how the Dirichlet problem for the Laplacian in Ω + ⊂ R 2 , or equivalently the Hilbert boundary value problem for analytic functions with prescribed real part on Σ, can be formulated in terms of invertibility of Theorem 4.4 allow us to solve this boundary value problem as a real linear singular integral equation in the space L 2 (Σ; C) as follows. Given the real valued Dirichlet datum g ∈ L 2 (Σ; R), we compute T + S − g = 1 2 g. To see that the equation T + S − N + E + S − h = T + S − g reduces to (8) using complex algebra, we note that T + S − N + S − h = 1 4 h since S anti-commute with T and N. Writing the Cauchy integral with complex algebra, we have Computing 4T + S − N + (ES − h), we obtain (8) with i replaced by j.
In the formula f = E + S − h for the solution to the boundary value problem, we need only to evaluate the ∧ ev R 2 part of f : The auxiliary anti-analytic function given by the ∧ 1 R 2 part will be trivial due to our choice of g. Thus we end up with the classical Cauchy integral of h for the solution u + jv.
Example 4.6. We saw in Example 2.2 how the Maxwell scattering problem in Ω − ⊂ R 3 against a perfect conductor Ω + can be formulated in terms of invertibility of N + : E − k L 2 → N + L 2 . Theorem 4.4 allows us to solve this linear equation as a singular integral equation in the space L 2 (Σ; C 4 ) as follows. Given the incoming electric and magnetic fields E 0 and H 0 , we compute the tangential Dirac data g = −N + (E 0 + * H 0 )| Σ . Note that g depends on the tangential part of E 0 and the normal part of H 0 . Compute the bivector field It is straightforward to compute that T + S + N + E − k S + h = T + S + g amounts to solving the singular integral equation for h ∈ L 2 (Σ; ∧ ev C 3 ), where M is the multiplier The solution to the Maxwell scattering problem is then This is the algorithm that we propose in this paper for solving the Maxwell scattering problem against a perfect conductor.

Maxwell scattering in piecewise constant media
We formulated in Example 4.6 a spin integral equation for solving the Maxwell scattering problem against a perfect conductor, which is a singular integral equation for four scalar functions. In this section, we fomulate a similar spin integral equation for solving more general scattering problems, for time-harmonic Maxwell's equations at frequency ω. We do not aim to present a complete solvability theory in this section, since it requires a solution of fundamental open problems. Instead we fomulate the algorithm and describe the future work that is needed.
We denote by N the number of bounded materials, and write Here Ω j are assumed to be bounded open sets, with Lipschitz regular boundaries ∂Ω j ⊂ Σ, j = 1, . . . , N, and Ω 0 is the complement of a bounded Lipschitz domain. The Lipschitz interface Σ is Σ = ∂Ω 1 ∪ . . . ∪ ∂Ω N . Write Σ i,j = Ω i ∩ Ω j . A unit normal vector at a boundary point x ∈ Σ, which is well defined almost everywhere, is denoted n = n(x). By n j = n j (x) at x ∈ ∂Ω j we mean the unit normal vector which is outward pointing relative Ω j . The region Ω j , j = 0, 1, . . . , N, we assume represent a homogeneous, linear and isotropic material, with electric permittivity ǫ j , magnetic permeability µ j and conductivity σ j as constant and scalar quantities. We formulate Maxwell's equations in Ω j as where the wave number is k j = ωα j /β j , with α j = (ǫ j + iσ j /ω) 1/2 and β j = µ −1/2 j . We choose to normalise the fields so that E denotes the geometric mean of the electric field and displacement, and H denotes the geometric mean of the magnetic field and intensity, so that the square of the fields has energy density as dimension.
In particular this means that at the interface Σ, we have jump conditions which require continuity of n · (β −1 H), n × (α −1 E), n × (βH), and n · (αE) if Maxwell's equations for the original electric and magnetic fields are to hold in distributional sense in all R 3 .
The data in the scattering problem we seek to solve are incoming electric and magnetic fields E inc and H inc solving (38) with k = k 0 in Ω 0 . Our problem is to solve for electric and magnetic fields E j and H j • solving Maxwell's equations (38) with wave number k j in Ω j , j = 0, 1, . . . , N, • with E 0 , E 0 satisfying the Silver-Müller radiation condition at infinity, see [16, eq. (4)], • and where E inc + E 0 + . . . + E N and H inc + H 0 + . . . + H N solve Maxwell's equations in distributional sense across Σ.
We note as before that in terms of the electromagnetic multivector field F = E + * H, we can write Maxwell's equations (38) as the Dirac equation in Ω j , and well posedness in L 2 = L 2 (Σ; ∧C 3 ) of the scattering problem descibed above follows from invertibility of the map Here E j L 2 denotes the image of E + k j L 2 (∂Ω j ) under the inclusion L 2 (∂Ω j ) ⊂ L 2 (Σ), where N j denotes the map and N j f (x) = 0 when x ∈ Σ \ ∂Ω j . Recall that T + is projection onto L 2 (Σ; ∧ ev C 3 ) and T − is projection onto L 2 (Σ; ∧ od C 3 ). We use the spin ansatz With this setup we obtain the following spin integral equation for solving the above Maxwell scattering problem.
Algorithm 5.1. Let E inc and H inc be the incoming rescaled electric and magnetic fields in Ω 0 , and define Solve the singular integral equation Then the solution to the Maxwell scattering problem in this section is given by Since S j are complementary projections in L 2 (Σ), it follows from Proposition 4.3 that the spin ansatz S Σ is an invertible map for any Im k ≥ 0. As we discussed in the introduction, for the spin integral equation to be computationally useful we also need to show that the Dirac scattering problem which we embed the Maxwell scattering problem into, is well posed. We conjecture that this is the case for Algorithm 5.1. It appears though that even for one bounded material, N = 1, this is beyond the currently available L 2 (Σ) solvability techniques, which are based on Rellich estimates like in Proposition 3.1, if we allow general Lipschitz interfaces. To show the main problem, consider the jump relations on Σ = ∂Ω 1 = ∂Ω 0 when N = 1. We want to show that (f 1 , f 0 ) → g is a Fredholm map in the L 2 (Σ) topology. To this end we note that where E = E 0 and K = (−1) j 1 2 T − (E k j − E) is a compact operator. Adding α + α − times (40) and (42) yields the estimate with a compact operator K ′ . This yields a Fredholm bound on T − f 1 + T − f 0 , provided λ = (α + + α − )/(α + − α − ) is outside the Fredholm spectrum of the rotation operator EN. Similarly, by adding β + β − times (39) and (41), we obtain a Fredholm bound on T + f 1 + T + f 0 , provided (β + + β − )/(β + − β − ) is outside the Fredholm spectrum of the rotation operator EN. As shown in [3], these conditions are equivalent to that (α + /α − ) 2 and (β + /β − ) 2 maps into the Fredholm resolvent set for the cosine operator 1 2 (EN + NE) by the map z → (z + 1)/(z − 1). Thus, if we allow arbitrary conductivity σ ≥ 0 and permittivity ǫ > 0, we need to know that the spectral radius of the cosine operator is ≤ 1. As we have seen in Example 2.3 and Proposition 4.1, the classical double layer potential operator embeds into this cosine operator, so in particular we need to know that the spectral radius of K is at most one. This well known spectral radius conjecture is to the authors knowledge still an open problem for general Lipschitz surfaces.
However, if furthermore Σ is smooth, then the spectrum of (EN + NE)/2 is contained in the unit disk. Indeed modulo compact operators is ν is smooth, since ab + ba = 2(a, b) with Clifford algebra. Here the double layer potential K acts componentwise on the multivector field f . Since E is a reflection operator it follows that E is a compact perturbation of a unitary operator. In particular we deduce that the Fredholm spectrum of the cosine operator is contained in [−1, 1].
To show injectivity of B Σ one can generalise the methods in Proposition 3.4. We omit the details and refer to [4].

Problems with the classical ansatz
We end this paper with an explicit computation on a conical domain, elaborating on Mellin transform techniques of Fabes, Jodeit and Lewis [13], that shows that the classical double layer potential equation may have a condition number which is significantly worse than that of the underlying boundary value problem. This in contrast to the spin integral equations proposed in Theorem 4.4, which typically is no worse than the boundary value problem numerically. By localising the result below, we obtain similar results for bounded domains which have corners.
Consider the double layer potential operator K in dimension 2 given by (1) in the Introduction, on a cone Ω = {z ∈ C ; 0 < arg z < θ}.
We are in particular interested in the limit as θ → 0 + . Recall that K is the composition of the restricted Cauchy projection (3) and the restricted real part projection (4), and that all these operators are considered as real linear only. The purpose of this example is to show by explicit computation, that (I + K) −1 ≈ 1/θ 2 , (N + : E + L 2 → N + L 2 ) −1 ≈ 1/θ, and (E + : N + L 2 → E + L 2 ) −1 ≈ 1/θ as θ → 0 + . This means that the interior Dirichlet problem which we intend to solve become increasingly illposed at the rate 1/θ. On the other hand, the operator I + K which is classically used to solve the Dirichlet problem become ill-conditioned at the faster rate 1/θ 2 . Note that the operators K, E and N are themselves uniformly bounded as θ → 0 + .
Since the computation uses the Fourier transform, it is natural to complexify these real linear operators, which we do as follows. The imaginary unit i used in the definition of E and N, we write as j and rather think of as the unit bivector j = e 1 ∧ e 2 as in Example 2.1, which squares to j 2 = −1 with the Clifford product. In the framework from Section 2, this means that N is reversion (rather than complex conjugation) and N + is projection onto scalars ∧ 0 C 2 ≈ C. Our operators act on functions taking values in the even subalgebra ∧ ev C 2 = {z + jw ; z, w ∈ C}, that is the commutative algebra of bicomplex numbers.
To this end, we use that the inverse of a matrix A = a b c d , where a, b, c, d in general are non-commuting operators, is given by Applying this to A equal to the Fourier multiplier of I + EN, it suffices for us to bound a −1 and (d − ca −1 b) −1 . We calculate a −1 = (1 − ij tanh(πξ)N) −1 = 1 + ij tanh(πξ)N 1 + tanh 2 (πξ) , which is bounded and independent of θ. With straightforward bicomplex algebra, and noting that reversion N commutes with i but anti-commutes with j, we also compute d − ca −1 b = 1 + ij tanh(2πξ)N + (cos(α) + j sin α) cosh(2αξ) − ij sinh(2αξ) cosh(2πξ) .
Doing the estimates, this yields |(d − ca −1 b) −1 | ≈ 1/θ, from which we deduce that the norms of the inverses of the restricted projections (3) and (4) are of the order 1/θ in this example.