Second-kind integral solvers for TE and TM problems of diffraction by open arcs

Second-kind integral solvers for TE and TM problems of diffraction by open arcs

Oscar P. Bruno and Stéphane K. Lintner

We present a novel approach for the numerical solution of problems of diffraction by open arcs in two dimensional space. Our methodology relies on composition of weighted versions of the classical integral operators associated with the Dirichlet and Neumann problems (TE and TM polarizations, respectively) together with a generalization to the open-arc case of the well known closed-surface Calderón formulae. When used in conjunction with spectrally accurate discretization rules and Krylov-subspace linear algebra solvers such as GMRES, the new second-kind TE and TM formulations for open arcs produce results of high accuracy in small numbers of iterations—for low and high frequencies alike.

1 Introduction

Problems of diffraction by infinitely thin open surfaces play central roles in a wide range of problems in science and engineering, with important applications to antenna and radar design, electronics, optics, etc. Like other wave scattering problems, open surface problems can be treated by means of numerical methods that rely on approximation of Maxwell’s Equations over volumetric domains (on the basis of, e.g., finite-difference or finite-element methods) as well as methods based on boundary integral equations. As a result of the singular character of the electromagnetic fields in the vicinity of open edges, open-surface scattering configurations present major difficulties for both volumetric and boundary integral methods. Boundary integral approaches, which require discretization of domains of lower dimensionality than those involved in volumetric methods, can generally be treated efficiently, even for high-frequencies, by means of accelerated iterative scattering solvers [3, 6, 24]. Unfortunately, integral methods for open surfaces do not give rise, at least in their classical formulations, to Fredholm integral operators of the second-kind, and can therefore prove to be computationally expensive—as the eigenvalues of the resulting equations accumulate at zero and/or infinity and, thus, iterative solution of these equations often requires large numbers of iterations.

This paper presents new Fredholm integral equations of second-kind and associated numerical algorithms for problems of scattering by two-dimensional open arcs (i.e., infinite cylinders of cross-section ) under either transverse electric (TE) or transverse magnetic (TM) incident fields. The new second-kind Fredholm equations (which result from composition of appropriately modified versions and of the classical single-layer and hypersingular integral operators and ) provide, for the first time, a generalization of the classical closed-surface Calderón formulas to the open-arc case. In particular, the new formulations possess highly favorable spectral properties: their eigenvalues are highly clustered, and they remain bounded away from zero and infinity, even for problems of very high frequency. When used in conjunction with spectrally accurate discretization rules and Krylov-subspace linear algebra solvers such as GMRES, the new open-arc formulations produce results of high accuracy in small numbers of iterations—for low and high frequencies alike.

The new second-kind formulation for the TM problem is particularly beneficial, as it gives rise to order-of-magnitude improvements in computing times over the corresponding weighted hypersingular formulation. Such gains do not occur in the TE case: although the new second-kind TE equation requires fewer iterations than the corresponding weighted first kind formulation, the total computational cost of the second-kind equation is generally higher in the TE case—since the application of the first-kind operator can be significantly less expensive than the application of the composite second-kind operator.

The difficulties that arise as integral equations are used to treat open surface scattering problems are of course well known, and many contributions have been devoted to their treatment. Like the present work, references [22] and [8] seek to tackle these problems by means of generalizations of the classical Calderón relations to the case of open surfaces. The first of these contributions establishes that the combination can be expressed in the form , where the kernel of the operator has local singularity of at most . This early result however does not take into account the singular edge behavior; the resulting operator is not compact (in fact it gives rise to extreme singularities at the edge [16]) and is therefore not a second-kind operator in any numerically meaningful functional space. When used in conjunction with boundary elements that vanish on the edges (by means of well-chosen projections) however, the combination does give rise to reduction of iteration numbers, as demonstrated in reference [8] through numerical examples for low frequency problems. This contribution does not include details on accuracy, and it does not utilize integral weights to resolve the solution’s edge singularity. A related but different method was introduced in [1] which exhibits, once again, low iteration numbers for low-frequency problems, but which does not resolve the singular edge behavior and for which no accuracy studies have been presented. Finally, high-order integration rules for the single-layer and hypersingular operators adapted to open arcs were introduced in a Galerkin framework in [29, 12, 28]. These methods have thus far only been applied for simple geometries and at low frequencies, and limited information is available on the actual convergence properties and performance of their computational implementations.

