Numerische Mathematik

Numer. Math. DOI 10.1007/s00211-014-0619-z

Boundary integral equation analysis on the sphere Felipe Vico · Leslie Greengard · Zydrunas Gimbutas

Received: 7 March 2013 © Springer-Verlag Berlin Heidelberg 2014

Abstract We present a systematic analysis of the integral operators of potential theory that arise when solving the Helmholtz or Maxwell equations in the exterior (or interior) of a sphere in the frequency domain. After obtaining expressions for the signatures of layer potentials in the spherical harmonic or vector spherical harmonic basis, we turn to a consideration of various integral equations that have been proposed in the literature for problems of acoustic and electromagnetic scattering. The selection of certain parameters in “combined field” and “Calderon-preconditioned” formulations is shown to have a significant impact on condition number, extending earlier work by Kress and others. Mathematics Subject Classification 35Q61 · 45B05

65R20 · 78M05 · 78A45 · 31B10 · 35Q60 ·

L. Greengard (B) · Z. Gimbutas New York University, New York, USA e-mail: [email protected] F. Vico Instituto de Telecomunicaciones y Aplicaciones Multimedia (ITEAM), Universitat Politècnica de València, Valencia, Spain e-mail: [email protected] Present Address Z. Gimbutas Information Technology Laboratory, National Institute of Standards and Technology, 325 Broadway, Mail stop 891.01, Boulder, CO 80305-3328, USA e-mail: [email protected]

123

F. Vico et al.

1 Introduction A number of problems in physics and engineering require the calculation of exterior acoustic and electromagnetic scattering, governed by the Helmholtz or Maxwell equations. Integral equation methods are natural candidates in such cases, because they discretize the boundary of the scatterer alone and impose the outgoing radiation condition at infinity by construction. We will restrict our attention in this paper to scattering from the unit sphere in the frequency domain, which is an extremely wellstudied and well-understood problem. Our purpose is not to develop a new method specific to spherical geometries, but to present a framework for easily studying the behavior of any integral equation formulation by analyzing its performance in the simplest possible setting. To do so, we systematically compute the action of a fairly comprehensive collection of integral operators on scalar and tangential vector fields (i.e. layer potentials) when expanded in spherical harmonics. These signatures permit the rapid analysis of invertibility and conditioning for virtually any proposed integral equation (or preconditioner). From both a theoretical and practical perspective, it is desirable to develop integral representations that lead to non-singular Fredholm equations of the second kind, with as low a condition number and as compact a spectrum as possible. Optimization of these latter two properties on a sphere can, hopefully, provide good heuristics for the general case. Our work is related to that of Kress [15], Hsiao and Kleinman [12], Contopanagos et al. [10], Bruno et al. [5], and Boubendir and Turc [4] who performed such an analysis for specific integral equations on the sphere, and to Chandler-Wilde et al. [6,7], who have carried out a remarkable analysis for more general geometries. For background material on integral equations and scattering, there are numerous excellent textbooks, such as [8,13,16,18,20]. The paper is structured as follows. We first introduce the necessary mathematical background on special functions (spherical harmonics, Bessel and Hankel functions, etc.), followed by a derivation of the diagonal form of a library of layer potentials. We then use these results to analyze various classical and recently proposed integral equations. Some of these equations involve coupling constants, and we study the optimal selection of these constants, which serve as tuning parameters in seeking to improve the conditioning and spectral properties of the formulation. 2 Special functions and partial wave expansions We first describe the classical separation of variables solution for Helmholtz and Maxwell potentials in spherical coordinates (also known as partial wave expansions in the former case and Lorenz–Debye–Mie expansions in the latter case) [1,8,13,20]. Letting Pn (x) denote the standard Legendre polynomial of degree n, the spherical harmonic of degree n and order m is given by  Ynm (θ, φ) =

123

2n + 1 4π



(n − |m|)! |m| P (cos θ )eimφ , (n + |m|)! n

(1)

Boundary integral equation analysis on the sphere

where the associated Legendre functions Pnm can be defined by the Rodrigues’ formula Pnm (x) = (−1)m (1 − x 2 )m/2

dm Pn (x). dxm

The functions Ynm are orthonormal with respect to the L 2 inner product on the unit sphere 2π π f (θ, φ)g(θ, φ) sin θ dθ dφ.

 f, g = 0

(2)

0 (1)

Definition 1 The spherical Bessel and Hankel functions jn (z), h n (z) are defined in terms of the usual Bessel and Hankel functions via  π Jn+1/2 (z) jn (z) = 2z  π (1) h (1) (z) = (z) H n 2z n+1/2 In particular, (1)

h 0 (z) =

ei z . iz

It is well-known that, in the interior of the unit sphere, any solution to the Helmholtz equation u(x) + k 2 u(x) = 0

(3)

can be expressed as u(x) =

∞  n 

i anm jn (kr )Ynm (θ, φ),

(4)

n=0 m=−n

