Abstract
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 a wider range of permittivities than other known formulations, 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, our Dirac integral equation performs equally well in three dimensions, as demonstrated in a companion paper.
Similar content being viewed by others
1 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 main novelty of our BIE is that it performs well for a very wide range of materials, including metamaterials with negative permittivities. We provide not only numerical evidence of such performance, but also rigorous proofs for general Lipschitz interfaces.
We refer to our new formulation as a Dirac integral equation, since our point of departure is to embed Maxwell’s equations into a Dirac-type equation by relaxing the constraints \({{\mathrm{div}}}E=0\) and \({{\mathrm{div}}}B=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. This idea to write Maxwell’s equations as a Dirac equation is old, and appears in the works of M. Riesz and D. Hestenes. Even Maxwell himself formulated his equations using Hamilton’s quaternions. Picard [33, 34] elliptifies Maxwell’s equations by embedding these into a Dirac equation, without using Clifford algebra but writing it as a div-curl-grad matrix system. This Picard’s extended Maxwell system is the basis for more recent works [38,39,40].
In harmonic analysis of non-smooth boundary value problems for Maxwell’s equations, BIEs based on the Dirac equation were developed by M. Mitrea, McIntosh et al., see [5, 31]. Building on this, the second author of the present paper further developed the theory of Dirac integral equations for solving Maxwell’s equations in his PhD work [1,2,3,4, 6] with Alan McIntosh. Applications of these methods to scattering from perfect electric conductors (PEC) are found in [35, 36], and the spin integral equations there are precursors of the Dirac integral equations presented here. More recent results on Dirac equations for Maxwell scattering problems with Lipschitz interfaces are also [26, 30], which deal with the \(L_p\) boundary topology, but only treat the case of equal wave numbers in the two domains.
The focus in the present paper is on BIEs which are numerically good for a very wide range of material combinations, including such where permittivity ratios are negative real numbers. In this case, the relevant function space is the trace space for the \(L_2\) topology in the domains, since this is where surface plasmons waves may appear, as we approach purely negative permittivity ratios. In particular, we seek BIEs without false eigenwavenumbers in this plasmonic regime. In the simpler cases of scattering from PEC objects or from dielectric objects with less permissible restrictions on their permittivities, and also on the genus of the domains, there are many good formulations available. See, for example, [15] for a BIE with six scalar unknowns containing only weakly singular integral operators and extending results from the pioneering work [40]. However, for metamaterials in plasmonic scattering, we are only aware of two BIEs without false eigenwavenumbers, namely the Dirac integral equations presented here and the BIE in [23]. Both these works seem to indicate that eight is the number of scalar densities needed to construct a Maxwell BIE in 3D that is free from false eigenwavenumbers for all complex-valued choices of material constants.
Turning to the description of 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
The PDE problem, Helmholtz or Maxwell, is to invert a map \((F^+,F^-)\mapsto g\), where g is boundary datum, and \(F^\pm \) are the solutions, the sought fields in the interior domain \(\Omega ^+\) and exterior domain \(\Omega ^-\), respectively. To obtain a BIE we need to choose an appropriate ansatz, that is an integral representation \(h\mapsto (F^+,F^-)\) of the fields in terms of a density function h on the interface \(\partial \Omega \) between \(\Omega ^+\) and \(\Omega ^-\). From this we obtain a BIE as the composed map \(h\mapsto 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 \({{\mathbf {R}}}^2\) as well as in \({{\mathbf {R}}}^3\), are formulated in Sect. 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_\pm \) are wavenumbers in \(\Omega ^\pm \), \({\mathrm{Im}}\,k_\pm \ge 0\), and \({\hat{k}}=k_+/k_-\).
-
Propositions 6.1 and 8.2 show sufficient conditions for the Helmholtz and Maxwell transmission problems to have uniqueness of solutions. For non-magnetic materials, our result is the hexagonal region shown in Fig. 1, which strictly contains the regions earlier obtained for the Helmholtz transmission problem in [27, 28]. See, however, the proof of [23, Prop. 3.1] for a comment on a minor flaw in [28]. When also the PDE problem defines a Fredholm map, something that can fail in the hexagon only when \(|\arg (k_+)-\arg (k_-)|=\pi /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 Sect. 2, is only one possible choice from a family of Dirac BIEs derived in Sects. 7 and 8. In \({{\mathbf {R}}}^2\), these BIEs depend on four Dirac parameters \(r,\beta ,\alpha ',\beta '\). 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 \({{\mathbf {R}}}^3\), there are two additional parameters \(\gamma ,\gamma '\).
The Dirac BIE is constructed in Propositions 7.3 and 8.5 using a dual PDE problem as ansatz. Our choice of dual PDE problem giving the Dirac BIEs in Sect. 2, has to effect that these BIEs are invertible whenever the PDE problem is well posed.
-
Propositions 6.2 and 8.3 show that the PDE problem defines a Fredholm map provided that the quotient \(\hat{\epsilon }\) of the permittivities is not in an interval \([-C, -1/C]\) for some \(C\ge 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, \({{\mathcal {H}}}_2\) or \({{\mathcal {H}}}_3\), which roughly speaking is a suitable mix of the function space \(H^{1/2}(\partial \Omega )\subset L_2(\partial \Omega )\) of traces to Sobolev \(H^1(\Omega )\) functions and its dual space \(H^{-1/2}(\partial \Omega )\supset L_2(\partial \Omega )\). In the case of pure metamaterials \(\hat{\epsilon }<0\), the PDE problem to be solved may fail to be a Fredholm map, and hence no integral equation for solving it can be Fredholm. However, when the PDE problem has unique solution we can solve it by an injective Dirac integral equation, and when Fredholmness fails we solve the PDE problem numerically for smooth boundary data by taking a limit from \({\mathrm{Im}}\,\hat{\epsilon }>0\). Numerical experiments indicate the existence of this limit. Theoretically it should be possible to prove its existence by analyzing the problem in slightly larger fractional Sobolev spaces.
More precisely, in \({{\mathbf {R}}}^n\), \(n\ge 2\), one can show that the PDE problem (4) defines a Fredholm map when \((\hat{\epsilon }+1)/(\hat{\epsilon }-1)\) is outside the essential spectrum of the double layer potential operator
$$\begin{aligned} K_\mathrm{d} f(x):=\int _{\partial \Omega }\langle \nu (y),(\nabla \Phi )(y-x) \rangle f(y)d\sigma (y), \qquad x\in \partial \Omega , \end{aligned}$$(1)where \(\Phi \) is the Laplace fundamental solution and \(\nu \) is the outward unit normal. Another basic operator for the Maxwell problem (14) in \({{\mathbf {R}}}^3\) is the magnetic dipole operator
$$\begin{aligned} K_\mathrm{m} f(x):= \nu (x)\times \int _{\partial \Omega }(\nabla \Phi )(y-x)\times f(y)d\sigma (y), \qquad x\in \partial \Omega , \end{aligned}$$(2)acting on tangential vector fields f. Compact perturbations of the operators \(K_\mathrm{d}\) and \(K_\mathrm{m}\), and their adjoints \(K_\mathrm{d}^*\) and \(K_\mathrm{m}^*\) with respect to the standard \(L_2\) pairing on \(\partial \Omega \), appear along the diagonals of our basic Cauchy integral operators (53) and (55). In particular, the (3:4,3:4) and (7:8,7:8) size \(2\times 2\) diagonal blocks in (55) are compact perturbations of \(-K_\mathrm{m}^*\).
It is important to note that in the energy norm, the essential spectra of \(K_\mathrm{d}\) and \(K_\mathrm{m}\) are subsets of \((-1,1)\). If considering instead the \(L_2(\partial \Omega )\) norm, then on domains in \({{\mathbf {R}}}^2\) with one corner, the essential spectrum of \(K_\mathrm{d}\) is a lying “figure eight”, centered around 0. See Fabes, Jodeit and Lewis [14] and the final comments and (137) in the present paper.
-
Although our main concern is wavenumbers \(k_\pm \ne 0\), it is important to have a good behaviour of the BIEs as frequency \(\omega \rightarrow 0\). Typical problems which can occur are dense-mesh low-frequency breakdown and topological low-frequency breakdown. We refer to [13, 41] 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 Sect. 2. One reason for this is that our four densities include the gradient of the field in \({{\mathbf {R}}}^2\) and our eight densities include the electric and the magnetic fields in \({{\mathbf {R}}}^3\), so that no numerical differentiation is needed. Moreover, if \(k_\pm \rightarrow 0\) with \(\hat{k}\) constant in (10) or (20), then \(P,P',N,N'\) are all constant whereas the Cauchy singular integral operators \(E_{k_\pm }\) converge to \(E_0\) in operator norm. That the limit Dirac BIE does not have any false eigenwavenumbers, for fixed \(0<|\hat{k}|<\infty \) as \(k_\pm \rightarrow 0\), is shown in [24, Sec. 7].
We remark that BIEs derived from the Picard system [39, 40] do not suffer from dense-mesh low-frequency breakdown, but are known to have false/near-false eigenwavenumbers. See [41, p. 163] for an excellent survey of such “charge-current formulations”. Note that the Picard equation [39, Eq. (45)] is the Dirac equation (41) written in matrix form, and that the BIEs [39, Eqs. (43),(44)] are \(8\times 8\) systems of the second kind with single and double potential type operator blocks. Thus there is much formal resemblance with the Dirac BIE, but the real difference lies in the precise choice of parameters made in the present paper.
The paper ends with Sect. 9, where the numerical performance of the Dirac integral equations is studied in \({{\mathbf {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 [23, 27] which is only applicable in \({{\mathbf {R}}}^2\), which involves only two unknown scalar densities, and which for certain \(k_\pm \) coincides with a 2D version of the classic Müller system [32, p. 319]. Not surprisingly this special-purpose system performs best, but the new \(4\times 4\) Dirac system is not far behind. In particular it performs well under computationally challenging plasmonic conditions, that is when \({\hat{\epsilon }}\) is negative and real.
In the companion paper [24] by the authors joint with A. Karlsson, we test the new \(8\times 8\) Dirac system for the Maxwell transmission problem in \({{\mathbf {R}}}^3\) using the numerical techniques in [20]. All BIE systems that we are aware of in the literature, with size less than \(8\times 8\), exhibit false eigenwavenumbers – primarily at the left middle corner point in Fig. 1. Our Theorem 2.3 gives very weak sufficient conditions for our Dirac integral equation (20) not to have false eigenwavenumbers. When these conditions fail, notably at the bottom middle corner point in Fig. 1, it is still possible to avoid false eigenwavenumbers with other choices of Dirac parameters in Sect. 8.
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 \(\Omega ^+\), separated by the interface \(\partial \Omega \) from the exterior domain \(\Omega ^-:= {{\mathbf {R}}}^n\setminus \overline{\Omega ^+}\). Fix a large \(R<\infty \) and define 
. For simplicity we assume throughout this paper that \(\Omega ^-\) 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 \(\Omega ^+\) and \(\Omega ^-\) respectively, and a jump parameter \({\hat{\epsilon }}\). We assume that \({\mathrm{Im}}\,k_\pm \ge 0\), but unless otherwise stated we do not assume any relation between \(k_+, k_-\) and \({\hat{\epsilon }}\). We denote the ratio between the wavenumbers by \({\hat{k}}= k_+/k_-\). Our main interest is in scattering for non-magnetic materials where \({\hat{\epsilon }}= {\hat{k}}^2\). In applications, the parameter \({\hat{\epsilon }}\) appears as the ratio \({\hat{\epsilon }}= \epsilon _+/\epsilon _-\) of the permittivities \(\epsilon _\pm \) in \(\Omega ^\pm \). Similarly for magnetic materials, we have a ratio \({\hat{\mu }}= \mu _+/\mu _-\) of the permeabilities \(\mu _\pm \) in \(\Omega ^\pm \). In this general case, we have the relation \({\hat{\epsilon }}= {\hat{k}}^2/{\hat{\mu }}\).
We remark that in what follows, \(E^\pm \) denotes the standard electric fields, but we have rescaled the magnetic fields \(B^\pm \) so that what we here call \(B^\pm \), in standard notation reads \(B^\pm /\sqrt{\epsilon _\pm \mu _\pm }\).
To formulate our results, we need the following conical subsets of \({{\mathbf {C}}}\). We use the argument \(-\pi <\arg (z)\le \pi \).
Definition 2.1
Let \(k_-, k_+\in {{\mathbf {C}}}\setminus \{0\}\) be such that \({\mathrm{Im}}\,k_-\ge 0\) and \({\mathrm{Im}}\,k_+\ge 0\). Define \(\phi _\pm := |\arg (k_\pm /i)|\). Let \({\mathrm{WP}}\,(k_-,k_+)\subset {{\mathbf {C}}}\setminus \{0\}\) be the set of \(z\in {{\mathbf {C}}}\setminus \{0\}\) satisfying
Our notation is to denote fields in the domains by \(F, U,\ldots \), with suitable superscript ±, and to write \(f,u,\ldots \) for respective boundary traces. On \(\partial \Omega \) we denote by \(\nu \) the unit normal vector pointing into the exterior domain \(\Omega ^-\). Furthermore \(\{\nu ,\tau \}\) and \(\{\nu ,\tau ,\theta \}\) denote positive ON-frames on \(\partial \Omega \) depending on dimension, so that \(\tau \) is the counter-clockwise tangent on curves \(\partial \Omega \subset {{\mathbf {R}}}^2\). On surfaces \(\partial \Omega \subset {{\mathbf {R}}}^3\), the theory we develop works for any choice of tangent ON-frame \(\{\tau ,\theta \}\). We denote the directional derivative in direction v by \(\partial _v\), and with slight abuse of notation, \(\partial _\nu u\) denotes normal derivative on \(\partial \Omega \) although this use values of U in a neighbourhood of \(\partial \Omega \).
Our main result for 2D scattering concerns the Helmholtz transmission problem
where \(\Omega ^+\subset {{\mathbf {R}}}^2\), \(n=2\), and with \(u^0\in H^{1/2}(\partial \Omega )\) being the trace of an incoming wave \(U^0\), and we want to solve for \(U^+\in H^1(\Omega ^+)\) and \(U^-\in H^1(\Omega ^-_R)\). We remark that for \({\mathrm{Im}}\,k_->0\), this not so well-known sufficient radiation condition of exponential growth entails a necessary condition of exponential decay. See [37, Prop. 9.3.6] for proofs.
Let \({\mathrm{diag}}\,D\) denote the square diagonal matrix with diagonal D, and let \(A^{\mathrm{T}}\) denote the transpose of a matrix A. Define
Theorem 2.2
Let \(\Omega ^+\subset {{\mathbf {R}}}^2\) be a bounded Lipschitz domain, with \(\Omega ^-\) being connected and notation as above. Then there exists a constant \(1\le C(\partial \Omega )<\infty \) depending on the Lipschitz constants of the parametrizations of \(\Omega ^+\) by smooth domains, so the following holds.
The transmission problem (4) is well posed if
For its solution, consider the Dirac integral equation
for four scalar functions \(h=[h_1\;h_2\;h_3\;h_4]^{\mathrm{T}}\), where \(E_k\) is the singular integral operator (53) which we introduce in Sect. 4, \(P,P', N, N'\) are the constant diagonal matrices (5)–(8), and
The operator in (10) is invertible on the energy trace space \({{\mathcal {H}}}_2\) from (63), introduced in Sect. 7, whenever (9) holds and \({\hat{k}}\in {{\mathbf {C}}}\setminus (-\infty ,0]\). Moreover, the solution to (4) in \(\Omega ^\pm \) is obtained from \(h^+= N' h\) and \(h^-= P'h\) as
and
with notation as in Sect. 4, where \({\tilde{K}}\) and \({\tilde{S}}\) denote layer potentials, \(\nu '\) and \(\tau '\) denote the normal and tangential vector at the point of integration \(y\in \partial \Omega \), and \(I=\begin{bmatrix} 1 &{}\quad 0 \\ 0 &{}\quad 1 \end{bmatrix}\) and \(J=\begin{bmatrix} 0 &{}\quad -1 \\ 1 &{}\quad 0 \end{bmatrix}\).
Proof
This result is derived in Sects. 6 and 7, where we arrive at equation (103) with \(\alpha ={{\hat{\epsilon }}}\). We precondition (103) by multiplying from the left by P and writing \({\tilde{h}}= P' h\) to obtain (10) with \(N=PM\), \(N'={\hat{k}}M'P'\) and \(P({\hat{k}}M'+M)P'=I\). We remark that the parameter ratios \({\hat{k}}<0\) and \({{\hat{\epsilon }}}=-1\) must be excluded for P and \(P'\) to be well-defined. However, the only problem when \({\hat{k}}<0\) is that we cannot precondition (103) with P and \(P'\) to achieve a non-integral term I in (10).
In applications of the Dirac problem (58) to the Helmholtz problem (4), as in Example 7.1, we have \(f^0= ik_- u^0 +\nabla u^0\), that is (11). Equations (12) and (13) for the fields are seen from the first row and the last two rows in (53) respectively, and the expressions for \(h^\pm \) are seen from the right and middle factors in (95). Note that \(h^-\) happens to coincide with the density \({\tilde{h}}\), which is introduced in (98) in a different context. \(\square \)
We remark that it is our experience that the precise choice of P and \(P'\) for Theorem 2.2, as well as for Theorem 2.3, is less important as long as they satisfy \(P({\hat{k}}M'+M)P'=I\). With the choices made, our aim has been to minimize the growth of P, \(P'\), N and \(N'\) as we vary |k|.
We next formulate the analogous result for the Maxwell transmission problem in \({{\mathbf {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^\pm \) and \(B^\pm \). All fields are assumed to be \(L_2\) integrable in a neighbourhood of \(\partial \Omega \).
Define matrices (where we write \({\hat{c}}=1/{\hat{k}}\) and \({\hat{\mu }}= {\hat{k}}^2/{\hat{\epsilon }}\))
Theorem 2.3
Let \(\Omega ^+\subset {{\mathbf {R}}}^3\) be a bounded Lipschitz domain, with \(\Omega ^-\) being connected and notation as above. Then there exists a constant \(1\le C(\partial \Omega )<\infty \) depending on the Lipschitz constants of the parametrizations of \(\Omega ^+\) by smooth domains, so the following holds.
The transmission problem (14) is well posed if
For its solution, consider the Dirac integral equation
for eight scalar functions \(h=[h_1\; h_2\; h_3\; h_4\; h_5\; h_6\; h_7\; h_8]^{\mathrm{T}}\), where \(E_k\) is the singular integral operator (55) which we introduce in Sect. 4, \(P,P', N, N'\) are the constant diagonal matrices (15)–(18), and
with field components in the frame \(\{\nu ,\tau ,\theta \}\). The operator in (20) is invertible on the energy trace space \({{\mathcal {H}}}_3\) from (64), introduced in Sect. 8, whenever (19), \({{\hat{k}}}^2/{\hat{\epsilon \in }} {{\mathbf {C}}}\setminus [-C(\partial \Omega ), -1/C(\partial \Omega )]\) and \({\hat{k}}\in {{\mathbf {C}}}\setminus (-\infty ,0]\) holds. Moreover, the solution to (14) in \(\Omega ^\pm \) is obtained from \(h^+= N' h\) and \(h^-= P'h\) as
and
with notation as in Sect. 4, where \({\tilde{K}}\) and \({\tilde{S}}\) denote layer potentials, \(\nu ',\tau ',\theta '\) denote the frame vectors at the point of integration \(y\in \partial \Omega \), \(I=\begin{bmatrix} 1 &{} 0 \\ 0 &{} 1 \end{bmatrix}\), and \(v\times \cdot \) denotes the map \(x\mapsto v\times x\).
Proof
This result is derived in Sect. 8, where we arrive at Eq. (133) with \(\alpha ={\hat{\epsilon }}\). We precondition (133) by multiplying from the left by P and writing \({\tilde{h}}= P' h\) to obtain (20) with \(N=PM\), \(N'={\hat{k}}^{-1}M'P'\) and \(P({\hat{k}}^{-1}M'+M)P'=I\). We remark that the parameter ratios \({\hat{k}}<0\), \({\hat{\mu }}= {{\hat{k}}}^2/{\hat{\epsilon }}=-1\) and \({{\hat{\epsilon }}}=-1\) must be excluded for P and \(P'\) to be well-defined. However, the only problem when \({\hat{k}}<0\) is that we cannot precondition (133) with P and \(P'\) to achieve a non-integral term I in (20).
In applications of the Dirac problem (58) to the Maxwell problem (14), or in multivector notation (107), we have \(F^0= E^0 +B^0\), that is (21). Equations (22) and (23) for the fields are seen from the second to fourth rows and last three rows in (55) respectively. \(\square \)
3 Multivector Algebra
This section contains the basics of multivector algebra which we need. For a complete account of the theory of multivector algebra and Dirac equations which we use, we refer to [37]. For our purposes, multivectors are the same objects as Cartan’s alternating forms, and multivector fields are the same as differential forms, and amounts to an algebra of not only the one-dimensional vectors but also j-dimensional algebraic objects for \(0\le j\le n\) in \({{\mathbf {R}}}^n\). Concretely, multivectors in \({{\mathbf {R}}}^n\) are the following type of objects, where our interest is in \(n=2,3\). Denote by \(\{e_1,\ldots , e_n\}\) the standard vector basis for \({{\mathbf {R}}}^n\). We write \(\wedge {{\mathbf {R}}}^n\) for the complex \(2^n\) dimensional space of multivectors in \({{\mathbf {R}}}^n\), which is spanned by basis multivectors \(e_s\), where \(s\subset {{\overline{n}}}= \{1,2,\ldots , n\}\). We write \(\wedge ^j{{\mathbf {R}}}^n\) for the subspace spanned by those \(e_s\) with j number of elements in the index set s, so that
We identify \(\wedge ^0{{\mathbf {R}}}^n={{\mathbf {C}}}\) and \(\wedge ^1{{\mathbf {R}}}^n={{\mathbf {C}}}^n\) with the scalars and vectors respectively. Objects in \(\wedge ^2 {{\mathbf {R}}}^n\) are referred to as bivectors.
In \({{\mathbf {R}}}^2\), a multivector is of the form
and so amounts to two scalars a, b and a vector \(v= v_1e_1+v_2e_2\). In \({{\mathbf {R}}}^3\), a multivector is of the form
and so amounts to two scalars a, b and two vectors \(v= v_1e_1+v_2e_2+v_3e_3\) and \(u= u_1e_1+u_2e_2+u_3e_3\) via Hodge duality (see below).
On this \(2^n\) dimensional space \(\wedge {{\mathbf {R}}}^n\) we use three products which, depending on dimension, generalize the classical vector operations: the exterior product \(u\mathbin {\scriptstyle {\wedge }}w\), the (left) interior product \(u\mathbin {\lrcorner }w\) and the Clifford product uw. The last is the standard short hand notation for the Clifford product, whereas the notation \(e_s\mathbin {\scriptscriptstyle {\triangle }}e_t\) was introduced in [37], for reasons which are clear from (29). Products of basis multivectors \(e_s\) correspond to signed unions, set differences and symmetric differences of index sets s respectively. More precisely, for index sets \(s,t\subset {{\overline{n}}}\) we define the permutation sign 
, that is the sign of the permutation that rearranges \(s\cup t\) in increasing order. Then
where \(s\triangle t= (s\cup t)\setminus (s\cap t)\) is the symmetric difference of sets. The exterior and Clifford products are associative, but the interior product is not. On the one hand, orthogonal vectors in \(\wedge ^1{{\mathbf {R}}}^n\) anticommute both with respect to the exterior and Clifford products. On the other hand, \(e_j\mathbin {\scriptstyle {\wedge }}e_j=0\) whereas \(e_je_j=1\), \(j=1,\ldots , n\). An important special case of the interior product is the Hodge star duality
The most important instance of the above three products is when the first factor is a vector, and not a general multivector. In this case, the three products are related as
A multivector \(w\in \wedge {{\mathbf {R}}}^n\) at a point \(x\in \partial \Omega \), with normal vector \(\nu \in \wedge ^1 {{\mathbf {R}}}^n\), is called tangential if \(\nu \mathbin {\lrcorner }w=0\) and is called normal if \(\nu \mathbin {\scriptstyle {\wedge }}w=0\).
Example 3.1
A general multivector in \({{\mathbf {R}}}^2\) is of the form \(w= a_1+v+*a_2\), where \(a_1,a_2\in \wedge ^0{{\mathbf {R}}}^2\) and \(v\in \wedge ^1{{\mathbf {R}}}^2\). The three basic multivector products with a vector \(u\in \wedge ^1{{\mathbf {R}}}^2\) can be written
Note that on vectors \(u\in \wedge ^1{{\mathbf {R}}}^2\), the Hodge star \(*u\) is counter clockwise rotation \(\pi /2\).
Example 3.2
A general multivector in \({{\mathbf {R}}}^3\) is of the form \(w= a_1+v_1+*v_2+*a_2\), where \(a_1,a_2\in \wedge ^0{{\mathbf {R}}}^3\) and \(v_1,v_2\in \wedge ^1{{\mathbf {R}}}^3\). The three basic multivector products with a vector \(u\in \wedge ^1{{\mathbf {R}}}^3\) can be written
We refer to [37, Chapters 2,3] for further details of multivector algebra.
In the same way that the scalar and vector products induce the divergence \(\nabla \cdot F\) and curl \(\nabla \times F\) in vector calculus, the exterior, interior and Clifford products induce the exterior, interior and Clifford/Dirac derivatives in multivector calculus, and we write these as
Replacing u and w by \(\nabla \) and F(x) in Example 3.2, it follows how these differential operators can be written in terms of divergence, curl and gradient in \({{\mathbf {R}}}^3\). We refer to [37, Chapters 7,8] for further details of multivector calculus.
The time-harmonic wave Dirac equation which we consider in this paper is
for multivector fields F, and wave number \(k\in {{\mathbf {C}}}\). The main applications are the Helmholtz and Maxwell equations. Given a scalar solution U to the Helmholtz equation \(\Delta U+ k^2 U=0\), we define the multivector field
It follows that \(F=F_{\text {helm}}\) satisfies (41) since \({{\mathbf {D}}}^2=\Delta \).
For 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 \(\wedge ^0{{\mathbf {R}}}^3\) and \(\wedge ^3{{\mathbf {R}}}^3\) parts are zero. In this formalism, the Dirac equation (41) for \(F=F_{\mathrm{em}}\,\) given by (43) coincides with the time-harmonic Maxwell’s equations, that is the PDE from (14). Indeed, the \(\wedge ^0{{\mathbf {R}}}^3\), \(\wedge ^1{{\mathbf {R}}}^3\), \(\wedge ^2{{\mathbf {R}}}^3\) and \(\wedge ^3{{\mathbf {R}}}^3\) parts of (41) are the Gauss, Maxwell–Ampère, Faraday and magnetic Gauss law respectively. See [37, Sec. 9.2].
4 The \({{\mathbf {R}}}^n\) Cauchy Integral
A main reason for using the Dirac framework is that it provides us with a Cauchy-type reproducing formula, which allows for a generalization of complex function theory to \(n\ge 2\) and \(k\ne 0\). See [37, Chapters 8,9] for further details. More precisely, if F satisfies (41) in a domain \(\Omega \) with boundary \(\partial \Omega \), then a Cauchy-type reproducing formula
holds. We write dy and \(\nu \in \wedge ^1{{\mathbf {R}}}^n\) for the standard measure and outward pointing unit normal on \(\partial \Omega \) respectively, and the integrand uses two Clifford products. The first factor
is a fundamental solution to the elliptic operator \({{\mathbf {D}}}+ik\), as \({{\mathbf {D}}}\Psi _k(x)+ik\Psi _k(x)= \delta (x)\). Here \(\Phi _k(x)\in \wedge ^0{{\mathbf {R}}}^n={{\mathbf {C}}}\) is the Helmholtz fundamental solution, and we use the normalization from [11] so that \((\Delta +k^2)\Phi _k= -2\delta (x)\). Hence the factor \(-1/2\) in (45). In dimension \(n=3\) we have
and in dimension \(n=2\) we have in terms of the Hankel function \(H^{(1)}_0(z)\) that
The classical theory of complex Hardy spaces generalizes from complex function theory to our Dirac setting. Our basic operator, acting on a suitable space \({{\mathcal {H}}}\) of functions \(h: \partial \Omega \rightarrow \wedge {{\mathbf {R}}}^n\) on \(\partial \Omega \), 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_k^2=I\), and
are two complementary projection operators, that is \((E_k^\pm )^2= E_k^\pm \). The operator \(E_k^\pm \) projects onto its range, the subspace of \({{\mathcal {H}}}\) which we denote \(E_k^\pm {{\mathcal {H}}}\) and consists of all traces \(F|_{\partial \Omega }\) of fields satisfying \({{\mathbf {D}}}F= ikF\) in \(\Omega ^\pm \). For the exterior domain \(\Omega ^-\) these fields also satisfy a Dirac radiation condition; See (58) 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\in \partial \Omega \). The single layer \(S^a_k\) is a weakly singular integral operator, and so is also \(K^v_k\) if \(\partial \Omega \) is smooth and \(v=\nu \) 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 \({{\mathcal {H}}}\), also for general Lipschitz interfaces \(\partial \Omega \).
Consider first dimension \(n=2\), with a curve \(\partial \Omega \) and \(\Phi _k\) given by (47). Here we use a positively oriented frame \(\{\nu ,\tau \}\) at \(x\in \partial \Omega \), with \(\nu =\nu (x)\) being the normal vector into \(\Omega ^-\) and \(\tau =\tau (x)\) the tangential vector counter clock-wise from \(\nu \). The corresponding frame at \(y\in \partial \Omega \) we write as \(\nu '=\nu (y)\) and \(\tau '=\tau (y)\). In the plane, the Clifford algebra \(\wedge {{\mathbf {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\in \partial \Omega \), we write
using instead the vector frame \(\{\nu ,\tau \}\). Here \(\nu \tau =e_{12}\in \wedge ^2{{\mathbf {R}}}^2\) does not depend on x, even though \(\nu \) and \(\tau \) do so. By writing out \(\nabla \Phi _k\) and the Clifford products in (48), we obtain
in the multivector frame \(\{1,\nu \tau , \nu ,\tau \}\).
Next consider dimension \(n=3\), with a surface \(\partial \Omega \) and \(\Phi _k\) given by (46). Here we use a positively oriented ON-frame \(\{\nu ,\tau ,\theta \}\) at \(x\in \partial \Omega \), with \(\nu =\nu (x)\) being the normal vector into \(\Omega ^-\) and \(\tau =\tau (x)\) and \(\theta =\theta (x)\) being tangential vector fields. The frame vectors at \(y\in \partial \Omega \), we write as \(\nu '=\nu (y)\), \(\tau '=\tau (y)\) and \(\theta '=\theta (y)\). A multivector field at \(x\in \partial \Omega \) we write as
Here \(\nu \tau \theta =e_{123}\in \wedge ^3{{\mathbf {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,\tau \theta ,\theta \nu ,\nu \tau ,\nu \tau \theta ,\nu ,\tau ,\theta \}\), we have
Finally we note that we similarly can write the Cauchy integral (44) for the fields in \(\Omega ^\pm \) , in matrix form. The only difference is the normalization factor 2 and the fact that we choose the frame \(\{e_1,\ldots , e_n\}\) at the field point \(x\in \Omega ^\pm \). We also allow for a general function \(h:\partial \Omega \rightarrow \wedge {{\mathbf {R}}}^n\) and not only a trace f of a solution to the Dirac equation in (44). By the associativity of the Clifford product, this still yields a field F solving the Dirac equation, but \(h\ne F|_{\partial \Omega }\) in general. In this case, when \(x\notin \partial \Omega \), we denote the layer potentials by \({\tilde{K}}^v_k\) and \({\tilde{S}}^a_k\), with v and a now only depending on \(y\in \partial \Omega \). The field evaluation formulas (13), (22) and (23) also use vector versions
of \({\tilde{K}}^v_k\) and \({\tilde{S}}^a_k\), where now \(A(y):{{\mathbf {R}}}^n\rightarrow {{\mathbf {R}}}^n\) is a matrix function and v(y) is a vector field.
5 The Energy Trace Space
We define in this section the appropriate norms and function spaces for the Dirac equation (41). Let \(\Omega ^+\subset {{\mathbf {R}}}^n\) be a bounded Lipschitz domain. The general Dirac transmission problem that we consider reads
This is our master transmission problem, into which we embed the Helmholtz and Maxwell transmission problems in Sects. 7 and 8, where the multiplier M will be specified. The radiation condition stated here for the Dirac equation reduces to the classical Sommerfeld and Silver–Müller conditions for Helmholtz and Maxwell solutions \(F_{\text {helm}}\) and \(F_{\mathrm{em}}\,\) respectively. See [37, Sec. 9.3].
For \(F=F^+\) in \(\Omega ^+\), we use the norm
It is important to note that for both Helmholtz and Maxwell fields \(F_{\text {helm}}\) and \(F_{\mathrm{em}}\,\), the second term \(\nabla \mathbin {\scriptstyle {\wedge }}F\) is bounded by the first term F, and the norm reduces to the \(L_2\) norm of F. In \(\Omega ^-\) we use the corresponding norm, but integrate only over the bounded set \(\Omega ^-_R\). Different choices of R yield equivalent norms for solutions to the Dirac equation which satisfy the radiation condition. See [37, Sec. 9.3].
It was shown in [4] that there exists a unique Hilbert space \({{\mathcal {H}}}\) (although it does not come with a canonical norm) of multivector fields on the Lipschitz interface \(\partial \Omega \), which is the trace space corresponding to the above norms on \(F^\pm \) in \(\Omega ^\pm \), with the following properties. It splits into closed subspaces in two ways as
where \({\scriptstyle \top }\) means tangential and \({\scriptstyle \perp }\) means normal. There is a one-to-one correspondence \(F^+\leftrightarrow f^+\) between Dirac solutions \(F^+\) in \(\Omega ^+\) and their traces \(f^+\in E_k^+{{\mathcal {H}}}\), with inverse given by the Cauchy integral (44). Similarly, there is a one-to-one correspondence \(F^-\leftrightarrow f^-\) between Dirac solutions \(F^-\) in \(\Omega ^-\) satisfying the radiation condition, and their traces \(f^-\in E_k^-{{\mathcal {H}}}\), with inverse given by the Cauchy integral.
The subspace \({{\mathcal {H}}}_{\scriptstyle \top }\) denotes the subspace of tangential multivector fields, and there is a bounded and surjective tangential trace map
of multivector fields F in a neighbourhood U of \(\partial \Omega \) with norm \((\int _U (|F(x)|^2+ |\nabla \mathbin {\scriptstyle {\wedge }}F(x)|^2) dx )^{1/2}<\infty \). The subspace \({{\mathcal {H}}}_{\scriptstyle \perp }\) denotes the subspace of normal multivector fields, and there is a bounded and surjective normal trace map
of multivector fields F in a neighbourhood U of \(\partial \Omega \) with norm \((\int _U (|F(x)|^2+ |\nabla \mathbin {\lrcorner }F(x)|^2) dx )^{1/2}<\infty \).
For smooth interfaces \(\partial \Omega \), the trace space \({{\mathcal {H}}}\) consists of all \(f\in H^{-1/2}(\partial \Omega ;\wedge {{\mathbf {R}}}^n)\) such that \(d_{\scriptstyle \top }f_{\scriptstyle \top }\in H^{-1/2}(\partial \Omega ;\wedge {{\mathbf {R}}}^n)\) and \(\delta _{\scriptstyle \top }(\nu \mathbin {\lrcorner }f)\in H^{-1/2}(\partial \Omega ;\wedge {{\mathbf {R}}}^n)\), with \(d_{\scriptstyle \top }\) and \(\delta _{\scriptstyle \top }\) denoting the tangential exterior and interior derivatives along \(\partial \Omega \). Here \(f_{\scriptstyle \top }\) denotes the tangential part of f and \(\nu \mathbin {\lrcorner }f\) is the normal part of f with the factor \(\nu \) removed.
To characterize \({{\mathcal {H}}}\) in terms of fractional Sobolev-type spaces on non-smooth Lipschitz interfaces \(\partial \Omega \) is a non-trivial problem. Such a characterization was achieved for polyhedra in [7, 8], for general Lipschitz interfaces in [9] in \({{\mathbf {R}}}^3\), and extended to general multivector traces in [42]. For the energy trace space space \({{\mathcal {H}}}\) on Lipschitz interfaces, this shows the following. For 2D domains, using the frame (52), our boundary function space is
It is for 3D domains that the non-trivial result from [9, Thm. 4.1] is needed, which shows that in the frame (54), we have
Here \(*\) is the 3D Hodge star and the space of tangential vector fields is defined as
where \({{\mathrm{curl}}}_{\scriptstyle \top }\) denotes tangential surface curl on \(\partial \Omega \) and \((\nu \times H^{1/2}(\partial \Omega )^3)^*\) denotes the dual space of \(\nu \times H^{1/2}(\partial \Omega )^3\). One should note that the test functions \(\nu \times H^{1/2}(\partial \Omega )^3\) will not be \(H^{1/2}\) smooth on general Lipschitz regular \(\partial \Omega \), since \(\nu \) is in general only measurable. To summarize, the conditions on \(h_3,h_4,h_7,h_8\) are that \(h_3\tau +h_4\theta \in H^{-1/2}({{\mathrm{curl}}},\partial \Omega )\) and \(h_7\tau +h_8\theta \in H^{-1/2}({{\mathrm{curl}}},\partial \Omega )\).
We remark that there is no canonical norm on \({{\mathcal {H}}}_2\) or \({{\mathcal {H}}}_3\), but this does not present a problem and it should be clear from the context in the estimates to come, which choice among equivalent norms is used when such needs to be fixed.
Throughout this paper \(X\lesssim Y\) means that there exists \(C<\infty \) independent of relevant variables so that \(X\le CY\), and \(X\approx Y\) means that \(X\lesssim Y\) and \(Y\lesssim X\).
6 Helmholtz Existence and Uniqueness
This section surveys basic solvability results for the Helmholtz transmission problem (4), which are valid in any dimension \(n\ge 2\). The proofs follow [4], with a translation from the Dirac to the Helmholtz framework. Consider first uniqueness of solutions \(U^\pm \).
Proposition 6.1
Let \(\Omega ^+\subset {{\mathbf {R}}}^n\) be a bounded Lipschitz domain with connected exterior \(\Omega ^-\), and consider a solution \(U^\pm \) to (4). Assume that the incoming wave vanishes so that \(u_0=0\). If \({\hat{\epsilon }}\) from (4) satisfies
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1, then \(U^+=U^-=0\) identically.
Proof
From the jump relations we have
Apply Green’s first identity for \(\Omega ^+\) and \(\Omega ^-_R\), use the Helmholtz equation for \(U^\pm \) and the radiation condition for \(U^-\), and multiply the equation by \({{\hat{\epsilon }}}k_-/i\) to obtain
as \(R\rightarrow \infty \). Denote the three integrals appearing in (68) by \(I^+\), \(I^-_R\) and \(I_R\) respectively, including \(i^{-1}\) in the first two, and set \(\phi _\pm :=|\arg (k_\pm /i)|\) as in Definition 2.1, let \(z:= {{\hat{\epsilon }}} k_-/k_+\), and define the sector
We note that
and similarly for \(I^-_R\).
Sectors and lines appearing, depending on \(\phi _\pm \)
We verify that the condition (66) implies \(U^\pm =0\), by examining (68) in the nine cases \(\phi _\pm =0\), \(0<\phi _\pm <\pi /2\) and \(\phi _\pm =\pi /2\) as follows.
-
(i)
Assume \(0<\phi _-<\pi /2\) and \(0<\phi _+<\pi /2\). Then \(U^-\) decays exponentially as \(R\rightarrow \infty \). Setting \(R=\infty \) in (68), we have the equation \(z I^++|{{\hat{\epsilon }}}|^2I^-_\infty =0\), where \(I^+\in S_{\phi _+}\) and \(I^-_\infty \in S_{\phi _-}\). If \(|\arg (z)|+\phi _-+\phi _+<\pi \), then this is possible only if \(I^+=I^-_\infty =0\). See Fig. 2i. Indeed, the sectors \(zS_{\phi _+}\) and \(-S_{\phi _-}\) intersect only at 0. We conclude in particular that \({{\,\mathrm{Re}\,}}I^+=0={{\,\mathrm{Re}\,}}I^-_\infty \), and therefore \(U^+=U^-=0\) according to (70).
If \(|\arg (z)|+\phi _-+\phi _+=\pi \), then the rotated sectors \(zS_{\phi _+}\) and \(-S_{\phi _-}\) touch and we observe that \(I^+\) and \(I^-_\infty \) lie on the boundary of respective sector \(S_{\phi _\pm }\). From (70), we conclude that \(\int _{\Omega ^+}(|\nabla U^+|^2-|k_+U^+|^2)dx=\pm \int _{\Omega ^+}(|\nabla U^+|^2+|k_+U^+|^2)dx\). This forces either \(U^+=0\) or \(\nabla U^+=0\). We conclude from the Helmholtz equation that \(U^+=0\). A similar argument shows that \(U^-=0\).
-
(ii)
Assume \(\phi _\pm <\pi /2\) and either \(\phi _+=0\) or \(\phi _-=0\). When \(|\arg (z)|+\phi _-+\phi _+<\pi \), the argument in (i) applies. When \(|\arg (z)|+\phi _-+\phi _+=\pi \) and \(\phi _+\ne 0\), the argument in (i) shows that \(U^+=0\). The jump condition then shows that \(U^-=0\). If instead \(|\arg (z)|+\phi _-+\phi _+=\pi \) and \(\phi _-\ne 0\), a similar argument shows that \(U^-=0\), from which \(U^+=0\) follows from the jump conditions.
-
(iii)
Assume \(\phi _-<\pi /2\) and \(\phi _+=\pi /2\). We now have the equation \(z I^++|{{\hat{\epsilon }}}|^2I^-_\infty =0\), with \({{\,\mathrm{Re}\,}}I^+=0\) and \(I^-_\infty \in S_{\phi _-}\). If \(\min (|\arg (z)|,|\arg (-z)|)< \pi /2-\phi _-\), then the line \(zS_{\phi _+}\) intersects the sector \(-S_{\phi _-}\) only at the origin, which shows that \(I^+=I^-_\infty =0\). See Fig. 2iii. From the exterior (\(\Omega ^-\)) analogue of (70), we conclude \(U^-=0\), and jump conditions imply \(U^+=0\).
If \(\min (|\arg (z)|,|\arg (-z)|)= \pi /2-\phi _-\) and if \(\phi _->0\), then it follows that \(I^-_\infty \) lies on the boundary of \(S_{\phi _-}\). As in (i), this shows that \(U^-=0\), and by jump conditions that \(U^+=0\).
-
(iv)
Assume \(\phi _-=\pi /2\). Then \({{\,\mathrm{Re}\,}}I^-_R=0\), and (68) reduces to \({{\,\mathrm{Re}\,}}(z I^+)+ \lim _{R\rightarrow \infty }I_R=0\). If \(|\arg (z)|+\phi _+\le \pi /2\), then \({{\,\mathrm{Re}\,}}(z I^+)\ge 0\) and \(\lim _{R\rightarrow \infty } I_R=0\) follows. See Fig. 2iv. By Rellich’s lemma this implies that \(U^-=0\). The jump relations then shows that \(U^+=0\).
If also \(\phi _+=\pi /2\) then \({{\,\mathrm{Re}\,}}I^+=0\), and we conclude that \(\lim _{R\rightarrow \infty } I_R=0\) also when \(z<0\), and can in the same way conclude that \(U^-=0=U^+\).
\(\square \)
Next consider the existence of solutions \(U^\pm \). The following result is essentially from [4], where more details and background can be found. For a short survey of the Fredholm theory that we apply, we refer to [37, Sec. 6.4].
Proposition 6.2
Let \(\Omega ^+\subset {{\mathbf {R}}}^n\) be a bounded Lipschitz domain with connected exterior \(\Omega ^-\). Then there exists \(1\le C(\partial \Omega )<\infty \) such that if
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1, then there exists a unique solution \(U^+\in H^1(\Omega ^+)\), \(U^-\in H^1(\Omega ^-_R)\) to the Helmholtz transmission problem (4) in \({{\mathbf {R}}}^n\), \(n\ge 2\), depending continuously on the datum \(u^0\in H^{1/2}(\partial \Omega )\).
Proof
-
(i)
We first use Fredholm theory to reduce the problem to an estimate of
$$\begin{aligned} \int _{\Omega ^+}|\nabla U^+|^2 dx+\int _{\Omega ^-_R}|\nabla U^-|^2 dx \end{aligned}$$(72)by the \({{\mathcal {H}}}\) norm of \(u^0\) and compact terms. To this end, it is convenient to consider the multivector fields \(F^\pm =\nabla U^\pm +ik_\pm U^\pm \) and \(F^0= \nabla U^0+ik_-U^0\) as in (42). For fixed \(k_\pm \), we define function spaces \({\widetilde{{{\mathcal {H}}}}}_{k_\pm }^{1,\pm }\) of such \(F^\pm \) in \(\Omega ^\pm \), with potentials \(U^\pm \) solving \(\Delta U^\pm + k_\pm ^2 U^\pm =0\), and \(U^-\) satisfying the radiation condition in \(\Omega ^-\). The norm of \(F^+\in {\widetilde{{{\mathcal {H}}}}}^{1,+}_{k_+}\) is \(\Vert F^+\Vert _{L_2(\Omega ^+)}\) and the norm of \(F^-\in {\widetilde{{{\mathcal {H}}}}}^{1,-}_{k_-}\) is \(\Vert F^-\Vert _{L_2(\Omega ^-_R)}\), with a fixed large R. For the data/right-hand sides, we define the subspace \({{\mathcal {H}}}^1_{k_-}\subset {{\mathcal {H}}}\) consisting of \(f= f_0+f_{1{\scriptstyle \top }}+ f_{1{\scriptstyle \perp }}\), with \(f_0\) scalar and \(f_{1{\scriptstyle \top }}, f_{1{\scriptstyle \perp }}\) tangential and normal vector fields respectively, satisfying the constraint \(\nabla _{\scriptstyle \top }f_0= ik_-f_{1{\scriptstyle \top }}\). The traces of incoming waves \(F^0\) all belong to \({{\mathcal {H}}}^1_{k_-}\).
The transmission problem (4) amounts to inverting a bounded linear map
$$\begin{aligned} T_{{\hat{\epsilon }}, k_+, k_-}: {\widetilde{{{\mathcal {H}}}}}_{k_+}^{1,+} \oplus {\widetilde{{{\mathcal {H}}}}}_{k_-}^{1,-}\rightarrow {{\mathcal {H}}}^1_{k_-}. \end{aligned}$$(73)Assume first that \(k_+=k_-\) and that \({\hat{\epsilon }}=1\). In this case the jump condition in (4) reduces to \(f^+-f^-=f^0\), and it follows that \(T_{{\hat{\epsilon }}, k_+, k_-}\) is invertible since the Cauchy integrals (44), with \(x\in \Omega ^\pm \), provide an explicit inverse.
Next consider general parameters \({\hat{\epsilon }}, k_+, k_-\) satisfying (71). It suffices to show that \(T_{{\hat{\epsilon }}, k_+, k_-}\) is a Fredholm operator with index zero. Indeed, Proposition 6.1 shows injectivity, from which surjectivity then follows. To prove index zero, we note that \({{\mathbf {C}}}\setminus [-C(\partial \Omega ), -1/C(\partial \Omega )]\) is an open connected set, and aim to apply Fredholm perturbation theory, perturbing \({\hat{\epsilon }}\) to 1 and \(k_+\) to \(k_-\). We prove in (ii) below that the operator \(T_{{\hat{\epsilon }}, k_+, k_-}\) is semi-Fredholm whenever \({\hat{\epsilon \notin }} [-C(\partial \Omega ), -1/C(\partial \Omega )]\), so it remains to verify that the operator and function spaces depend continuously on the parameters.
From Hodge decompositions of the space \({{\mathcal {H}}}\), see [4], it follows that for \(k_-\ne 0\) there are projections \({{\mathcal {H}}}\rightarrow {{\mathcal {H}}}^1_{k_-}\) onto these subspaces, depending continuously on \(k_-\). This uses the Hodge projections not for the exterior and interior derivatives, but for zero order perturbations of these defined by the wave number \(k_-\) (c.f. (74) below). The crucial observation is that for all \(k_-\ne 0\), the cohomology for this perturbed Hodge decomposition vanishes. This implies that the perturbed Hodge projections depend continuously on \(k_-\).
For the domain spaces \({\widetilde{{{\mathcal {H}}}}}_{k_\pm }^{1,\pm }\), we note that the trace map and the Cauchy integral give an isomorphism between \(F^\pm \in {\widetilde{{{\mathcal {H}}}}}_{k_\pm }^{1,\pm }\) and \(f^\pm \in {{\mathcal {H}}}_{k_\pm }^{1,\pm }\subset E_k^\pm {{\mathcal {H}}}\subset {{\mathcal {H}}}\). Like for \({{\mathcal {H}}}^1_{k_-}\), it follows from Hodge decompositions of \({{\mathcal {H}}}\) that for \(k_\pm \ne 0\) there are projections \({{\mathcal {H}}}\rightarrow {{\mathcal {H}}}_{k_\pm }^{1,\pm }\) onto these trace spaces for Helmholtz fields, depending continuously on \(k_\pm \). Since \(T_{{\hat{\epsilon }}, k_+, k_-}\) clearly depend continuously on \({\hat{\epsilon }}\), perturbation theory applies to show index zero, provided we show the estimate (72).
-
(ii)
To establish the estimate of (72), we construct certain auxiliary potentials \(V^\pm \) to the gradient vector fields \(\nabla U^\pm \). In \(\Omega ^-_R\), we simply use the given scalar potential \(V^-=U^-\). In \(\Omega ^+\), we find a bivector field \(V^+:\Omega ^+\rightarrow \wedge ^2{{\mathbf {R}}}^n\) and a vector field \(\tilde{V}^+:\Omega ^+\rightarrow \wedge ^1{{\mathbf {R}}}^n\) such that
$$\begin{aligned} {\left\{ \begin{array}{ll} \nabla U^+= \nabla \mathbin {\lrcorner }V^++ ik_+ {\tilde{V}}^+,\\ ik_+ U^+= \nabla \mathbin {\lrcorner }{\tilde{V}}^+. \end{array}\right. } \end{aligned}$$(74)(Modulo the term \({\tilde{V}}^+\), 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(\Omega ^+)\rightarrow L_2(\Omega ^+)^2: U^+\mapsto (V^+,{\tilde{V}}^+)\) follows from Hodge decompositions as in [6], after translation from the spacetime framework using [37, Sec. 9.1]. To complete the construction of \(V^\pm \), we extend the potentials \(V^\pm \) to compactly supported functions on \({{\mathbf {R}}}^n\), with \(\nabla V^-\) and \(\nabla \mathbin {\lrcorner }V^+\) belonging to \(L_2({{\mathbf {R}}}^n)\).
Pairing the jump relations with \(v^\pm \), we have
$$\begin{aligned} {\left\{ \begin{array}{ll} \int _{\partial \Omega }\langle \nu \mathbin {\scriptstyle {\wedge }}\nabla u^+,v^+ \rangle d\sigma = \int _{\partial \Omega }(\langle \nu \mathbin {\scriptstyle {\wedge }}\nabla u^-,v^+ \rangle +\langle \nu \mathbin {\scriptstyle {\wedge }}\nabla u^0,v^+ \rangle ) d\sigma ,\\ \int _{\partial \Omega }\langle v^-,\nu \mathbin {\lrcorner }\nabla u^+ \rangle d\sigma =\overline{{\hat{\epsilon }}} \int _{\partial \Omega }(\langle v^-,\nu \mathbin {\lrcorner }\nabla u^- \rangle + \langle v^-,\nu \mathbin {\lrcorner }\nabla u^0 \rangle )d\sigma . \end{array}\right. } \end{aligned}$$(75)Using the general Stokes theorem, see [37, Sec. 7.3], we have modulo compact terms
$$\begin{aligned} \int _{\partial \Omega }\langle \nu \mathbin {\scriptstyle {\wedge }}\nabla u^+,v^+ \rangle d\sigma\approx & {} \int _{\Omega ^+}|\nabla U^+|^2 dx, \end{aligned}$$(76)$$\begin{aligned} \int _{\partial \Omega }\langle v^-,\nu \mathbin {\lrcorner }\nabla u^- \rangle d\sigma\approx & {} -\int _{\Omega ^-_R}|\nabla U^-|^2 dx, \end{aligned}$$(77)and
$$\begin{aligned}&\int _{\partial \Omega }\langle v^-,\nu \mathbin {\lrcorner }\nabla u^+ \rangle d \sigma - \int _{\partial \Omega }\langle \nu \mathbin {\scriptstyle {\wedge }}\nabla u^-,v^+ \rangle d\sigma \nonumber \\&\quad \approx \int _{\Omega ^+}\langle \nabla V^-,\nabla U^+ \rangle dx+ \int _{\Omega ^-_R}\langle \nabla U^-,\nabla \mathbin {\lrcorner }V^+ \rangle dx\nonumber \\&\quad \approx \int _{{{\mathbf {R}}}^n}\langle \nabla V^-,\nabla \mathbin {\lrcorner }V^+ \rangle dx=0. \end{aligned}$$(78)Therefore, adding the equations (75) yields an estimate of (72) whenever we do not have \({{\hat{\epsilon }}}<0\). When \({{\hat{\epsilon }}}\) is negative and close to \(\infty \), we instead subtract the equations to conclude. When \({{\hat{\epsilon }}}\) is negative and close to 0, we can also obtain an estimate of (72) by instead starting with a bivector potential \(V^+\) in \(\Omega ^+\), and a scalar potential \(V^-\) in \(\Omega ^-\).
\(\square \)
7 The 2D Dirac Integral
In this section, we derive the Dirac integral equation (10) in 2D by combining two Helmholtz problems, and using a duality ansatz. We start by recalling how 2D Maxwell transmission problems reduce to Helmholtz transmission problems. We end by optimizing the Dirac parameters \(r,\beta ,\alpha ',\beta '\), which is a main step in the construction of the Dirac BIE (10).
Example 7.1
We consider applications of the Dirac equation to the scattering of EM fields as in (43), 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_1e_1+V_2e_2\), the Dirac equation \({{\mathbf {D}}}F=ikF\) amounts to \(V=\nabla U\) and \({{\mathrm{div}}}V= -k^2 U\), that is the Helmholtz equation \(\Delta U +k^2 U=0\) for the scalar function U. We have
With this setup for both domains \(\Omega ^\pm \), jump relations for the electromagnetic field specify the jump matrix
for \(F^\pm \) in the frame \(\{1,\nu \tau , \nu ,\tau \}\). 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.
Proposition 7.2
Consider the Dirac transmission problem (58) for a bounded Lipschitz domain \(\Omega ^+\subset {{\mathbf {R}}}^2\) with connected exterior \(\Omega ^-\). Let
with parameters \(\beta ,\alpha \in {{\mathbf {C}}}\setminus \{0\}\). Assume that
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1. Then the operator \(B_2: E_{k_+}^+{{\mathcal {H}}}_2\oplus E_{k_-}^-{{\mathcal {H}}}_2 \rightarrow {{\mathcal {H}}}_2\) given by
is invertible, where we as in Sect. 5 denote the spaces of traces of solutions on \(\partial \Omega \) by \(E_{k_\pm }^\pm {{\mathcal {H}}}_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^\pm _2\) solve the Helmholtz equation as a consequence of \(F^\pm \) solving the Dirac equation, with wavenumbers \(k_\pm \) respectively. Moreover, the vector part of \({{\mathbf {D}}}F^\pm = ik_\pm F^\pm \) shows that
$$\begin{aligned} \nabla F^\pm _2= (\nabla F^\pm _0-ik_\pm F^\pm _1)e_{12}. \end{aligned}$$(90)From the assumed jump relations for \(f^\pm \) it therefore follows that \(f_2^+= {\hat{k}}f^-_2/\beta \) and \(\partial _\nu f_2^+= \beta \partial _\nu ({\hat{k}}f^-_2/\beta )\), and hence Proposition 6.1 with \({{\hat{\epsilon }}}=\beta \) shows that \(f_2^\pm =0\).
-
(ii)
Next consider the scalar functions \(F^\pm _0\), which also solve the Helmholtz equation. Since \(F^\pm _2=0\) by (i), we have \(\nabla F^\pm _0=ik_\pm F^\pm _1\) and obtain jump relations \(f_0^+= \hat{k}f^-_0\) and \(\partial _\nu f_0^+= \alpha \partial _\nu ({\hat{k}}f^-_0)\). Again Proposition 6.1 applies, now with \({{\hat{\epsilon }}}=\alpha \), and shows that \(f_0^\pm =0\). From (90) we conclude \(f^\pm _1=0\) and in total \(f^\pm =0\).
To show existence, it suffices by perturbation theory for Fredholm operators to prove an estimate
$$\begin{aligned} \Vert f^+\Vert _{{{\mathcal {H}}}_2}+ \Vert f^-\Vert _{{{\mathcal {H}}}_2}\lesssim \Vert f^0\Vert _{{{\mathcal {H}}}_2}. \end{aligned}$$(91)This follows as in steps (i) and (ii) but instead appealing to Proposition 6.2. In this case, we obtain in step (i) that \(\Vert f^+_2\Vert _{{{\mathcal {H}}}_2}+\Vert f^-_2\Vert _{{{\mathcal {H}}}_2}\lesssim \Vert f^0\Vert _{{{\mathcal {H}}}_2}\), which in step (ii) is used to estimate \(f^\pm _0\) and \(f^\pm _1\).
\(\square \)
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.
Proposition 7.3
Assume the hypothesis of Proposition 7.2, and further assume that
with parameters satisfying
Then the operator
is invertible on \({{\mathcal {H}}}_2\), for any \(r\in {{\mathbf {C}}}\setminus \{0\}\).
Proof
We factorize (94) via \(E_{k_+}^+{{\mathcal {H}}}_2\oplus E_{k_-}^-{{\mathcal {H}}}_2\) as
By Proposition 7.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 \({{\mathcal {H}}}_2\), where \({\widehat{w}}= \sum _j (-1)^j w_j\) denotes the involution of a multivector \(w=\sum _j w_j\), \(w_j\in \wedge ^j {{\mathbf {R}}}^n\). It is readily verified that
where \({\widetilde{M}}={\mathrm{diag}}\,\begin{bmatrix} 1/{\hat{k}}&\quad (1/{\hat{k}})/\beta '&\quad \alpha '&\quad 1 \end{bmatrix}\). This shows that in the natural way \(E^\mp _{k_\pm }{{\mathcal {H}}}_2\) is the dual space of \(E^\pm _{k_\pm }{{\mathcal {H}}}_2\). Hilbert space duality theory shows that in (95), invertibility of the right factor is equivalent to invertibility of the left factor, with M replaced by \({\widetilde{M}}\), and \(k_-\) and \(k_+\) swapped. \(\square \)
Consider the Dirac integral equation
involving the operator from (94), with right-hand side specified by the jump condition in (58) and auxiliary density \({\tilde{h}}\) which we precondition as \({\tilde{h}}= P'h\) in the proof of Theorem 2.2. We optimize (98) by choosing the parameters \(r,\beta ,\alpha ',\beta '\). Recall that for EM fields \(\alpha = {\hat{\epsilon }}\), where our main interest is \({\mathrm{Im}}\,{\hat{\epsilon \ge }} 0\) and \({\hat{k}}= \sqrt{{\hat{\epsilon }}}\), so that \({{\,\mathrm{Re}\,}}{\hat{k}}\ge 0\). We therefore consider \(\alpha \) as having a prescribed value. Clearly
where
Let \(K^v\) and \(S^a\) be the static double and single layer type operators, that is \(K^v= K^v_0\), and \(S^a\) is \(S^a_k\) without the factor ik at \(k=0\). Modulo operators of the form \(K^v_{k_\pm }-K^v\), \(S^a_{k_+}-ik_-{\hat{k}}S^a\) and \(S^a_{k_-}-ik_-S^a\), the integral operator \(T=rE_{k_+} M'-M E_{k_-}\) from (99) is the entry-wise product of
and \(E_k\), with diagonal and off-diagonal \(2\times 2\) blocks replaced by \(K^v\) and \(ik_- S^a\) respectively. Indeed, under these approximations
-
Our first choice is to set \(r={\hat{k}}\). This gives cancellation in the (1,2) and (4,3) elements of the operator T, which for a smooth domain \(\Omega \) yields \(T^2=0\) on \({{\mathcal {H}}}_2\) modulo compact operators. In particular the essential spectrum of T is \(\{0\}\).
-
Our choices for \(\beta , \beta ', \alpha '\) are to set \(\beta = {\hat{k}}/|{\hat{k}}|\) and \(\beta '=\alpha '= \overline{{\hat{k}}}/|{\hat{k}}|\). These choices make \(\beta /{\hat{k}}>0\), \(\beta ' {\hat{k}}>0\) and \(\alpha ' {\hat{k}}>0\). Therefore, if \(\alpha \in {{\mathbf {C}}}\setminus [-C(\partial \Omega ), -1/C(\partial \Omega )]\) and \(\alpha /{\hat{k}}\in {\mathrm{WP}}\,(k_-,k_+)\), then Propositions 7.2 and 7.3 guarantee invertibility of \(E_{k_+}^+M'+M E_{k_-}^-\).
Furthermore, the choice of \(\beta , \beta ',\alpha '\) yields diagonal (1,1), (2,2) and (4,4) elements in (94) which are compact perturbations of invertible operators on any Lipschitz domain, whenever \({\hat{k}}\in {{\mathbf {C}}}\setminus (-\infty ,0]\). Indeed, \(K^\nu _{k_\pm }-K^\nu _0\) is compact and the essential spectrum of \(K^\nu _0\) is contained in \((-1,1)\). Moreover, when \({{\,\mathrm{Re}\,}}{\hat{k}}\ge 0\) then the normalization \(|\beta |= |\beta '|=|\alpha '|=1\) gives spectral points \(\lambda \) for \(K^\nu _0\) on the imaginary axis with \(|\lambda |\ge 1\), since \(\lambda = (1+z)/(1-z)\) maps \(\{|z|=1,{{\,\mathrm{Re}\,}}z\ge 0\}\) onto \(\{{{\,\mathrm{Re}\,}}\lambda =0, |\lambda |\ge 1\}\).
To summarize, for solving the Helmholtz/TM Maxwell transmission problem as described above, we have obtained the 2D Dirac integral equation
with
8 The 3D Dirac Integral
In Sects. 6 and 7, we derived an integral equation for solving Dirac transmission problems in \({{\mathbf {R}}}^2\), which applies to Helmholtz/TM Maxwell scattering. We here derive the completely analogous integral equation in \({{\mathbf {R}}}^3\), with applications to scattering for the full Maxwell equations and not only the Helmholtz equation. This Dirac integral equation in 3D combines one Maxwell problem and two Helmholtz problems, and uses a duality ansatz. We end this section by optimizing the Dirac parameters \(r,\beta ,\gamma ,\alpha ',\beta ',\gamma '\), which is a main step in the construction of the Dirac BIE (20).
Example 8.1
Maxwell’s equations correspond to an electromagnetic field F with \(F_0=0=F_3\) as in (43). The energy norm that we consider is simply the \(L_2\) norm of \(F^\pm \) in \(\Omega ^+\) and \(\Omega ^-_R\), respectively, and the corresponding function space \({{\mathcal {H}}}_3\) on \(\partial \Omega \) is \((H^{-1/2}(\partial \Omega ))^6\), with tangential curls of E and B belonging to \(H^{-1/2}(\partial \Omega )\).
In the 3D Dirac transmission problem (58), Maxwell scattering for the field \(F=F_{\mathrm{em}}\,\) from (43) specifies the jump matrix
in the frame (54), by the jump relations for the electric and magnetic field. The parameters a and b can be chosen freely since \(F_0=0= F_3\).
Consider the Maxwell transmission problem (14), which we write in multivector notation, with \(F_1^\pm = E^\pm \) and \(F_2^\pm \) being the Hodge dual of \(B^\pm \), as follows.
We require the following Maxwell versions of the results in Sect. 6.
Proposition 8.2
Consider the Maxwell transmission problem (107) for a bounded Lipschitz domain \(\Omega ^+\subset {{\mathbf {R}}}^3\) with connected exterior \(\Omega ^-\). Assume that the incoming wave vanishes so that \(f^0=0\). If \({\hat{\epsilon }}\) from (107) satisfies
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1, then \(F^+=F^-=0\) identically.
Proof
Similar to the proof of Proposition 6.1, we use the jump relations to obtain
We then apply a Stokes theorem for \(\Omega ^+\) and \(\Omega ^-_R\), to obtain
as \(R\rightarrow \infty \). We here used that \(\nabla \mathbin {\lrcorner }F^\pm _2= ik_\pm F^\pm _1\), \(\nabla \mathbin {\scriptstyle {\wedge }}F^\pm _1= ik_\pm F^\pm _2\), and by the radiation condition that \(\langle F_1^-,\nu \mathbin {\lrcorner }F_2^- \rangle \approx |F^-_1|^2\approx |F^-_2|^2\approx \tfrac{1}{2}|F^-|^2\) on \(|x|=R\). Using (110), the result follows similarly to the proof of Proposition 6.1. \(\square \)
Proposition 8.3
Let \(\Omega ^+\subset {{\mathbf {R}}}^3\) be a bounded Lipschitz domain with connected exterior \(\Omega ^-\). Then there exists \(1\le C(\partial \Omega )<\infty \) such that if \({\hat{\epsilon }}\) from (107) satisfies
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1, then there exists a unique solution \(F^+\in L_2(\Omega ^+)\), \(F^-\in L_2(\Omega ^-_R)\) to the Maxwell transmission problem (107), depending contiuously on the trace \(f^0\in {{\mathcal {H}}}_3\) of the incoming electromagnetic wave \(F^0\).
Proof
-
(i)
The proof is similar to that of Proposition 6.2, replacing the scalar function \(ik_\pm U^\pm \) by the vector field \(F^\pm _1\), and the vector field \(\nabla U^\pm \) by the bivector field \(F^\pm _2\). We first use Fredholm theory to reduce the problem to an estimate of
$$\begin{aligned} \int _{\Omega ^+}|F^+|^2 dx+\int _{\Omega ^-_R}|F^-|^2 dx \end{aligned}$$(112)by the norm of \(f^0\) and compact terms. We define function spaces \({\widetilde{{{\mathcal {H}}}}}^{2,\pm }_{k_\pm }\) of Maxwell fields \(F^\pm = F^\pm _1+F^\pm _2\) in \(\Omega ^\pm \) (\(F^-\) satisfying the radiation condition in \(\Omega ^-\)). The norm of \(F^+\in \widetilde{{\mathcal {H}}}^{2,+}_{k_+}\) is \(\Vert F^+\Vert _{L_2(\Omega ^+)}\) and the norm of \(F^-\in {\widetilde{{{\mathcal {H}}}}}^{2,-}_{k_-}\) is \(\Vert F^-\Vert _{L_2(\Omega ^-_R)}\), with a fixed large R. Furthermore we define the closed subspace \({{\mathcal {H}}}^2_{k_-}\subset {{\mathcal {H}}}_3\), consisting of \(f= f_{1{\scriptstyle \top }}+ f_{1{\scriptstyle \perp }}+ f_{2{\scriptstyle \top }}+ f_{2{\scriptstyle \perp }}\), with \(f_{1{\scriptstyle \top }}, f_{1{\scriptstyle \perp }}\) being tangential and normal vector fields and \(f_{2{\scriptstyle \top }}, f_{2{\scriptstyle \perp }}\) being tangential and normal bivector fields, satisfying the constraints \(d_{\scriptstyle \top }f_{1{\scriptstyle \top }}= ik_- f_{2{\scriptstyle \top }}\) and \(\delta _{\scriptstyle \top }(\nu \mathbin {\lrcorner }f_{2{\scriptstyle \perp }})= -ik_- \nu \mathbin {\lrcorner }f_{1{\scriptstyle \perp }}\). The traces of incoming Maxwell fields \(F^0\) all belong to \({{\mathcal {H}}}^2_{k_-}\).
The transmission problem (107) defines a bounded linear map
$$\begin{aligned} T_{{\hat{\epsilon }}, k_+, k_-}: {\widetilde{{{\mathcal {H}}}}}_{k_+}^{2,+} \oplus {\widetilde{{{\mathcal {H}}}}}_{k_-}^{2,-}\rightarrow {{\mathcal {H}}}^2_{k_-}. \end{aligned}$$(113)Like in Proposition 6.2, we can continuously perturb the operator and spaces to the case \(k_+=k_-\) and that \({\hat{\epsilon }}=1\), and, given (ii) below, conclude by Fredholm perturbation theory that \(T_{{\hat{\epsilon }}, k_+, k_-}\) is a Fredholm operator with index zero and hence an isomorphism whenever (111) holds.
-
(ii)
To establish the estimate of (112), in \(\Omega _R^-\) we construct potentials such that
$$\begin{aligned} {\left\{ \begin{array}{ll} F^-_2= \nabla \mathbin {\scriptstyle {\wedge }}U^-_1,\\ F^-_1= \nabla \mathbin {\scriptstyle {\wedge }}U^-_0+ ik_- U^-_1, \end{array}\right. } \end{aligned}$$(114)and in \(\Omega ^+\) we construct potentials such that
$$\begin{aligned} {\left\{ \begin{array}{ll} F^+_2= \nabla \mathbin {\lrcorner }U^+_3+ ik_+ U^+_2,\\ F^+_1= \nabla \mathbin {\lrcorner }U^+_2, \end{array}\right. } \end{aligned}$$(115)where subscript j refers to a \(\wedge ^j {{\mathbf {R}}}^3\) valued field. (In traditional terminology, \(U^-_0\) is the scalar potential and \(U^-_1\) is the vector potential to the electromagnetic field \(F^-\).) The existence and estimates of such potentials follow from the Hodge decompositions in [6], by writing out the homogeneous parts of the multivector fields.
We consider first \(F^\pm _1\), where we pair the jump relation for \(\nu \mathbin {\scriptstyle {\wedge }}f_1^\pm \) with \(u_2^+\). From the jump relation for \(\nu \mathbin {\lrcorner }f_2^\pm \) and Maxwell’s equations, we obtain
$$\begin{aligned} \nu \mathbin {\lrcorner }f^+_1={\hat{\epsilon }}^{-1} \nu \mathbin {\lrcorner }(f^-_1+f^0_1), \end{aligned}$$(116)which we pair with \(u_0^-\). We get
$$\begin{aligned} {\left\{ \begin{array}{ll} \int _{\partial \Omega }\langle \nu \mathbin {\scriptstyle {\wedge }}f_1^+,u_2^+ \rangle d\sigma = \int _{\partial \Omega }(\langle \nu \mathbin {\scriptstyle {\wedge }}f_1^-,u_2^+ \rangle d\sigma + \langle \nu \mathbin {\scriptstyle {\wedge }}f_1^0,u_2^+ \rangle )d\sigma ,\\ \int _{\partial \Omega }\langle u_0^-,\nu \mathbin {\lrcorner }f_1^+ \rangle d\sigma =\overline{1/\hat{\epsilon }} \int _{\partial \Omega }(\langle u_0^-,\nu \mathbin {\lrcorner }f_1^- \rangle + \langle u_0^-,\nu \mathbin {\lrcorner }f_1^0 \rangle ) d\sigma . \end{array}\right. } \end{aligned}$$(117)Next consider \(F^\pm _2\), where we pair the jump relation for \(\nu \mathbin {\lrcorner }f_2^\pm \) with \(u_1^-\). From the jump relation for \(\nu \mathbin {\scriptstyle {\wedge }}f_1^\pm \) and Maxwell’s equations, we obtain
$$\begin{aligned} \nu \mathbin {\scriptstyle {\wedge }}f^+_2={\hat{k}}^{-1} \nu \mathbin {\scriptstyle {\wedge }}(f^-_2+f^0_2), \end{aligned}$$(118)which we pair with \(u_3^+\). We get
$$\begin{aligned} {\left\{ \begin{array}{ll} \int _{\partial \Omega }\langle \nu \mathbin {\scriptstyle {\wedge }}f_2^+,u_3^+ \rangle d\sigma =(1/{\hat{k}}) \int _{\partial \Omega }(\langle \nu \mathbin {\scriptstyle {\wedge }}f_2^-,u_3^+ \rangle + \langle \nu \mathbin {\scriptstyle {\wedge }}f_2^0,u_3^+ \rangle )d\sigma ,\\ \int _{\partial \Omega }\langle u_1^-,\nu \mathbin {\lrcorner }f_2^+ \rangle d\sigma =\overline{{\hat{k}}/{\hat{\epsilon }}} \int _{\partial \Omega }(\langle u_1^-,\nu \mathbin {\lrcorner }f_2^- \rangle + \langle u_1^-,\nu \mathbin {\lrcorner }f_2^0 \rangle ) d\sigma . \end{array}\right. } \end{aligned}$$(119)Applying the general Stokes theorem, see [37, Sec. 7.3], to (117) and (119) and summing the so obtained estimates yield an estimate of (112), similar to the proof of Proposition 6.2. A main idea is that the potentials \(U^\pm _j\) depend compactly on the fields \(F^\pm _j\), and for details of the estimates we refer to [4, Lem. 4.9, 4.17].
\(\square \)
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 Sect. 7.
Proposition 8.4
Consider the Dirac transmission problem (58) for a bounded Lipschitz domain \(\Omega ^+\subset {{\mathbf {R}}}^3\) with connected exterior \(\Omega ^-\). Let
with parameters \(\alpha ,\beta ,\gamma \in {{\mathbf {C}}}\setminus \{0\}\). Assume that
recalling that \({\hat{k}}= k_+/k_-\) and Definition 2.1. Then the operator \(B_3:E_{k_+}^+{{\mathcal {H}}}_3\oplus E_{k_-}^-{{\mathcal {H}}}_3 \rightarrow {{\mathcal {H}}}_3\) given by
is invertible, where we as in Sect. 5 denote the spaces of traces of solutions on \(\partial \Omega \) by \(E_{k_\pm }^\pm {{\mathcal {H}}}_3\).
Proof
The proof is similar to that of Proposition 7.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.
-
(i)
The scalar functions \(F_0^\pm \) solve the Helmholtz equation as a consequence of \(F^\pm \) solving the Dirac equation, which also shows
$$\begin{aligned} \nabla F^\pm _0= ik_\pm F^\pm _1-\nabla \mathbin {\lrcorner }F^\pm _2. \end{aligned}$$(124)From this and (120), we conclude that \(\partial _\nu f_0^+=\beta \partial _\nu ({\hat{k}} f^-_0/(\alpha \beta ))\) and \(f_0^+= {\hat{k}} f^-_0/(\alpha \beta )\) and hence Proposition 6.1 with \({\hat{\epsilon }}= \beta \) shows that \(F_0^\pm =0\).
-
(ii)
Writing \(F^\pm _3= U^\pm _3 e_{123}\), as in (i) the scalar functions \(U^\pm _3\) solve the Helmholtz equation. From the Dirac equation we have
$$\begin{aligned} \nabla \mathbin {\lrcorner }F^\pm _3= ik_\pm F^\pm _2-\nabla \mathbin {\scriptstyle {\wedge }}F^\pm _1. \end{aligned}$$(125)From this and (120), we conclude that \(\partial _\nu u_3^+=\gamma \partial _\nu ( u^-_3/\gamma )\) and \(u_3^+= u^-_3/\gamma \) and hence Proposition 6.1 with \({\hat{\epsilon }}= \gamma \) shows that \(F_3^\pm =0=U^\pm _3\).
-
(iii)
From (i) and (ii) we conclude that \(F^\pm \) solve Maxwell’s equations (107) with \({\hat{\epsilon }}=\alpha \), and Proposition 8.2 applies to show that \(F^\pm =0\).
To show existence, it suffices by perturbation theory for Fredholm operators to prove an estimate
$$\begin{aligned} \Vert f^+\Vert _{{{\mathcal {H}}}_3}+ \Vert f^-\Vert _{{{\mathcal {H}}}_3}\lesssim \Vert f^0\Vert _{{{\mathcal {H}}}_3}. \end{aligned}$$(126)This follows similarly to steps (i)–(iii) by instead appealing to Propositions 6.2 and 8.3.
\(\square \)
Proposition 8.5
Assume the hypothesis of Proposition 8.4, and further assume that
with parameters satisfying
Then the operator
is invertible on \({{\mathcal {H}}}_3\), for any \(r\in {{\mathbf {C}}}\setminus \{0\}\).
Proof
This follows by duality from Proposition 8.4, entirely analogously to the proof of Proposition 7.3. \(\square \)
Consider the 3D Dirac integral equation
involving the operator from (129), with right-hand side specified by the jump condition in (58) and auxiliary density \({\tilde{h}}\) which we precondition as \({\tilde{h}}= P'h\) in the proof of Theorem 2.3. Similar to what we did for the 2D twin (98), we optimize (130) by choosing the parameters \(r,\beta ,\gamma ,\alpha ',\beta ',\gamma '\). Recall that for EM fields \(\alpha ={\hat{\epsilon }}\). We therefore consider \(\alpha \) as having a prescribed value. Writing (129) as in (99), we have in \({{\mathbf {R}}}^3\) that
Modulo operators of the form \(K^v_{k_\pm }-K^v\), \(S^a_{k_+}-ik_-\hat{k}S^a\) and \(S^a_{k_-}-ik_-S^a\), the integral operator \(T=rE_{k_+} M'-M E_{k_-}\) from (99) is the entry-wise product of
and \(E_k\), with diagonal and off-diagonal \(2\times 2\) blocks replaced by \(K^v\) and \(ik_- S^a\) respectively. It is the diagonal \(4\times 4\) blocks in (132) which are our main concern, within which we note that the diagonal \(2\times 2\) blocks are weakly singular operators on smooth domains.
-
Our first choice is to set \(r=1/{\hat{k}}\), \(\alpha \beta ={\hat{k}}^2\) and \(\alpha '\beta '=1/{\hat{k}}^2\). This gives cancellation in the (1:2,3:4) and (7:8,5:6) size \(2\times 2\) blocks of T, which for a smooth domain \(\Omega \) yields a nilpotent operator T with essential spectrum \(\{0\}\), if \({{\mathcal {H}}}_3\) is replaced by a function space of fixed regularity.
-
It remains to choose \(\gamma , \alpha ',\gamma '\). The choice of \(\alpha '\) concerns the diagonal elements (1,1), (7,7) and (8,8). We set \(\alpha '= 1/{\hat{k}}\), and so \(\beta '=1/{\hat{k}}\). The choices of \(\gamma , \gamma '\) concern the diagonal elements (5,5) and (2,2). We set \(\gamma = {\hat{k}}/|{\hat{k}}|\) and \(\gamma '= \overline{{\hat{k}}}/|{\hat{k}}|\).
These choices make \(\beta /{\hat{k}}=(\alpha /{\hat{k}})^{-1}\), \(\alpha '{\hat{k}}=\beta ' {\hat{k}}=1\), \(\gamma ' {\hat{k}}>0\) and \(\gamma /{\hat{k}}>0\). Therefore, if \(\alpha , {\hat{k}}^2/\alpha \in {{\mathbf {C}}}\setminus (-\infty ,0]\) and \(\alpha /{\hat{k}}\in {\mathrm{WP}}\,(k_-,k_+)\), then Propositions 8.4 and 8.5 guarantee invertibility of \(rE_{k_+}^+M'+M E_{k_-}^-\). Note that for non-magnetic materials \({\hat{k}}^2/\alpha =\beta ={\hat{\mu }}=1\). Similar to the situation in \({{\mathbf {R}}}^2\), we also obtain good invertibility properties of the (1,1), (2,2), (5,5) diagonal elements, as well as the (7:8,7:8) diagonal block, on any Lipschitz domain in this way. Furthermore, with \({\hat{\mu }}=1\) we also obtain good invertibility properties for the (3,3) and (4,4) diagonal elements, so that it is only the (6,6) element which we do not control in the sense that we may hit its essential spectrum for \({\hat{\epsilon }} <0\).
To summarize, for solving the Maxwell transmission problem (14), we have obtained the 3D Dirac integral equation
with
9 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\times 4\) system of integral equations [23, Eq. (12.4)], along with its field representation formulas [23, Eqs. (12.2) and (12.3)]. For simplicity, we refer to the system (10) as “Dirac” and to the system [23, Eq. (12.4)] as “HK 4-dens”. The reason for comparing with “HK 4-dens” is that this system, just like “Dirac”, has a \(8\times 8\) counterpart in 3D which applies to Maxwell’s equations.
We also compare to a state-of-the-art \(2\times 2\) system of integral equations [23, Eq. (12.7)] of Kleinman–Martin type [27] and to the 2D version of the classic Müller system [32, p. 319], which also is a \(2\times 2\) system [23, Sec. 14.1], 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”.
We stress, at this point, that the purpose of the numerical tests is to verify that “Dirac” has the properties claimed in the theoretical sections of this paper. We perform these tests in 2D simply because we have not yet access to a high-order accurate solver for “Dirac” in 3D. For example, we monitor condition numbers under wavenumber sweeps in order to detect eigenwavenumbers. For scattering problems, we investigate achievable field accuracy and the convergence of iterative solvers. We are not trying to show that “Dirac” is more efficient than “best KM-type” in 2D, because it is not. Rather, we demonstrate that “Dirac” is almost as efficient as “best KM-type” in 2D and often more efficient than “HK 4-dens”. This is of interest because “Dirac” is also applicable in 3D, while “best KM-type” is not.
All systems of integral equations are discretized using Nyström discretization and underlying composite 16-point Gauss–Legendre quadrature. The reason for using Nyström discretization, rather than the more common Galerkin discretization (BEM or MoM), has to do with accuracy and speed. High achievable accuracy and fast execution is, in our opinion, more easily obtained in a Nyström setting. This applies not only in 2D, but also to scattering problems involving rotationally symmetric interfaces \(\partial \Omega \) in 3D. See, for example, [21, Sec. VII.A], where double-precision results, obtained with a high-order Fourier–Nyström method, agree to almost 14 digits with semi-analytical solutions given by Mie theory for an electromagnetic transmission eigenproblem on the unit sphere. See also [13, 29] for other recent examples where authors prefer Nyström schemes over Galerkin for electromagnetic scattering. The discussion of Nyström versus Galerkin discretization of integral equations in computational electromagnetics in [29, Sec. 1] is very informative.
On smooth \(\partial \Omega \) we use a variant of the “Nyström scheme B” in [19]. In the presense of singular boundary points on \(\partial \Omega \) 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) [18]. We note, in passing, that RCIP is applicable also to Fourier–Nyström discretization near sharp edges and singular boundary points on rotationally symmetric \(\partial \Omega \) in 3D [20, 25]. Accurate evaluation of singular operators on \(\partial \Omega \) and accurate field evaluation of layer potentials at field points close to \(\partial \Omega \) are accomplished using panel-wise product integration. See [22, Sec. 4], and references therein, for details. See also [22, Sec. 9.1] for a thorough test of the implementation of the integral operators needed in “Best KM-type” and “2D Müller” via Calderón identites and see [19, Sec. 9] for a high-wavenumber test of the implementation of the corresponding layer potentials needed for field evaluations. Large discretized linear systems are solved iteratively using GMRES (without restart).
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.
Condition numbers of the operators in “Dirac”, in “HK 4-dens”, in “best KM-type”, and in “2D Müller”: a the starfish-like interface \(\partial \Omega \); b the positive dielectric case; c the plasmonic case; d the reverse plasmonic case
9.1 The Operators
We compute condition numbers of the discretized system matrices in the systems under study. The main purpose is to detect false eigenwavenumbers. Another purpose is to compare the conditions number of the system matrices with each other. The interface \(\partial \Omega \) is the smooth starfish-like curve [22, Eq. (92)], shown in Fig. 3a and originally suggested for scattering problems in [16]. A number of 976 discretization points are placed on \(\partial \Omega \). We study three cases of jumps \({\hat{k}}\) where we vary \(k_-\) or \(k_+\):
-
The positive dielectric case. The exterior wavenumber is positive real, \(0<k_-\le 20\), and \({\hat{k}}=k_+/k_-= 1.5\) so that \(0<k_+\le 30\) and \({\hat{\epsilon }}=2.25\). This corresponds to the lower left corner point in Fig. 1, where Theorem 2.2 guarantees that no true eigenwavenumbers exist, as well as no false eigenwavenumbers for “Dirac”.
Figure 3b 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_-\le 20\), but \({\hat{k}}=k_+/k_-= i\sqrt{1.1838}\) so that \(0<k_+/i\le 20\sqrt{1.1838}\), \(k_+\) is imaginary, and \({\hat{\epsilon }}=-1.1838\) is rather close to the essential spectrum \(\{-1\}\). This corresponds to the left middle corner point in Fig. 1, where Theorem 2.2 guarantees that no true eigenwavenumbers exist, as well as no false eigenwavenumbers for “Dirac”.
Figure 3c 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_+\le 20\), and \({\hat{k}}=k_+/k_-= (i\sqrt{1.1838})^{-1}\) so that \(0<k_-/i\le 20\sqrt{1.1838}\), \(k_-\) is imaginary, and \({\hat{\epsilon }}=-1/1.1838\). This corresponds to the lower middle corner point in Fig. 1, which does not belong to the hexagon, and indeed Fig. 3d 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”.
9.2 Field Computations
Field computation in the positive dielectric case; a the scattered fields \({{\,\mathrm{Re}\,}}U^\pm \); b \(\log _{10}\) of the estimated absolute error using “HK 4-dens”; c same using “Dirac”; d same using “2D Müller”
We solve the Dirac system (10) for h, compute interior and exterior fields \(U^\pm \) via (12), and compare with results from the other systems. Gradient fields \(\nabla U^\pm \) are not computed, but we remark that the representation formula (13) for \(\nabla U^\pm \) uses layer potentials with the same type of (near-logarithmic and near-Cauchy) singular kernels as the representation formula (12). The \(2\times 2\) systems, on the other hand, have accompanying representation formulas for \(\nabla U^\pm \) with layer potentials that contain near-hypersingular kernels. Evaluation of these layer potentials at field points near \(\partial \Omega \) may cause additional loss of accuracy.
The curve \(\partial \Omega \) is chosen as the boundary of the one-corner drop-like object [22, Eq. (92)], discernible in Fig. 4a. 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)/\sqrt{2}}\), and 800 discretization points are placed on \(\partial \Omega \). When \({\hat{\epsilon }}=-1.1838\), then \(\pm ({\hat{\epsilon }}+1)/({\hat{\epsilon }}-1)\) is in the essential \(H^{1/2}(\partial \Omega )\)-spectrum of (1) and a homotopy-based numerical procedure is adopted where \({\hat{\epsilon }}=-1.1838\) is approached from above in the complex plane [17, Sec. 6.3].
Field computation in the plasmonic case; a the scattered fields \({{\,\mathrm{Re}\,}}U^\pm \); b \(\log _{10}\) of the estimated absolute error using “HK 4-dens”; c same using “Dirac”; d same using “best KM-type”
Field computation in the reverse plasmonic case; a the scattered fields \({{\,\mathrm{Re}\,}}U^\pm \); b \(\log _{10}\) of the estimated absolute error using “HK 4-dens”; c same using “Dirac”; d same using “2D Müller”
-
Figure 4 covers the positive dielectric case. The real parts of the scattered fields \(U^\pm \) are shown in Fig. 4a. There are propagating waves in both \(\Omega ^+\) and \(\Omega ^-\). The remaining images show \(\log _{10}\) of the estimated absolute pointwise error in \({{\,\mathrm{Re}\,}}U^\pm \), computed from the different systems. “Dirac” loses one digit of accuracy in some regions near \(\partial \Omega \) compared to the other systems. Most likely, this is because (12) does not exploit null-fields in the near-field evaluation. See [23, 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.
-
Fig. 5 covers the plasmonic case and is organized as Fig. 4. There are propagating waves in \(\Omega ^-\), exponentially decaying waves into \(\Omega ^+\), and a surface plasmon wave along \(\partial \Omega \). “Dirac” here performs almost on par with “best KM-type” and gives two more accurate digits than “HK 4-dens”.
The number of GMRES iterations needed are 266, 173, and 143 for “HK 4-dens”, “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 6 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 \(\Omega ^+\) and exponentially decaying waves into \(\Omega ^-\). “Dirac” performs almost on par with “2D Müller” and gives two more accurate digits than than “HK 4-dens”.
9.3 Densities and Function Spaces
We show asymptotics of the density \(h=[h_1\; h_2\; h_3\; h_4]^\mathrm{T}\) obtained from the Dirac integral equation (10). The computations relate to the three examples for the drop-like object in Sect. 9.2.
The densities \(h_1,h_2,h_3,h_4\) as functions of the arc length distance to the corner vertex in the positive dielectric case: a along the lower part of \(\partial \Omega \); b along the upper part of \(\partial \Omega \)
Same as Fig. 7, but for the plasmonic case
Same as Fig. 7, but for the reverse plasmonic case
-
The positive dielectric case: Here the hypothesis on \({\hat{\epsilon }}\) and \({\hat{k}}\) in Theorem 2.2 is satisfied. Then h belongs to the energy function space \({{\mathcal {H}}}_2\) from (63), meaning that \(h_1,h_2\in H^{1/2}(\partial \Omega )\) and \(h_3,h_4\in H^{-1/2}(\partial \Omega )\). The result is shown in Fig. 7, where indeed \(h_1\) and \(h_2\) are seen to be continuous at the corner vertex (note that \(h_2\approx 0\)). The only singular density is \(h_3\), which is related to that it is only the (3, 3) diagonal element in the \(4\times 4\) block-operator of (10) which we do not control by the choices of parameters \(r,\beta ,\alpha ',\beta '\) in Sect. 7. Using the automated eigenvalue method of [18, Sec. 14], the asymptotic behaviour of \(h_3\) near the corner is determined to be
$$\begin{aligned} h_3(t)\propto t^\eta , \end{aligned}$$(136)where \(\eta =-0.12319432456634\) and t is the arc length distance to the corner vertex. So \(h_3\) is in fact in \(L_2(\partial \Omega )\).
-
The plasmonic case: Here the hypothesis on \({\hat{\epsilon }}\) and \({\hat{k}}\) in Theorem 2.2 is not satisfied since \({\hat{\epsilon }}<0\) makes \(\pm ({\hat{\epsilon }}+1)/({\hat{\epsilon }}-1)\) hit the essential \(H^{-1/2}(\partial \Omega )\)-spectrum of (1). Nevertheless, the RCIP-accelerated Nyström scheme manages to produce the limit solution h shown in Fig. 8. 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\notin H^{-1/2}(\partial \Omega )\). More precisely, its asymptotics near the corner are as in (136) with \(\eta =-1.00000000000000-i1.57105873276994\). So \(h_3\in H^{-s}(\partial \Omega )\) for any \(s>1/2\).
-
The reverse plasmonic case: the results, shown in Fig. 9, are very similar to those of the plasmonic case. The asymptotics of \(h_3\) are as in (136) with \(\eta =-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 \({{\mathcal {H}}}_2\) from (63). 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\le s\le 1\). In the plasmonic cases, our computed densities h belong to the larger spaces \(s<1/2\). For \(0\le s<1\), the corresponding norms of the fields are weighted Sobolev norms using \(\int |\nabla U^\pm (x)|^2 \text {dist}(x)^{1-2s}dx\), where \(\text {dist}(x)\) denotes distance from x to \(\partial \Omega \), 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 [12]. The bounds on double layer potential operators for the endpoints \(s=0,1\) require harmonic analysis and the Coifman–McIntosh–Meyer theorem [10].
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
where \(\delta =\theta /\pi -1\) on a domain with a corner of opening angle \(\theta \).
The point of departure for the investigations reported on in this paper was the spin integral equation proposed in [36, 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 this spin integral equation makes use of non-diagonal matrices \(P,P',N, N'\) in (10) which mix \(H^{\pm 1/2}(\partial \Omega )\) and is limited to the function space \(L_2(\partial \Omega )^4\), which is different from \({{\mathcal {H}}}_2\). As we have seen, surface plasmon waves appear in the function space \({{\mathcal {H}}}_2\) or, in the case of pure meta materials, in the larger function spaces \(s<1/2\). The numerical algorithm used here fails for the spin integral equation when the “figure eight” of (137) is approached. When \({\hat{\epsilon }}= -1.1838\) is approached from above in the complex plane, this happens near \({\hat{\epsilon }}=-1.1838+i0.2168\).
References
Axelsson, A.: Oblique and normal transmission problems for Dirac operators with strongly Lipschitz interfaces. Commun. Partial Differ. Equ. 28(11–12), 1911–1941 (2003)
Axelsson, A.: Transmission problems for Dirac’s and Maxwell’s equations with Lipschitz interfaces. Ph.D. thesis, The Australian National University, 2003. Available at https://openresearch-repository.anu.edu.au/handle/1885/46056
Axelsson, A.: Transmission problems and boundary operator algebras. Integral Equ. Oper. Theory 50(2), 147–164 (2004)
Axelsson, A.: Transmission problems for Maxwell’s equations with weakly Lipschitz interfaces. Math. Methods Appl. Sci. 29(6), 665–714 (2006)
Axelsson, A., Grognard, R., Hogan, J., McIntosh, A.: Harmonic., analysis of Dirac operators on Lipschitz domains. In: Clifford analysis and its applications (Prague 2000), vol. 25 of NATO Sci. Ser. II Math. Phys. Chem. Kluwer Acad. Publ. Dordrecht 2001, 231–246 (2000)
Axelsson, A., McIntosh, A.: Hodge decompositions on weakly Lipschitz domains. In: Advances in analysis and geometry, Trends Math. Birkhäuser, Basel, pp. 3–29 (2004)
Buffa, A., Ciarlet, P.J.: On traces for functional spaces related to Maxwell’s equations. I. An integration by parts formula in Lipschitz polyhedra. Math. Methods Appl. Sci. 24(1), 9–30 (2001)
Buffa, A., Ciarlet, P.J.: On traces for functional spaces related to Maxwell’s equations. II. Hodge decompositions on the boundary of Lipschitz polyhedra and applications. Math. Methods Appl. Sci. 24(1), 31–48 (2001)
Buffa, A., Costabel, M., Sheen, D.: On traces for \(H(curl,\Omega )\) in Lipschitz domains. J. Math. Anal. Appl. 276(2), 845–867 (2002)
Coifman, R.R., McIntosh, A., Meyer, Y.: L’intégrale de Cauchy définit un opérateur borné sur \(L^{2}\) pour les courbes lipschitziennes. Ann. Math. 116(2), 361–387 (1982)
Colton, D., Kress, R.: Integral Equation Methods in Scattering Theory, 1st edn. Wiley, New York (1983)
Costabel, M.: Boundary integral operators on Lipschitz domains: elementary results. SIAM J. Math. Anal. 19(3), 613–626 (1988)
Epstein, C., Greengard, L., O’Neil, M.: A high-order wideband direct solver for electromagnetic scattering from bodies of revolution. J. Comput. Phys. 387, 205–229 (2019)
Fabes, E., Jodeit, M., Lewis, J.: Double layer potentials for domains with corners and edges. Indiana Univ. Math. 26(1), 95–114 (1977)
Ganesh, M., Hawkins, S.C., Volkov, D.: An all-frequency weakly-singular surface integral equation for electromagnetism in dielectric media: reformulation and well-posedness analysis. J. Math. Anal. Appl. 412(1), 277–300 (2014)
Hao, S., Barnett, A.H., Martinsson, P.G., Young, P.: High-order accurate methods for Nyström discretization of integral equations on smooth curves in the plane. Adv. Comput. Math. 40(1), 245–272 (2014)
Helsing, J.: The effective conductivity of arrays of squares: large random unit cells and extreme contrast ratios. J. Comput. Phys. 230(20), 7533–7547 (2011)
Helsing, J.: Solving integral equations on piecewise smooth boundaries using the RCIP method: a tutorial. arXiv:1207.6737v9 [physics.comp-ph] (revised 2018)
Helsing, J., Holst, A.: Variants of an explicit kernel-split panel-based Nyström discretization scheme for Helmholtz boundary value problems. Adv. Comput. Math. 41(3), 691–708 (2015)
Helsing, J., Karlsson, A.: Determination of normalized electric eigenfields in microwave cavities with sharp edges. J. Comput. Phys. 304, 465–486 (2016)
Helsing, J., Karlsson, A.: Resonances in axially symmetric dielectric objects. IEEE Trans. Microw. Theory Tech. 65(7), 2214–2227 (2017)
Helsing, J., Karlsson, A.: On a Helmholtz transmission problem in planar domains with corners. J. Comput. Phys. 371, 315–332 (2018)
Helsing, J., Karlsson, A.: An extended charge-current formulation of the electromagnetic transmission problem. SIAM J. Appl. Math. 80(2), 951–976 (2020)
Helsing, J., Karlsson, A., Rosén, A.: Comparison of integral equations for the Maxwell transmission problem with general permittivities. Preprint at arXiv:2007.12260 [physics.comp-ph]
Helsing, J., Perfekt, K.-M.: The spectra of harmonic layer potential operators on domains with rotationally symmetric conical points. J. Math. Pures Appl. 118, 235–287 (2018)
Hernández-Herrera, A.: Higher dimensional transmission problems for Dirac operators on Lipschitz domains. J. Math. Anal. Appl. 478(2), 499–525 (2019)
Kleinman, R., Martin, P.: On single integral equations for the transmission problem of acoustics. SIAM J. Appl. Math. 48(2), 307–325 (1988)
Kress, R., Roach, G.: Transmission problems for the Helmholtz equation. J. Math. Phys. 19(6), 1433–1437 (1978)
Lai, J., O’Neil, M.: An FFT-accelerated direct solver for electromagnetic scattering from penetrable axisymmetric objects. J. Comput. Phys. 390, 152–174 (2019)
Marmolejo-Olea, E., Mitrea, I., Mitrea, M., Shi, Q.: Transmission boundary problems for Dirac operators on Lipschitz domains and applications to Maxwell’s and Helmholtz’s equations. Trans. Am. Math. Soc. 364(8), 4369–4424 (2012)
McIntosh, A., Mitrea, M.: Clifford algebras and Maxwell’s equations in Lipschitz domains. Math. Methods Appl. Sci. 22(18), 1599–1620 (1999)
Müller, C.: Foundations of the mathematical theory of electromagnetic waves. Revised and enlarged translation from the German. Die Grundlehren der mathematischen Wissenschaften, Band 155. Springer, Berlin (1969)
Picard, R.: On the low frequency asymptotics in electromagnetic theory. J. Reine Angew. Math. 354, 50–73 (1984)
Picard, R.: On a structural observation in generalized electromagnetic theory. J. Math. Anal. Appl. 110(1), 247–264 (1985)
Rosén, A.: A spin integral equation for electromagnetic and acoustic scattering. Appl. Anal. 96(13), 2250–2266 (2017)
Rosén, A.: Boosting the Maxwell double layer potential using a right spin factor. Integral Equ. Oper. Theory 91, 3 (2019)
Rosén, A: Geometric multivector analysis. Birkhäuser Advanced Texts: Basler Lehrbücher. [Birkhäuser Advanced Texts: Basel Textbooks]. Birkhäuser, Cham (2019). From Grassmann to Dirac
Schulz, E., Hiptmair, R.: First-kind boundary integral equations for the Dirac operator in 3D Lipschitz domains. Preprint at arXiv:2012.11994 [math.AP]
Taskinen, M., Vänskä, S.: Current and charge integral equation formulations and Picard’s extended Maxwell system. IEEE Trans. Antennas Propag. 55(12), 3495–3503 (2007)
Taskinen, M., Ylä-Oijala, P.: Current and charge integral equation formulation. IEEE Trans. Antennas Propag. 54(1), 58–67 (2006)
Vico, F., Greengard, L., Ferrando, M.: Decoupled field integral equations for electromagnetic scattering from homogeneous penetrable obstacles. Commun. Part. Differ. Equ. 43(2), 159–184 (2018)
Weck, N.: Traces of differential forms on Lipschitz boundaries. Analysis (Munich) 24(2), 147–169 (2004)
Acknowledgements
We thank Anders Karlsson for many useful discussions.
Funding
Open access funding provided by University of Gothenburg.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This work was supported by the Swedish Research Council under contract 2015-03780.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Helsing, J., Rosén, A. Dirac Integral Equations for Dielectric and Plasmonic Scattering. Integr. Equ. Oper. Theory 93, 48 (2021). https://doi.org/10.1007/s00020-021-02657-1
Received:
Revised:
Published:
DOI: https://doi.org/10.1007/s00020-021-02657-1