A second class of methods include those proposed in references [2, 20, 13]. The contribution [2], some aspects of which are incorporated in our method, treats the Dirichlet problem for Laplace’s equation by means of second kind equations. The basis of this approach lies in the observation that the cosine basis has the dual positive effect of diagonalizing the logarithmic potential for a straight arc and removing the singular edge behavior—so that the inverse of the logarithmic potential can be easily computed and used as a preconditioner to produce a second kind operator for a general arc. The approach [20], which also uses a cosine basis, treats the Neumann problem for the non-zero frequency Helmholtz equation with spectral accuracy by means of first kind equations. The contribution [13], finally, treats, just like [2], the Laplace problem by means of second-kind equations resulting from inversion of the straight arc logarithmic potential; like [20], further, it produces spectral accuracy through use of the cosine transforms. The second-kind integral approach developed in [2] and later revisited in [13] seems essentially limited to the specific problem for which it was proposed: neither an extension to the Neumann problem nor to the full three dimensional problem seem straightforward. And, more importantly, this approach does not lead to adequately preconditioned equations for non-zero frequencies: a simple experiment conducted in Section 3.4 shows that a direct generalization of the algorithm [2] to the Helmholtz problem generally requires significantly more linear algebra iterations than are necessary if the operator alone is used.

The remainder of this paper is organized as follows: after recalling in Section 2 the classical boundary integral formulations for the TE and TM open-arc scattering problems, in Sections 3.1 and 3.2 we present our new weighted operators and as well as certain periodized counterparts and (which are obtained by considering a sinusoidal changes of variables for source and observation points). Our main result, the generalized Calderón formula, is presented in Section 3.3. Section 3.4 then presents results of numerical evaluation of eigenvalues for a non-trivial open arc problem, illustrating the spectral properties of previous open-arc operators as well as the second kind operators introduced in this paper. Theoretical considerations concerning the new second-kind equations are presented in Section 4, including a succinct but complete proof of the open-arc Calderón formulae; a more detailed theoretical discussion, including full mathematical technicalities, can be found in [16]. The high-order quadrature rules we use for evaluation of the new integral operators are described in Section 5. Numerical results, finally, are presented in Section 6 for a wide range of frequencies and for various geometries (including a brief study of resonant open cavities) demonstrating the uniformly well conditioned character of the integral formulations proposed in this paper.

2 Preliminaries

Let denote a smooth open arc in the plane. The TE and TM problems of scattering by the open arc amount to two-dimensional problems for the Helmholtz equation with Dirichlet and Neumann boundary conditions on , respectively:


Here n denotes the normal to , and and are given in terms of the incident electric excitation : and on .

2.1 Classical open-arc equations for TE and TM problems

The (unique) solutions and of the TE and TM problems above can be expressed, for , in terms of single- and double-layer potentials of the form and ; see e.g. [29, 30, 26]. Here is a unit vector normal to at the point (we assume, as we may, that is a smooth function of ), and, calling the first Hankel function of order zero,


and Denoting by and the single-layer and hypersingular operators




further, the densities and are the unique solutions of the integral equations


It is known that the eigenvalues of the integral operators in Equations (5) accumulate at zero and infinity, respectively (see e.g. Figure 1). As a result (and as illustrated in Section 6), solution of these equations by means of Krylov-subspace iterative solvers such as GMRES generally require large numbers of iterations. In addition, as discussed in Section 2.3, the solutions and of equations (5) are not smooth at the end-points of , and, thus, they give rise to low order convergence (and require high discretization of the densities for a given accuracy) if standard quadrature methods are used. Details in these regards are provided in what follows.

2.2 Calderón relation for closed and open arcs.

As is well known [9, 21, 10] second-kind formulations for closed surfaces (arcs) result from: (1) The classical jump relations associated with singular integral operators, as well as (2) application of the Calderón formula


which shows that the composition of the closed-surface hypersingular and single-layer operators and equals a compact perturbation of the identity.