where (r, θ, φ) denote the spherical coordinates of x, for some sequence of coeffii }. In the exterior of the unit sphere, any solution satisfying the Sommerfeld cients {anm radiation condition ∂u − iku = o(r −1 ), r := |x| → ∞ ∂r can be expressed as u(x) =

∞  n 

o anm h n (kr )Ynm (θ, φ)

(5)

n=0 m=−n

123

F. Vico et al. o }. When the context is clear, we will use the for some sequence of coefficients {anm shorthand



↔

m,n

n ∞  

.

n=0 m=−n

Finally, in Sect. 4, we will need a basis in which to expand tangential vector fields on the unit sphere. For this, we define Unm (θ, φ) = √

1 ∇s Ynm (θ, φ), Xnm (θ, φ) = rˆ × Unm (θ, φ), n(n + 1)

(6)

where ∇s is the surface gradient ∇s u = θˆ

1 ∂u ∂u + φˆ ∂θ sin θ ∂φ

(7)

and rˆ denotes the unit outward normal. The functions Unm and Xnm are known as vector spherical harmonics. They are orthonormal with respect to the inner product 2π π F, G =

F(θ, φ) · G(θ, φ) sin θ dθ dφ. 0

(8)

0

For fields V that are divergence-free and satisfy the vector Helmholtz equation V(x) + k 2 V(x) = 0,

(9)

the Lorenz–Debye–Mie analog of (4) for the interior of the unit sphere is   1 i i V(x) = √ anm ∇ × ∇ × (ˆrJn (kr )Ynm ) + bnm ∇ × (ˆrJn (kr )Ynm ) , n(n + 1) m,n (10) where Jn are the Riccati–Bessel functions defined by Jn (x) = x jn (x).

(11)

Assuming V satisfies the (outgoing) Sommerfeld–Silver–Müller radiation condition,

  1 lim ∇ × V × rˆ − ikV = o , r →∞ r

123

(12)

Boundary integral equation analysis on the sphere

V can be expressed in the exterior of the unit sphere as V(x) = √

 1 o o anm ∇ × ∇ × (ˆrHn (kr )Ynm ) + bnm ∇ × (ˆrHn (kr )Ynm ) n(n + 1) m,n (13)

where Hn are the Riccati–Hankel functions defined by Hn (x) = xh n (x).

(14)

3 Layer potentials: the scalar case In this section, we derive the diagonal forms of single and double layer potentials and their normal derivatives. These layer potentials, Sk [σ ](x) and Dk [σ ](x), respectively, are defined for points x off the unit sphere as  Sk [σ ](x) =

gk (x − y)σ (y) d A y |y|=1

 Dk [σ ](x) = |y|=1

∂gk (x − y)σ (y) d A y . ∂n y

(15)

Here, gk is the outgoing free-space Green’s function gk (x) =

eik|x| 4π |x|

(16)

and ∂n∂ y denotes the outward normal derivative at the point y. For a point y on the sphere, Sk is weakly singular and the limit is well-defined; Sk σ [y ] = lim Sk [σ ](x). x→y

(17)

For the double layer, the limit from the interior and exterior are different and denoted by Dk+ σ [y ] = lim Dk [σ ](x) x→y

x >1

Dk− σ [y ]

= lim Dk [σ ](x) x→y

(18)

x <1

123

F. Vico et al.

The normal derivatives of the layer potentials are given by the following expressions Sk + σ [y ]

∂ = lim x→y ∂n y

x >1

∂ Sk − σ [y ] = lim x→y ∂n y

x <1

∂ Dk [σ ](y ) = z ∂n y



 gk (x − y)σ (y) d A y |y|=1



gk (x − y)σ (y) d A y |y|=1

|y|=1

∂gk (y − y)σ (y) d A y ∂n y

(19)

Note that, for the normal derivatives, it is the single layer potential which takes on a different limit from the interior and exterior. (Dk is hypersingular and defined in the Hadamard finite part sense, but takes on the same limit from the interior and exterior.) These properties can be summarized and made more precise through the following well-known jump conditions [8,13,16,18,20].

Sk σ = ∂S σ k = ∂n

Dk σ = ∂ D σ k = ∂n

Sk+ σ − Sk− σ = 0 Sk + σ − Sk − σ = −σ,

(20)

Dk+ σ − Dk− σ = σ Dk + σ − Dk − σ = 0.

(21)

Here, [ f ] denotes the jump in the quantity f when crossing the boundary of the scatterer (the unit sphere).

3.1 Diagonal form of the single layer potential Let us now consider the single layer potential Sk σ induced by the source density σ on the surface of the unit sphere. From (4) and (5), we may write  Sk σ =



i m m,n anm jn (kr )Yn (θ, φ) o m m,n anm h n (kr )Yn (θ, φ)

for r < 1 for r > 1.

(22)

Likewise, we may expand the source density σ in terms of spherical harmonics σ (θ, φ) =

 m,n

123

σnm Ynm (θ, φ).

(23)

Boundary integral equation analysis on the sphere

Setting r = 1, and using the orthogonality of the spherical harmonics as well as (22), (23), and (20), we have o i anm h n (k) − anm jn (k) = 0 o i kh n (k) − anm k jn (k) = −σnm anm

(24)

Solving this 2 × 2 linear system yields −σnm jn (k) = ik jn (k)σnm k( jn (k)h n (k) − jn (k)h n (k)) −σnm h n (k) = ik h n (k)σnm = k( jn (k)h n (k) − jn (k)h n (k))

o anm = i anm

(25)

The denominator was simplified using the Wronskian identity [1] jn (x)h n (x) − jn (x)h n (x) =

i . x2

(26)

Substituting (25) in (22), we obtain the diagonal form of the operator Sk  Sk Ynm (θ, φ) = ik jn (k)h n (k)Ynm (θ, φ)

(27)

Taking the limit k → 0 and using the appropriate limits for the Bessel and Hankel functions, we get the corresponding expression for the single layer potential with the Green’s function for the Laplace equation  S0 Ynm (θ, φ) =

1 Y m (θ, φ) 2n + 1 n

(28)

Taking the normal derivative of (22) and using (25), we obtain the diagonal forms for the normal derivative of the single layer potential on either side of the surface  Sk + Ynm (θ, φ) = ik 2 jn (k)h n (k)Ynm (θ, φ) (29) and  Sk − Ynm (θ, φ) = ik 2 jn (k)h n (k)Ynm (θ, φ)

(30)

Letting k → 0, it is straightforward to derive the results  n+1 m S0 + Ynm (θ, φ) = − Y (θ, φ) 2n + 1 n

(31)

 S0 − Ynm (θ, φ) =

(32)

and n Y m (θ, φ). 2n + 1 n

123

F. Vico et al.

3.2 Diagonal form of the double layer potential For the double layer potential Dk σ induced by the source density σ , we write  Dk σ =



i m m,n anm jn (kr )Yn (θ, φ)

for r < 1

o m m,n anm h n (kr )Yn (θ, φ)

for r > 1.

(33)

Matching coefficients and using (21), we have o i h n (k) − anm jn (k) = σnm anm o i anm kh n (k) − anm k jn (k) = 0

(34)

or o anm = ik 2 jn (k)σnm i anm = ik 2 h n (k)σnm .

(35)

This yields the diagonal forms for the double layer potential  Dk+ Ynm (θ, φ) = ik 2 jn (k)h n (k)Ynm (θ, φ)  Dk− Ynm (θ, φ) = ik 2 h n (k) jn (k)Ynm (θ, φ).

(36)

For the Laplace equation (k = 0), we have  D0+ Ynm (θ, φ) =

n Y m (θ, φ) 2n + 1 n  n+1 m D0− Ynm (θ, φ) = − Y (θ, φ) 2n + 1 n

(37)

Taking the normal derivative of (33) and using (35), we obtain the diagonal form for Dk  Dk Ynm (θ, φ) = ik 3 jn (k)h n (k)Ynm (θ, φ).

(38)

Finally, taking care with the limiting procedure  −n(n + 1) m Y (θ, φ) D0 Ynm (θ, φ) = (2n + 1) n

(39)

Note that Sk + , Dk− and Sk − , Dk+ have the same signature (which can also be seen from duality arguments).

123

Boundary integral equation analysis on the sphere

4 Layer potentials: the vector case In this section, we derive expansions for the various integral operators which arise in solving the Maxwell equations. Unfortunately, the analysis here is somewhat more involved than for the scalar case, and there is not a single standard formalism in the literature. Most treatments, however, make use of the vector potential in the Lorenz gauge, for which the source density is a tangential vector field (the electric current J). In particular, if we write Maxwell’s equations in the form ∇ × H = −ikE ∇ × E = ikH.

(40)

it is standard to express the scattered electric and magnetic fields (E and H) as H = ∇ ×A −1 ∇ ×∇ ×A E= ik

(41) (42)

where  A[J](x) =

gk (x − y)J(y)d A y 

and  denotes the surface of the scatterer. With a slight abuse of notation, we will denote the vector potential A by A[J](x) = Sk [J](x) by analogy with the scalar case. It is straightforward to verify that the vector potential satisfies the following jump conditions:

n × ∇ × Sk (J) = n × ∇ × Sk (J)+ − n × ∇ × Sk (J)− = J

n × ∇ × ∇ × Sk (J) = n × ∇ × ∇ × Sk (J)+ − n × ∇ × ∇ × Sk (J)− = 0. (43) Here, n is the unit outward normal at some point on . When the scatterer (here, the unit sphere) is a perfect conductor, two homogeneous conditions must be satisfied on the boundary [13,20] n × Etot = 0, n · Htot = 0.

(44)

It is also well-known that, on the surface, we must have n × Htot = J, n · Etot = ρ,

(45)

123

F. Vico et al.

when the scattered fields are induced by a physical surface current J and a corresponding charge ρ with ∇s · J = ikρ,

(46)

where ∇s · denotes the surface divergence. Here, Etot (x) = Ein (x) + E(x) and Htot (x) = Hin (x) + H(x) denote the total electric and magnetic fields, Ein and Hin are known incoming fields and E, H are the scattered fields. The first condition in (44) yields the electric field integral equation (EFIE) n × ∇ × ∇ × Sk [J] = −n × Ein . −ik

(47)

The first condition in (45) yields the magnetic field integral equation (MFIE) n × Hin + n × ∇ × Sk [J]+ = J.

(48)

Using the jump relations (43), we may rewrite the MFIE in the form − n × ∇ × Sk [J]− = n × Hin .

(49)

We will denote by Mk the integral operator on the left-hand side of (49) and by Ek the integral operator on the left-hand side of (47). Suppose now that we expand the field V(x) = ∇ × Sk [J] in the interior and exterior of the unit sphere using (10) and (13):  V(x) =

√ 1 n(n+1) √ 1 n(n+1)



i m,n anm ∇ o m,n anm ∇



i ∇ × (ˆ × ∇ × (ˆrJn (kr )Ynm ) + bnm rJn (kr )Ynm )

×∇

o ∇ × (ˆrHn (kr )Ynm ) + bnm

× (ˆrHn (kr )Ynm )

if r < 1 if r > 1.

(50) Taking the curl of (50), we obtain  ∇ ×V=

√ 1 n(n+1) √ 1 n(n+1)

 

i 2 m,n anm k ∇ o 2 m,n anm k ∇

i ∇ × ∇ × (ˆ × (ˆrJn (kr )Ynm ) + bnm rJn (kr )Ynm )

if r < 1

o ∇ × (ˆrHn (kr )Ynm ) + bnm

if r > 1.

×∇

× (ˆrHn (kr )Ynm )

(51) Suppose now that expand the source J in terms of the orthogonal basis of vector spherical harmonics: J=



U m X m Jnm Un (θ, φ) + Jnm Xn (θ, φ).

(52)

m,n

We wish to compute the action of Mk , Ek and n × ∇ × Sk+ in this basis. This will require the following relations. (When restricting our attention to the sphere, we will generally use rˆ in place of the more general outward normal n.)

123

Boundary integral equation analysis on the sphere

 rˆ × ∇ × (ˆrJn (kr )Ynm ) = Jn (kr )∇s Ynm = Jn (kr )Unm n(n + 1)  rˆ × ∇ × ∇ × (ˆrJn (kr )Ynm ) = kJ n (kr )ˆr × ∇s Ynm = kJ n (kr )Xnm n(n + 1)  rˆ × ∇ × (ˆrHn (kr )Ynm ) = Hn (kr )∇s Ynm = Hn (kr )Unm n(n + 1)  rˆ × ∇ × ∇ × (ˆrHn (kr )Ynm ) = kH n (kr )ˆr × ∇s Ynm = kH n (kr )Xnm n(n + 1) (53) Lemma 1 The operators Ek , Mk and n × ∇ × Sk+ have the following diagonal forms:

  

n × ∇ × ∇ × Sk Unm Unm −J n (k)H n (k)Xnm = = Xnm Xnm Jn (k)Hn (k)Unm −ik

m

m −

 Un Un −iJn (k)H n (k)Unm = −n × ∇ × Sk = Mk Xnm Xnm iJ n (k)Hn (k)Xnm

m +  iJ n (k)Hn (k)Unm Un = n × ∇ × Sk . Xnm −iJn (k)H n (k)Xnm

Ek

(54)

Proof From (50), (51), (53) and the orthogonality of the vector spherical harmonics, we can write the jump relations (43) componentwise on the unit sphere.

rˆ × ∇ × Sk (J) = J ⇒

rˆ × ∇ × ∇ × Sk (J) = 0 ⇒

 

o kH (k) − a i kJ (k) = J X anm n n nm nm o H (k) − bi J (k) bnm n nm n

U = Jnm

o k 2 H (k) − a i k 2 J (k) = 0 anm n n nm o kH (k) − bi kJ (k) = 0 bnm n n nm

.

(55)

o , a i and bo , bi ; This takes the form of decoupled 2×2 systems of equations for anm nm nm nm

 

o kH (k) − a i kJ (k) = J X anm n n nm nm o k 2 H (k) − a i k 2 J (k) = 0 anm n n nm o H (k) − bi J (k) bnm n nm n

U = Jnm

o kH (k) − bi kJ (k) bnm n n nm

=0

.

(56)

Using the identity H n (x)Jn (x) − J n (x)Hn (x) = i

(57)

for Riccati–Bessel and Riccati–Hankel functions, we obtain Jn (k) ik U = i Jnm J n (k)

o X anm = Jnm o bnm

Hn (k) ik U = i Jnm H n (k).

i X anm = Jnm i bnm

(58)

123

F. Vico et al.

Combining (51), (53), and (58), it is straightforward to derive the diagonal form for Ek . Using (50), (53), and (58), it is straightforward to obtain the diagonal forms for

Mk and n × ∇ × Sk+ . Note that the signatures for the vector case are block diagonal (with 2 × 2 blocks), rather than being strictly diagonal. Note also that the limit of the Ek operator does not exist as k → 0 (one component goes to zero and one becomes infinite). Nevertheless, one can easily derive an asymptotic expansion for small values of k that permits stable numerical evaluation, namely

Ek

Unm Xnm



 ≈

n(n+1) m ik(2n+1) Xn ik − 2n+1 Unm

 (59)

For Mk and n × ∇ × Sk+ , the static limit (k → 0) does exist.

n × ∇ × S0

M0

Unm Xnm



= −n × ∇ × S0

Unm Xnm Unm Xnm

+ −

 =  =

n+1 m 2n+1 Un n m 2n+1 Xn n m 2n+1 Un n+1 m 2n+1 Xn

  .

(60)

4.1 Nonstandard operators For scattering from a perfect conductor, the MFIE and EFIE are the starting point in the classical treatments of electromagnetic scattering. Many implementations are based on a resonance-free (non-singular) linear combination of the two called the combined field integral equation (CFIE). In recent years, however, there has been a resurgence of interest in other integral representations that make use of other integral operators. In this section, we introduce the relevant operators and compute their diagonal forms. Lemma 2 The operators n × ∇ Sk and ∇ · Sk have the following diagonal forms:   (61) n × ∇ Sk Ynm = ik jn (k)h n (k) n(n + 1)Xnm   m  √ m Un −ik jn (k)h n (k) n(n + 1)Yn (62) ∇ · Sk = Xnm 0 Proof These formulas follow from (27), the identity ∇ · Sk (J) = Sk (∇s · J), (52) and (6).

123



Boundary integral equation analysis on the sphere

Note that n ×∇ Sk maps scalar vector fields to tangential vector fields and that ∇ · Sk maps tangential vector fields to scalar ones. Taking the limit k → 0, we obtain the following static limits: √  m n(n + 1) m n × ∇ S0 Yn = (63) Xn 2n + 1   √  m − n(n+1) m Un 2n+1 Yn = ∇ · S0 (64) Xnm 0 Lemma 3 The operator n × Sk has the diagonal form:   m  i((n + 1) jn (k)h n−1 (k) + n jn+1 (k)h n (k) − k jn+1 (k)h n−1 (k))Xnm Un n × Sk = Jn (k)Hn (k) m Xnm Un ik

(65) Proof The result follows from direct computation and a modest amount of algebra to obtain the signature in the form (65), which avoids catastrophic cancellation.

In the static limit, we obtain  n × S0

Unm Xnm



 =

4n 2 +4n+3 Xm 8n 3 +12n 2 −2n−3 n −1 m 2n+1 Un

 (66)

We will eventually have reason to consider the tangential components of the single layer operator acting on a normally oriented vector field, namely Sk (nρ). In order to find a suitable representation for this function, which is not necessarily divergencefree, we must supplement the representation used in (50) for the unit sphere. ⎧  i i ∇ × (ˆ rJn (kr )Ynm ) anm ∇ × ∇ × (ˆrJn (kr )Ynm ) + bnm ⎪ ⎪ m,n ⎪ i ∇( j (kr )Y m (θ, φ)) √ 1 ⎨ +cnm if r = |x| < 1 n n n(n+1) . Sk [ˆrρ](x) =   o m ) + bo ∇ × (ˆ m) ⎪ a ∇ × ∇ × (ˆ r H (kr )Y r H (kr )Y n n ⎪ nm n nm n m,n ⎪ ⎩ +co ∇(h (kr )Y m (θ, φ)) √ 1 if r = |x| > 1 n nm n n(n+1)

(67) o , ci can be computed from the fact that The values of cnm nm

∇ · Sk (nρ) = −Dk (ρ)

(68)

 o cnm = i jn (k)ρnm n(n + 1)  i = i h n (k)ρnm n(n + 1) cnm

(69)

and (36), yielding

123

F. Vico et al.

To find the remaining coefficients, we make use of the jump conditions

n × Sk (nρ) = n × Sk (nρ)+ − n × Sk (nρ)− = 0

n × ∇ × Sk (nρ) = n × ∇ × Sk (nρ)+ − n × ∇ × Sk (nρ)− = 0

(70)

Lemma 4 The tangential components of the single layer potential acting on a normally oriented vector field have the following diagonal form:    rˆ × Sk (ˆrYnm ) = i n(n + 1) jn (k)h n−1 (k) − jn+1 (k)h n (k) Xnm .

(71)

The normal component of the single layer potential acting on tangential vector fields has the following diagonal form:  rˆ · Sk

Unm Xnm

 =

  √ i n(n + 1)( jn (k)h n−1 (k) − jn+1 (k)h n (k))Ynm . 0

(72)

Proof The derivation of formula (71) is analogous to that of (54) in Lemma 1. Formula (72) follows from duality arguments (or can be computed directly).

Taking the static limit, √ 4 n(n + 1) Xm (2n − 1)(2n + 3)(2n + 1) n  √  m  4 n(n+1) Un Ynm (2n−1)(2n+3)(2n+1) = n · S0 . Xnm 0

n × S0 (nYnm ) =

(73) (74)

5 Summary of signatures In this section, we collect all the preceding results in tabular form. Here, S ∗ and D ∗ denote the principal value operators acting as functions on the boundary of the unit sphere. As pseudodifferential operators, S, S ∗ and D ∗ are of order -1, S + , S − , D + , D − are of order 0, and D is of order 1. In the vector case, suppose that we are given a block diagonal operator D such that  m  m  nm nm nm nm Un + anm d21 Xn . D anm Unm + bnm Xnm = anm d11 + bnm d12 + bnm d22 Then, if we expand the current in a spherical harmonic basis, the nm-block of D will be denoted ny

Dnm =

123

nm d nm d11 12 nm d nm d21 22

 .

(75)

Boundary integral equation analysis on the sphere

If D maps tangential vector fields to scalar fields,   D anm Unm + bnm Xnm = anm d1nm + bnm d2nm Ynm ,

(76)

then the nm-block will be denoted by  Dnm = d1nm d2nm .

(77)

Finally, if D maps scalar fields to tangential vector fields, D(Ynm ) = d1nm Unm + d2nm Xnm , then the nm-block will be denoted by

Dnm =

d1nm d2nm

 .

(78)

6 Integral operators in scattering theory We turn now to the analysis of various classical and more recently proposed integral equations for problems of acoustic and electromagnetic scattering. 6.1 Sound hard scatterers For “sound-hard” scatterers, the boundary condition to be enforced is the homogeneous Neumann condition ∂u tot = 0, ∂n where u tot denotes the total pressure. Assuming the incoming field is denoted by u in , that u tot = u in + u, and that u is represented as u(x) = Sk [σ ](x) + iηDk [S0 [σ ]] (x), we obtain the integral equation C0 [σ ] = Sk + [σ ] + iηDk [S0 [σ ]] = −

∂u in . ∂n

(79)

This equation was proposed by Panich [19] and the signature of C0 can be easily obtained from Tables 1 and 2: [S + ]

[D ]

[C0 ] [S0 ] λnm = λnmk + iηλnmk λnm .

(80)

123

F. Vico et al. Table 1 Signatures of Helmholtz layer potentials

Operator

Signature (λnm )

Sk

ik jn (k)h n (k)

Sk +

ik 2 jn (k)h n (k)

Sk − Sk ∗ Dk+ Dk− Dk∗ Dk

Table 2 Signatures of Laplace layer potentials

ik 2 jn (k)h n (k) ik 2 ( j (k)h (k) + n n 2 ik 2 jn (k)h n (k) ik 2 jn (k)h n (k) ik 2 ( j (k)h (k) + n n 2 ik 3 jn (k)h n (k)

jn (k)h n (k))

jn (k)h n (k))

Operator

Signature (λnm )

S0

1 2n+1 n+1 − 2n+1 n 2n+1 1 − 4n+2 n 2n+1 n+1 − 2n+1 1 − 4n+2 −n(n+1) (2n+1)

S0 + S0 − S0 ∗

D0+ D0− D0∗ D0

A variant with u represented as

u(x) = Sk [σ ](x) + iηDk S02 [σ ] (x),

(81)

is analyzed in [8], with a straightforward proof of invertibility for arbitrary geometries. Both are known to be non-resonant (invertible) Fredholm equations of the second kind. The signatures here permit the direct calculation of the eigenvalues, singular values and condition number of the integral equation. A variation of (79) with a different regularizing operator was suggested in [5]: Cik [σ ] = Sk + [σ ] + iηDk [Sik [σ ]] = −

∂u in . ∂n

(82)

It is also known to be second kind and non-resonant (see [5]) and its signature is given by: [S + ]

[D ]

Cik ] ik ] λ[nm = λnmk + iηλnmk λ[S nm

(83)

It is noted in [5] that the iteration count using Cik is noticeably better than using C0 and recommended as the integral equation of choice for sound-hard scattering.

123

Boundary integral equation analysis on the sphere Fig. 1 Condition number for the regularized combined source integral equation in solving the sound-hard (scalar Neumann) scattering problem using either C0 or Cik

Condition Number

400

CSIE Sik CSIE S0

300

200

100

0

0

20

40

60

80

100

k (wavenumber)

5

Eigenvalues CSIE Sik

0.4

Eigenvalues CSIE S0

4

0.2

Imag

Imag

3 2 1

0.2

0

1 2 10

0

0.4 8

6

4

2

0

1

Real

0.8

0.6 Real

0.4

0.2

Fig. 2 Eigenvalues of C0 and Cik for sound-hard scattering with k = 30. Note that the axes are different in the two plots

The condition number of the two integral equations is shown in Fig. 1, which verifies the assertion that the high-frequency behavior of Cik is much better than for C0 . We also carry out a more detailed study of the spectral properties of the two operators. Figure 2 shows that Cik has more compactly supported eigenvalues than does C S0 . (We omit a comparison with (81) whose spectral properties are less desirable.)

6.2 Electromagnetic scattering from a perfect conductor We turn now to the solution of Maxwell’s equations in the exterior of a perfect conductor. Colton and Kress [9] suggested the representation E(x) = ∇ × Sk [J](x) + iη∇ × ∇ × Sk [n × S02 (J)].

(84)

This is a variant of the combined source integral equation, where J is not a physical current, but an arbitrary tangential vector field on the surface of the scatterer. Using the homogeneous boundary condition on Etot from (44), one obtains

123

F. Vico et al. Table 3 Signatures of vector layer potentials Signature Dnm

Operator

 0 iJ n (k)Hn (k) 0 −iJn (k)Hn (k) 

0 −iJn (k)H n (k) 0 iJ n (k)Hn (k) 

0 Jn (k)Hn (k)  −Jn (k)Hn (k) √ 0 −ik jn (k)h n (k) n(n + 1) 0 

0√ ik j (k)h n (k) n(n + 1) ⎛ n Jn(k) Hn(k) ⎞ 0 ik  ⎜ i (n + 1) j (k)h ⎟ ⎜ ⎟ n n−1 (k)+ ⎝ ⎠ +n jn+1 (k)h n (k)− 0   −k jn+1 (k)h n−1 (k) 0   √ i n(n + 1) jn (k)h n−1 (k) − jn+1 (k)h n (k)  √    i n(n + 1) jn (k)h n−1 (k) − jn+1 (k)h n (k) 0

n × ∇ × Sk+ (J) Mk (J) = −n × ∇ × Sk− (J) Ek (J) =

n×∇×∇×Sk (J) −ik

∇ · Sk (J) n × ∇ Sk (ρ)

n × Sk (J)

n × Sk (nρ) n · Sk (J)

Table 4 Signatures of static vector layer potentials

Vector Laplace layer potentials Operator

Signature (λnm ) 

n × ∇ × S0 (J)+ M0 (J) Ek→0 (J) ∇ · S0 (J)



n 2n+1

0

 0 n 2n+1 

0

n+1 2n+1

ik 0 − 2n+1 n(n+1) 0  ik(2n+1)  √ − n(n+1) 0  2n+1  √

 n × S0 (J) 

n · S0 (J)

0



n × ∇ S0 (ρ)

n × S0 (nρ)

n+1 2n+1



0

n(n+1) 2n+1

0 4n 2 +4n+3 8n 3 +12n 2 −2n−3

0 √

−1  2n+1

0 

4 n(n+1)

 (2n−1)(2n+3)(2n+1) √

4 n(n+1) (2n−1)(2n+3)(2n+1) 0

+ K 2 inc SC 0 (J) = n × ∇ × Sk (J) + iηn × ∇ × ∇ × Sk (n × S0 (J)) = −n × E .



(85)

A proof of invertibility for any frequency k > 0 can be found in [9]. A variant (which uses less regularization) and also leads to a Fredholm equation of the second kind was proposed in [3]. Letting

123

Boundary integral equation analysis on the sphere Fig. 3 Condition number of regularized combined source integral equations for electromagnetic scattering from a perfect conductor

300

CSIE Sik CSIE S0

Condition Number

250 200 150 100 50 0

0

20

40

60

80

100

k (wavenumber)

E(x) = ∇ × Sk [J](x) + iη∇ × ∇ × Sk [n × S0 (J)]

(86)

leads to S0 (J) = n × ∇ × Sk+ (J) + iηn × ∇ × ∇ × Sk (n × S0 (J)) = −n × Einc (87) As in the scalar case, it was suggested and shown in [3] that replacing S0 with the Yukawa single layer potential Sik in (86) improves the condition number. The integral equation takes the form Sik (J) = n × ∇ × Sk+ (J) + iηn × ∇ × ∇ × Sk (n × Sik (J)) = −n × Einc

(88)

The signature of these equations can be obtained by making use of Table 3. [n×∇×Sk+ ]

S0 ] = Dnm D[nm

[n×∇×∇×Sk ] [n×S0 ] + iηDnm · Dnm

(89)

[n×∇×∇×Sk ] [n×Sik ] + iηDnm · Dnm .

(90)

and [n×∇×Sk+ ]

[Sik ] Dnm = Dnm

Figure 3 shows a remarkable difference in terms of conditioning at high frequencies between S0 and Sik . Figure 4 shows that the eigenvalues of Sik are much more clustered than those of S0 as well (Table 4). We should note that a variety of other local regularizers have been used successfully in the high frequency regime, such as [2]. 6.3 The generalized Debye representation for electromagnetic fields A recent formulation for electromagnetic scattering that avoids the problem of lowfrequency breakdown was proposed in [11]. For simply connected domains (such as

123

F. Vico et al. 2

0.6

Eigenvalues CSIE Sik 1

0.4

Imag

Imag

0 1 2 3 4

Eigenvalues CSIE S0 0

2

4

6

8

10

0.2 0 0.2 0.4 0.2

0.4

0.6

Real

0.8

1

1.2

Real

Fig. 4 Eigenvalues of S0 and Sik for electromagnetic scattering at k = 30

the sphere) it is based on the following representation for the electric and magnetic fields: ∇ × ∇ × Sk (J) −ik ∇ × ∇ × Sk (K) H = ∇ × Sk (J) + −ik E = −∇ × Sk (K) +

(91) (92)

where   −1 J = ik ∇s · −1 s ρ + n × ∇s · s σ ,

(93)

s denotes the surface Laplacian, and K = n × J. The boundary conditions imposed in this formulation are non-standard but derived directly from the two conditions in (44): − n · Hinc = n · H   S0 −∇s · (n × n × Einc ) = S0 (∇s · (n × n × E))

(94)

This is an invertible Fredholm equation of the second kind in terms of the generalized Debye sources ρ and σ . On the sphere, Eq. (94) is decoupled, so the first equation determines σ and the second equation determines ρ. That is, (94) takes the form λ1nm σnm = P(n · Hinc )nm   λ2nm ρnm = P(S0 ∇s · (n × n × Einc ) )nm

123

(95)

Boundary integral equation analysis on the sphere 200

Condition Number

Fig. 5 Condition number of the generalized Debye formulation for electromagnetic scattering with S0 and Sik regularization

gDEBYE Sik gDEBYE S0

150

100

50

0

0

20

40

60

80

100

k (wavenumber)

where P( f )nm denotes the projection of a function f on the sphere onto the spherical harmonic Ynm and [n×∇×Sk+ ]

[Ek ] + Dnm = (Jn (k) + iJ n (k))Hn (k) 11 12  [n×∇×S + ]  1 [S0 ] [Ek ] k Dnm = −ik ikλnm − Dnm (Jn (k) − iJ n (k))H n (k). 22 21 2n + 1

λ1nm = Dnm λ2nm =

(96) The formulas in (96) are straightforward to compute from Tables 1, 2 and 3 and the fact that the signature of the surface Laplacian s is −n(n + 1). It was conjectured in [11] that replacing S0 with Sik should improve conditioning. We investigate that possibility here, replacing (94) with − n · Hinc = n · H   Sik −∇s · (n × n × Einc ) = Sik (∇s · (n × n × E)) .

(97)

The signature in this case is the following: [n×∇×Sk+ ]

[Ek ] + Dnm = (Jn (k) + iJ n (k))Hn (k) 11 12  [n×∇×S + ]  [Sik ] [Ek ] k Dnm ikλnm − Dnm 22 21 2 ik jn (ik)h n (ik)(Jn (k) − iJ n (k))H n (k).

λ1nm = Dnm λ2nm = =

(98)

Figure 5 compares the condition number of the formulations G0 in (94) and Gik in (97). As above, the formulation using the Yukawa regularizing operator has a much better condition number in the high frequency regime. As shown in Fig. 6, the eigenvalues are also more clustered—improving the performance of iterative methods such as GMRES.

123

F. Vico et al. 10

1

Eigenvalues gDEBYE S0

Eigenvalues gDEBYE Sik 0.5 0

10

Imag

Imag

0

20

0.5 1

30 40

1.5 1

0

1

2

2

3

Real

1

0

1

2

3

Real

Fig. 6 Spectrum of the generalized Debye formulation for electromagnetic scattering at k = 50. The lefthand plot is for G0 and the right-hand plot is for Gik (the second equations in (94) and (97) are multiplied by a factor of 2 to make the cluster points coincident)

6.4 Transmission boundary value problems for the Helmholtz equation In many acoustic problems, an incoming wave is partially transmitted to the interior of the scatterer. In such settings, the Helmholtz coefficient for the interior and the exterior are typically different. We will denote by k0 the coefficient for the exterior region and k the coefficient for the interior region. The scattered field in the interior will be denoted by u and in the exterior by u 0 . The Helmholtz transmission boundary value problem [8] consists in finding the scattered field that satisfies the jump conditions − μ0 u + 0 − μu = f ∂u + 0 ∂n

−

∂u − ∂n

=g

(99)

for some given material parameters μ0 and μ and some functions f and g on the boundary. In [17] and [21], it was suggested that one use the representation u 0 (x) = Dk0 [σ ](x) + μ0 Sk0 [ρ](x) u(x) = Dk [σ ](x) + μSk [ρ](x),

(100) (101)

leading to the integral equation    μ0 Dk+0 − μDk− σ + μ20 Sk0 − μ2 Sk ρ = f     , − ρ = g Dk 0 − Dk σ + μ0 Sk + − μS k 0 where f is the jump in u and g is the jump in

∂u ∂n ,

that is:

− μ0 u + 0 − μu = f . ∂u 0 ∂u − ∂n − = g ∂n +

123

(102)

(103)

Boundary integral equation analysis on the sphere 200

800

150

DST

Condition Number

Condition Number

DST

DST SST

100

50

0

0

2

4

6

8

600

SST 400 200

0

10

k (wavenumber)

20

40

60

80

100

k (wavenumber)

Fig. 7 Condition numbers of single and double source formulations for the Helmholtz transmission problem. The plots on the left correspond to a contrast μμ = 1.42 > 1. The plots on the right correspond to a 0

1 <1 contrast μμ = 1.42 0

The signature of the operator, in matrix notation, is given by: ⎛ ] D[DST = nm

[Dk+ ] [Dk− ] 0 μ λ − μλ 0 nm nm ⎝ [Dk ] [D ] λnm 0 − λnmk

⎞

[Sk ] [Sk ] μ20 λnm0 − μ2 λnm ⎠. + [Sk ] [S − ] μ0 λnm0 − μλnmk

(104)

The superscript DST is used because this is a “double source” formulation, involving unknowns σ and ρ. While invertible for any k > 0, assuming μ, μ0 ∈ R+ , Eq. (102) behaves poorly with increasing frequency. Following a scaling suggestion in [15], we let ρ → k0 ρ and divide the second equation by k0 :    μ0 Dk+0 − μDk− σ + k0 μ20 Sk − μ2 Sk ρ = f     . + − 1 ρ = kg0 k0 Dk0 − Dk σ + μ0 Sk0 − μSk

(105)

The signature is now ⎛

[Dk+ ] 0 λ μ 0 nm ⎜

[DST ] Dnm =⎝

[D − ] − μλnmk

[Dk ] 1 0 k0 (λnm

[D ]

− λnmk )



[Sk ] [Sk ] k0 μ20 λnm0 − μ2 λnm [Sk + ] [S − ] μ0 λnm0 − μλnmk

⎞ ⎟ ⎠

(106)

The condition numbers for both equations are plotted in Fig. 7. When the right hand side of (103) comes from a scattering problem, f and g are in not independent since f = −μ0 u in and g = − ∂u∂n . Under these conditions, various “single source” formulations have been derived. We consider the formulation of [14] Dk+



   μ0  + Dk0 σ − ik0 Sk0 σ − Sk Dk 0 σ − ik0 Sk + σ = h, 0 μ

(107)

123

2

1

1

0

Eigenvalues SST.

Imag

Imag

F. Vico et al.

0

1

2

1

2 1

0

1

2

3

3

0

1

2

3

4

Real

Real

Fig. 8 Spectrum of the single source (SST ) and double source (DST ) integral equations for the Helmholtz 1 < 1 and k = 50 transmission problem with contrast μμ = 1.42 0

which follows from the representation: u 0 (x) = Dk0 [σ ](x) − ik0 Sk0 [σ ](x),

(108)

for u 0 (x) in the outer region and the representation formula: u(x) = Dk [u − ](x) − Sk

!

" ∂u (x), ∂n −

(109)

for u(x) in the inner region. The right hand side here depends on u inc , h=

Dk+



inc   μ0 inc ∂u , (u ) − Sk μ ∂n +

(110)

but we leave the derivation to the original paper. A Calderon identity shows that (107) is a Fredholm equation of the second kind. The signature of the single source formulation is given by: λ[SS] nm

=

[D + ] λnmk



  [Dk ] [Sk + ] μ0 [Dk+0 ] [Sk0 ] [Sk ] 0 0 . (111) (λnm − ik0 λnm ) − λnm λnm − ik0 λnm μ

Figure 7 shows a comparison between the three formulations—the double source integral equation (with and without the scaling factor) and the single source integral equation. Figure 8 shows that the spectrum {λ[SS] nm } is bounded away from zero for [DST ] contrast ratios smaller than one, while {λnm } has eigenvalues close to the origin. 7 Conclusions In this paper, we have analyzed and tabulated the behavior of a large library of layer potential operators when acting on scalar or vector density functions on the unit sphere.

123

Boundary integral equation analysis on the sphere

These tables permit the rapid analysis of the condition number and spectral properties of integral equations used in acoustic or electromagnetic scattering. The analytic expressions allow for the optimization of tunable parameters that appear in many formulations—both in regularizing operators (as investigated above) and in the selection of coupling constants [such as the choice of η in Eq. (89) or (90)]. We have concentrated here on acoustic and electromagnetic problems, but the same formulae can be used for analyzing integral equations in fluid dynamics or elasticity. Acknowledgments This work was supported in part by the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under NSSEFF Program Award FA9550-10-1-0180 and by the Department of Energy under contract DEFGO288ER25053. This work was supported also by the Spanish Ministry of Science and Innovation (Ministerio de Ciencia e Innovacion) under the projects CSD200800068 and TEC2010-20841-C04-01.

References 1. Abramowitz, M., Stegun, A.: Handbook of Mathematical Functions. Dover, New York (1964) 2. Borel, S., Levadoux, D., Alouges, F.: A new well-conditioned integral formulaiton for Maxwell equations in three dimensions. IEEE Trans. Antennas Propag. 9, 2995–3004 (2005) 3. Bruno, O., Elling, T., Paffenroth, R., Turc, C.: Electromagnetic integral equations requiring small numbers of Krylov-subspace iterations. J. Comput. Phys. 228, 6169–6183 (2009) 4. Boubendir, Y., Turc, C.: Wave-number estimates for regularized combined field boundary integral operators in acoustic scattering problems with Neumann boundary conditions. IMA J. Numer. Anal. (2013). doi:10.1093/imanum/drs038 (published online: March 7) 5. Bruno, O., Elling, T., Turc, C.: Well-conditioned high-order algorithms for the solution of threedimensional surface acoustic scattering problems with Neumann boundary conditions (preprint) 6. Chandler-Wilde, S.N., Graham, I.G., Langdon, S., Lindner, M.: Condition number estimates for combined potential boundary integral operators in acoustic scattering. J. Integr. Equ. Appl. 21, 229–279 (2009) 7. Chandler-Wilde, S.N., Graham, I.G., Langdon, S., Spence, E.A.: Numerical-asymptotic boundary integral methods in high-frequency acoustic scattering. Acta Numer. 21, 89–305 (2012) 8. Colton, D., Kress, R.: Integral Equation Methods in Scattering Theory. Wiley, New York (1983) 9. Colton, D., Kress, R.: Inverse Acoustic and Electromagnetic Scattering Theory. Springer, Berlin (1992) 10. Contopanagos, H., Dembart, B., Epton, M., Ottusch, J., Rokhlin, V., Visher, J., Wandzura, S.: Wellconditioned boundary inte- gral equations for three-dimensional electromagnetic scattering. IEEE Trans. Antennas Propag. 50, 1824–1830 (2002) 11. Epstein, C.L., Greengard, L.: Debye sources and the numerical solution of the time harmonic Maxwell equations. Commun. Pure Appl. Math. 63, 0413–0463 (2010) 12. Hsiao, G., Kleinman, E.: Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics. IEEE Trans. Antennas Propag. 45, 316–328 (1997) 13. Jackson, J.D.: Classical Electrodynamics. Wiley, New York (1975) 14. Kleinman, R., Martin, P.: On single integral equations for the transmission problem of acoustics. SIAM J. Appl. Math. 48, 307–325 (1988) 15. Kress, R.: Minimizing the condition number of boundary integral operators in acoustic and electromagnetic scattering. SIAM J. Appl. Math. 48, 307–325 (1988) 16. Kress, R.: Linear Integral Equations. Springer, Heidelberg (1999) 17. Kress, R., Roach, G.: Transmission problems for the Helmholtz equation. J. Math. Phys. 19, 1433–1437 (1978) 18. Nédélec, J.-C.: Acoustic and Electromagnetic Equations. Springer, New York (2001) 19. Panich, I.: On the question of the solvability of the exterior boundary problem for the wave equation and Maxwell’s equations. Uspekhi Mat. Nauk. 20, 221–226 (1965) 20. Papas, C.H.: Theory of Electromagnetic Wave Propagation. Dover, New York (1988) 21. Rokhlin, V.: Solution of acoustic scattering problems by means of second kind integral equations. Wave Motion 5, 257–272 (1983)

123