In the open-surface context however, jump relations cannot be exploited since, according to the boundary conditions (1), the same limits must be reached on both sides of the open surface . Further, use of the combined operator does not lead to well-posed equations: as shown in [16], the image under of the constant function is highly singular, with edge asymptotics , , where denotes the Euclidian distance to the edge. The combination is also problematic in the functional framework set forth in [29, 30, 27], since the image under of the Sobolev space (as defined in [27]) is larger than the domain of definition of (as defined in [27]) of the operator .

2.3 Regularity and singular behavior at the edge

The singular character of the solutions of equations (5) is well documented [27, 17]: and can be expressed in the forms


where denotes the distance to the edge, and are smooth cut-off functions, and where the functions and are somewhat smoother than and . More recently furthermore, it was established [17] that expressions such as (7) can be viewed as part of an expansion in powers of which can be carried out as long as the degree of smoothness of the surface and the right hand side of the equations allow it. For example, if the curve itself and the functions and in (5) are infinitely differentiable, it follows that


where and are infinitely differentiable functions throughout , up to and including the endpoints [17]. Thus the singular character of these solutions is fully characterized by the factors and in equation (8).

3 Well-posed Second-Kind Integral Equation Formulations

3.1 Weighted operators

In view of the regularity results (8), we introduce a positive integral weight with asymptotic edge behavior at the edge (by which it is implied that the quotient is infinitely differentiable up to the edge), and we define the weighted operators


so that, in view of the discussion of Section 2.3, for smooth functions and , the solutions of the equations


are also smooth.

Without loss of generality, we choose a smooth parameterization of defined in the interval , and we define the weight canonically as As a result, the operators and give rise to the parameter-space operators




defined on functions and of the variable , ; in these equations we have set . Clearly, for and we have and

3.2 Periodized spaces and regularity

Introducing the changes of variables and , and defining , we obtain the periodic weighted single-layer and hypersingular operators



Clearly then, the solutions to the periodic equations


where and , are related to the solutions of equations (10) by the relations


In view of (15), the solutions to equations (14) are smooth and periodic, and it is therefore natural to study the properties of and in the the Sobolev spaces of periodic and even functions cf. [31, 5].

Definition 1.

Let . The Sobolev space is defined as the set of -periodic functions of the form for which the -norm is finite.

Clearly the set is a basis of the Hilbert space for all .

Figure 1: Eigenvalue distribution for the spiral-shaped arc displayed in Figure 3, with frequency , for the various operators under consideration. Left . Left-center . Right-center . Right . Note the important difference of scale between the four plots.

3.3 Generalized Calderón Formula

Given the smoothness of the solutions of the equations arising from the weighted periodic operators and , and in light of equation (6), it is reasonable to consider the composition as possible basis for solution of open-arc problems. As stated in the following theorem and demonstrated below, this combination does indeed give rise to a second-kind integral formulation. The statement of the theorem makes use of the concept of bicontinuity: an operator between two Hilbert spaces is said to be bicontinuous if and only if is invertible and both and are continuous operators.

Theorem 1.
  • For all the operators and define bicontinuous mappings and . Thus, the composite operator is bicontinuous from into for all .

  • The generalized Calderón formula


    holds, where is a compact operator, and where is a bicontinuous operator independent of the frequency .

  • The set of eigenvalues of equals the union of the discrete set and a certain open set which is bounded away from zero and infinity. The sets are nested, they form a decreasing sequence, () and they satisfy , where denotes the closure of .

It follows that the open-arc TE and TM scattering problems (1) can be solved by means of the second-kind integral equations


respectively. The smooth and periodic solutions and of these equations are related to the singular solutions and of equations (5) via and We thus see that, through introduction of the weight and use of spaces of even and -periodic functions, a picture emerges for the open-surface case that resembles closely the one found for closed-surface configurations: the generalized Calderón relation (16) is analogous to the Calderón formula (6), and mapping properties in terms of the complete range of Sobolev spaces are recovered for and , in a close analogy to the framework presented in [14, 21] for the operators and .

3.4 Eigenvalue Distributions

In order to gain additional insights on the character of the various open-arc operators under consideration (namely, , , as well as a generalization to non-zero frequencies of the operator introduced in reference [2]), we consider their corresponding eigenvalue distributions. In Figure 1 we thus display the eigenvalues associated with these operators, for the spiral-shaped arc displayed in Figure 3 and described in Section 6, as they were produced by means of the quadrature rules presented in Section 5 and subsequent evaluation of matrix eigenvalues. The frequency was chosen to ensure a size to wavelength ratio , where is the length of the arc and the wavelength of the incident wave. As expected, the eigenvalues of tend slowly to zero, the eigenvalues of are large, while the eigenvalues of are bounded away from zero and infinity, and they accumulate at .

The generalization of the equation [2], whose operator eigenvalues are displayed in the center-right portion of Figure 1, is obtained from right-multiplication of the single layer operator in equation (10) by the inverse of the flat-arc zero-frequency single-layer operator (defined in equation (39) below); the resulting equation is given by


This equation, which can be re-expressed in the form , is a second-kind Fredholm integral equation: the operator is compact. Unfortunately, the spectrum of the operator in equation (19) is highly unfavorable at high frequencies, as illustrated in the center-right image in Figure 1. Such poor spectral distributions translate into dramatic increases, demonstrated in Table 6, in the number of iterations required to solve (19) by means of Krylov subspace solvers as the frequency grows. In fact, a direct comparison with Table 3 shows that the second-kind integral equation (19) may require many more iterations at non-zero frequencies than the original first-kind equation (14).

4 Theoretical Considerations

This section presents a proof of Theorem 1. The argument proceeds by consideration of the flat-arc zero-frequency case (Sections 4.1 through 4.3) and subsequent extension to the general case (Sections 4.4,4.5 and 4.6). A more detailed discussion, including full mathematical technicalities, can be found in [16].

4.1 Flat arc at zero-frequency

In the particular case where the frequency vanishes () and the curve under consideration is the flat strip (), the operator reduces to Symm’s operator


whose well-known diagonal property [31, 5] in the cosine basis ,


plays a central role in our analysis. Note that, in particular, equation (21) establishes the bicontinuity of the operator from into . The corresponding property of bicontinuity (from into ) for the weighted flat-strip, zero-frequency hypersingular operator

results, in turn, from our study of the operator in Section 4.2.

An integration by parts argument presented in reference [20] gives rise to the factorization where




As a result, for the flat-arc at zero-frequency case, the composite operator , is given by


4.2 The operator

It is easy to evaluate the action of the operator on the cosine basis: in view of equation (24), using (21) and expressing as a linear combination of cosines (which results from simple trigonometric manipulations) we obtain the relation


Clearly then,


where the operator is defined by


The right hand side of equation (26) may appear to bear a direct connection with the classical closed-curve Calderón formula (6). Yet, unlike the operator in equation (6), the operator is not compact–see Section 4.3. It is easy, however, to verify that can also be expressed in the form


and is therefore closely related to the Césaro operator Since is bounded (but not compact) from into  [4] (where denotes the space of square-integrable functions in the interval ), it follows that is a continuous operator for all , and, therefore, is also a continuous operator.

Composing the operator with from the right on one hand, and from the left, on the other, and applying the two resulting composite operators to the basis , shows that is the inverse of , and that is a continuous operator from into . It follows that is a bicontinuous operator. As indicated in Section 4.1, finally, the bicontinuity of from into follows directly from equation (24) and the corresponding bicontinuity properties of and .

4.3 Eigenvalues of

Re-expressing (25) in the form


and making use the well-known expansion [11, 19]


we obtain


where the diagonal elements are given by


In view of (29), the operator can be viewed as an infinite upper-triangular matrix. The diagonal terms of this matrix are eigenvalues of ; the corresponding eigenvectors , in turn, can be expressed in terms of a finite linear combination of the first cosine basis functions: . This establishes that the set defined in Theorem 1 is contained in the spectrum for all .

As is known, the diagonal elements of an infinite upper-triangular matrix are a subset of all the eigenvalues of the corresponding operator; for instance, the point spectrum of the upper-triangular bounded operator (the adjoint of the discrete Cesàro operator ) is given [4, Th. 2] by the entire disc . A similar situation arises for the operator : searching for sequences such that satisfies leads to the relations




The spectrum of equals the set of all values for which the relations (33) and (34) admit solutions satisfying . This set is the union of the discrete set (see point (iii) in Theorem 1), and the open set defined by


It follows in particular that the operator is not compact from into : if it were, in view of equation (26), the spectrum of would be discrete, in contradiction with (35)

We thus see that, even though is not compact (so that, in particular, the decomposition (26) does not present as the sum of a multiple of the identity and a compact operator), it follows from (32) and (35) that the eigenvalues of are all bounded away from zero and infinity, and that they are tightly clustered around —in close analogy to the eigenvalue distribution associated with the closed-surface Calderón operator (6).

TE case, TE case, TM case, TM case,
Table 1: Scattering by a spiral-shaped arc of size : far-field errors.

4.4 Properties of the operator

The study of in the general case (possibly curved arc, ) hinges on the decomposition


of the kernel , which results from the expression [10, p. 64]


where is the Bessel function and is analytic. Clearly, the functions and are smooth functions of and , and the singular behavior of thus resides entirely in a logarithmic term. In particular, the continuity of the operator follows directly from the corresponding continuity of the operator .

To establish the invertibility of the operator and the continuity of its inverse , on the other hand, we consider the decomposition


where the operator , defined by


differs from only in the additional multiplicative term Since the operator is compact (in view of classical embedding results and the increased regularity of the kernel of the operator ) and since the operator is injective (as it follows from results established in [29, 30, 27] for the operator (3)), the bicontinuity of follows from the invertibility of (and, thus, of ) together with an application of the Fredholm theory to the term in parenthesis in equation (38).

4.5 Properties of the operator

As is known, for closed curves, the operator (the normal derivative of the double-layer potential), can be expressed in terms of a composition of a single layer potential and derivatives tangential to the curve  [9, p. 57]. While for open surfaces the conversion of normal derivatives into tangential derivatives for the operator gives rise to highly singular boundary terms (which arise from the integration-by-parts calculation that is part of the conversion process), the presence of the boundary-vanishing weight in our operator eliminates all singular terms, and we obtain where


and where


see [20]. Using the changes of variables and in equations (40) and (41), it follows that




and where, defining ,


Note for future use that equation (44) can be re-expressed in the form



In view of equation (42), the operator equals the sum of (which, like , maps into ) and . But, from (45), we see that can be expressed as the sum of and a compact perturbation. Since, like , the operator maps into , it follows that is a bounded operator from into , and, thus, is a bounded operator as well.

The bicontinuity properties of can now be established by an argument similar to the one applied to in the previous section using, this time, the identity

Figure 2: Sample diffraction patterns. Left: TE scattering (total field) by an infinitely thin strip of size for horizontal left-to-right incidence. (The TM total field at this incidence equals the incident field.) Right: TM solution for a half-circle of size under normal incidence from below. Note the trailing shadow in the wake of the strip, and the caustics inside the circular reflector.

4.6 Generalized Calderón formula

The generalized Calderón formula is obtained by expressing the combination in the form


where the operator , given by is compact in view of the classical Sobolev embedding theorems and the (increased) smoothness properties of the kernels of the operators and respectively. In view of the relation and have the same spectrum: is an eigenvalue-eigenvector pair for if and only if is an eigenvalue-eigenvector pair for . Clearly, has the same mapping properties and regularity as ( is bicontinuous from into ), and therefore, the decomposition (46) shows that equals the sum of an invertible bicontinuous operator and a compact operator; that is, is a second-kind Fredholm operator.

5 High-Order Numerical Methods

In this section we present spectral quadrature rules for the operators , and which give rise to an efficient and accurate solver for the general open arc diffraction problem (1).

5.1 Spectral discretization for

Use of the nodes , , gives rise [23, eq. (5.8.7),(5.8.8)] to a spectrally convergent cosine representation for smooth, -periodic and even functions :


Thus, applying equation (21) to each term of expansion (47), we obtain the well-known spectral quadrature rule for the logarithmic kernel