Regularized combined field integral equations for acoustic transmission problems

# Regularized combined field integral equations for acoustic transmission problems

Yassine Boubendir, Víctor Domínguez, David Levadoux, Catalin Turc
New Jersey Institute of Technology, Universidad Pública de Navarra, Spain,
ONERA France, New Jersey Institute of Technology
###### Abstract

We present a new class of well conditioned integral equations for the solution of two and three dimensional scattering problems by homogeneous penetrable scatterers. Our novel boundary integral equations result from suitable representations of the fields inside and outside the scatterer as combinations of single and double layer potentials acting on suitably defined regularizing operators. The regularizing operators are constructed to be suitable approximations of the admittance operators that map the transmission boundary conditions to the exterior and respectively interior Cauchy data on the interface between the media. The latter operators can be expressed in terms of Dirichlet-to-Neumann operators. We refer to these regularized boundary integral equations as Generalized Combined Source Integral Equations (GCSIE). The ensuing GCSIE are shown to be integral equations of the second kind in the case when the interface of material discontinuity is a smooth curve in two dimensions and a smooth surface in three dimensions.
Keywords: transmission problems, Combined Field Integral Equations, regularizing operators, Dirichlet-to-Neumann operators.
MSC2010: 35J05, 47G10, 45A05

## 1 Introduction

Numerical methods based on integral equation formulations for the solution of scattering problems, when applicable, have certain advantages over those that use volumetric formulations, largely owing to the dimensional reduction and the explicit enforcement of the radiation conditions. A crucial requirement of reformulating a linear, constant-coefficient PDE in terms of boundary integral equation formulations is that the latter are well-posed. In the case when the boundary of the scatterer is regular enough, there is a myriad of possibilities to derive well-posed boundary integral equations. Amongst those, the most widely used methodology of deriving well-posed integral equations for solution of scattering problems relies on Combined Field Integral Equations (CFIE) [6, 8, 18, 13]. Here in what follows we assume that the boundary of the scatterer is regular enough; the case of less regular interfaces (e.g. Lipschitz) is much less understood. In the scalar case, this methodology seeks scattered fields in terms of suitable linear combinations of single and double layer potentials so that the enforcement of the boundary conditions leads to boundary integral equations (CFIE) whose underlying operators are Fredholm in suitable boundary trace spaces of scattering problems. The well-posedness of CFIE is then settled via uniqueness arguments for the Helmholtz equation with certain boundary conditions (e.g. impedance boundary conditions, transmission boundary conditions).

Solvers based on integral equation formulations of scattering problems lead to systems of linear equations that involve dense matrices that can be very large in the high-frequency regime. Due to the large size of the underlying matrices, the solution of these linear algebra problems relies on Krylov subspace iterative solvers and is greatly facilitated by the availability of fast algorithms to perform matrix vector products. Thus, the efficiency of scattering solvers based on integral equations hinges a great deal on the spectral properties of the integral operators that enter the integral formulations, spectral properties themselves that influence the speed of convergence of iterative solvers. Thus, boundary integral equations of the second kind (e.g. equations whose operators are compact perturbations of identity in appropriate functional spaces) are preferable for the solution of scattering problems. The combined field strategy delivers boundary integral equations of the second kind for scalar scattering problems with Dirichlet boundary conditions [6] and transmission boundary conditions [18, 13, 11, 19]. However, it is typically the case that CFIE formulations, although extremely reliable, do not necessarily possess the best spectral properties amongst all well-posed formulations possible. For instance, in the case of scattering problems with other boundary conditions such as Neumann and impedance boundary conditions, various regularization procedures can deliver integral equations of the second kind [6, 8, 1, 2, 14, 7, 5] with superior spectral properties. The regularization procedure is a means to derive systematically well-conditioned boundary integral formulations of scattering problems. The techniques are quite general and can be applied, in principle, to any boundary value problem for linear, constant-coefficient PDEs or systems of PDEs.

All the regularization procedures for scalar scattering problems rely on the same main steps: (1) the fields (i.e. the solutions of the Helmholtz equations) are represented in each domain of interest with the aid of Green’s formulas in terms of the Cauchy data (i.e. Dirichlet and Neumann traces) on the boundary of each domain; (2) an abstract operator that maps the given boundary conditions to all the Cauchy data needed (e.g: in the case of Neumann boundary conditions, the operator maps Neumann traces to Dirichlet traces and is thus a Neumann-to-Dirichlet operator)—the operator is defined in terms of Dirichlet-to-Neumann operators and possibly their inverses; (3) the fields are represented using Green’s formulas where the Cauchy data is represented as an operator that approximates the operator acting on some unknown sources/densities. The enforcement of the boundary conditions on this representation leads to Generalized Combined Source Integral Equations (GCSIE) / Regularized Combined Field Integral Equations (CFIER). We note at this stage that by construction, if the operator were used instead of in the GCSIE, then these would consist of the identity operator. In the next step (4) the degree of approximation of the operator (that is the degree of smoothing of the difference operator ) is established so that second kind Fredholm GCSIE are obtained in appropriate boundary trace spaces of scalar scattering problems. The desired degree of smoothing is achieved provided that (5) the operator is constructed via suitable approximations of Dirichlet-to-Neumann operators. In order to meet the additional requirement that the ensuing GCSIE operators are injective (and thus invertible with continuous inverses), the aforementioned approximations of Dirichlet-to-Neumann operators are constructed through complexification of boundary integral operators [7, 5], through complexification of the wavenumbers in the definition of Dirichlet-to-Neumann operators for simple geometries such as the half-planes/half-spaces corresponding to tangent lines/planes to the boundary of the scatterer [1, 2]—in these cases the Dirichlet-to-Neumann maps can be defined by Fourier multipliers, or through boundary integral operators corresponding to the wavenumbers of the Helmholtz equations under consideration and quadratic partitions of unities [14, 15, 16]. Calderón’s identities are a crucial ingredient in the calculus in part (5).

We present in this paper novel integral equation formulations of two and three dimensional scalar transmission scattering problems that apply the five-step program outlined above. These integral equation formulations are actually systems of integral equations whose unknowns are certain densities defined on the interface of material discontinuity. There are two possibilities in terms of the functional spaces in which we seek those densities and in which we aim to construct GCSIE operators that are Fredholm of the second kind: (i) we assume that both densities belong to the same boundary Sobolev space and (ii) we assume that one of the densities has one more order of regularity than the other. We note that numerical methods based on GCSIE with property (i) are more accurate, and more amenable to an error analysis. We note that this distinction plays an important role in the construction of appropriate approximations of the operator in step (2) which maps the difference of exterior and interior Dirichlet and Neumann traces on the interface of material discontinuity to the Cauchy data of transmission problems. The matrix operator is expressed in terms of compositions of exterior and interior Dirichlet-to-Neumann operators corresponding to two different domains and wavenumbers and inverses of operators that involve linear combinations of those. Depending on the case (i) and (ii), the degree of smoothing we require on the difference operator is different. The approximating operators , in turn, are constructed per step (5) above via approximations of both exterior and interior Dirichlet-to-Neumann operators. The degree of these latter approximations is different according to case (i) and case (ii) respectively, and the dimension of the ambient space in which we solve the transmission problem. We rely on Calderón’s calculus to construct approximations of the Dirichlet-to-Neumann operators in terms of normal derivative of double layer operators corresponding to complex wavenumbers in case (i) and two dimensional ambient space and in case (ii) and three dimensional ambient space; or linear combinations of those with compositions of the normal derivative of double layer operators and double layer potentials corresponding to imaginary wavenumbers in case (ii) and three dimensional ambient space. The positivity of the imaginary parts of the former operators allows us to establish the injectivity and thus the invertibility of the GCSIE operators in both case (i) and case (ii). As it was illustrated in [3], solvers based on the GCSIE formulations, on account of the superior spectral properties of these formulations, outperform solvers based on classical integral formulations of transmission problems [13, 11, 19].

The paper is organized as follows: in Section 2 we present the acoustic transmission problems and review the basic properties of scattering boundary integral operators; in Section 3 we define and compute the admittance operator of transmission problems and we set up regularized integral equations in the form of Generalized Combined Source Integral Equations (GCSIE) that are based on regularizing operators that approximate the operators ; in Section 4 we derive sufficient conditions on the regularizing operators so that they lead to Fredholm second kind GCSIE; in Section 5 we construct approximations of Dirichlet-to-Neumann operators for each medium that lead via Caldéron’s calculus to constructions of regularizing operators ; in Sections 67, and 8 we establish the Fredholm properties of the GCSIE for various choices of regularizing operators in two and three dimensions; finally, in Section 9 we establish the well-posedness of the GCSIE.

## 2 Integral Equations acoustic transmission

### 2.1 Acoustic transmission problem

We consider the problem of evaluating the time-harmonic fields and that result as an incident field impinges upon the boundary of a homogeneous penetrable scatterer which occupies a bounded region in . The frequency domain acoustic transmission problem is formulated in terms of finding fields and that are solutions to the Helmholtz equations

 Δu2+k22u2=0in D2, (1)
 Δu1+k21u1=0in D1=Rd∖D2, (2)

given an incident field that satisfies

 Δuinc+k21uinc=0in D1, (3)

where the wavenumbers are the wavenumbers corresponding to the domains respectively. In addition, the fields , , and are related on the boundary by the the following boundary conditions

 γ1Du1+γ1Duinc = γ2Du2on Γ γ1Nu1+γ1Nuinc = νγ2Nu2on Γ. (4)

In equations (2.1) and what follows denote exterior and respectively interior Dirichlet traces, whereas denote exterior and respectively interior Neumann traces taken with respect to the exterior unit normal on . We assume in what follows that the wavenumbers are positive and that the density ratio is also positive. We note that in the case , equations (1)-(3) can also model electromagnetic scattering by two-dimensional penetrable obstacles , in which case (TE case) or (TM case). We assume in what follows that the boundary is a closed and smooth curve in and a closed and smooth surface in . We furthermore require that satisfies Sommerfeld radiation conditions at infinity:

 lim|r|→∞r(d−1)/2(∂u1/∂r−ik1u1)=0. (5)

Under the assumption that , , and are real and positive, it is well known that the systems of partial differential equations (1)-(3) together with the boundary conditions (2.1) and the radiation condition (5) has a unique solution [13, 12]. The results in this text can be extended to the case of complex wavenumbers , provided we assume uniqueness of the transmission problem and its adjoint.

### 2.2 Layer integral potentials and operators

A variety of integral equations for the transmission problem (1)-(2.1) exist [13, 10, 12]. The starting point in the derivation of direct integral equations for transmission problems is the Green’s identities. Hence, let

 Gk(x):=i/4(k|x|−1/2π)(d−2)/2H(1)(d−2)/2(k|x|)

the free space Green’s functions corresponding to the Helmholtz equation with wavenumber . For the sake of a simpler exposition, from now on we will commit a slight abuse of notation and denote

 Gj(x)=Gkj(x),j=1,2.

(The context will avoid any possible confusion).

Next we define the associated single and double layer potential

 [SLkφ](z):=∫ΓGk(z−y)φ(y)dσ(y),[DLkψ](z):=∫Γ∂Gk(z−y)∂n(y)ψ(y)dσ(y)

for . As before, denotes the layer potentials for the wavenumbers , with .

We have then the representation formulas for the exterior and interior domain

 u1=DL1(γ1Du1)−SL1(γ1Nu1),u2=−DL2(γ2Du2)+SL2(γ2Nu2). (6)

In addition to Green’s identities (6), trace formulas of the single and double layer potential are needed in the derivation of integral equations for transmission problems. For a given wavenumber , the traces on of the single and double layer potentials corresponding to the wavenumber and densities and are given by

 γ1DSLk(φ)=γ2DSLk(φ)=SkφγjNSLk(φ)=(−1)j2φ+K⊤kφj=1,2γjDDLk(ψ)=(−1)j+12ψ+Kkψj=1,2γ1NDLk(ψ)=γ2NDLk(ψ)=Nkψ. (7)

In equations (7) the operators and are the double layer and the adjoint of the double layer operator defined for a given wavenumber and density as

 (Kkφ)(x) := ∫Γ∂Gk(x−y)∂n(y)φ(y)dσ(y), x on Γ, (K⊤kφ)(x) := ∫Γ∂Gk(x−y)∂n(x)φ(y)dσ(y), x on Γ.

Furthermore,

 (Nkφ)(x):=FP∫Γ∂2Gk(x−y)∂n(x)∂n(y)φ(y)dσ(y)

is the so-called hypersingular operator (FP stands for the Hadamard Finite Part Integral). We point out

 (Nkφ)(x)=k2∫ΓGk(x−y)(n(x)⋅n(y))φ(y)dσ(y)+PV∫Γ∂s(x)Gk(x−y)∂s(y)φ(y)dσ(y),

when , and

 (Nkφ)(x)=k2∫ΓGk(x−y)(n(x)⋅n(y))φ(y)dσ(y)+PV∫Γ−−→curlxΓ Gk(x−y)⋅−−→curlyΓ φ(y)dσ(y),

when , where PV denotes the Cauchy Principal value of the integral, the tangential derivative and . Finally, the single layer operator is defined as

 (Skφ)(x):=∫ΓGk(x−y)φ(y)dσ(y), x on Γ.

Again, we will use and for for denoting the layer operator associated to the wavenumbers .

Having recalled the definition of the scattering boundary integral operators, we present next their mapping properties in appropriate Sobolev spaces of functions defined on the manifold  [17, 3]:

###### Theorem 2.1

For a smooth curve/surface the mappings

are continuous for all .

We will use throughout the text the following results about the smoothing properties of differences of boundary integral operators corresponding to different wavenumbers

###### Theorem 2.2

Let and be such that and . Then,

 Sκ1−Sκ2:Hs(Γ)→Hs+3(Γ),Nκ1−Nκ2:Hs(Γ)→Hs+1(Γ) (8)

are continuous for any .

Moreover, for ,

 Kκ1−Kκ2,K⊤κ1−K⊤κ2:Hs(Γ)→Hs+2(Γ) (9)

are also continuous for all .

Proof. Let and set . We assume to be sufficiently large so that . Fix and denote

 v1:=SLκ1φ,v2:=SLκ2φ.

Observe that for cf [17, Cor 6.14]

 ∥vj∥Hs+3/2(DR)+∥vj∥Hs+3/2(D2)≤C∥φ∥Hs(Γ),

with independent of . Then satisfies

 Δω=f,in DR∪D2,γ1Nω−γ2Nω=0,γ1Dω−γ2Dω=0,

with . Observe that, due to (7), has a jump in the normal derivative across , no matter how smooth is. Therefore, has a limited regularity in . However, we can apply [17, Theorem 4.20] to prove that and that there exists , again independent of , so that

 ∥ω∥Hs+7/2(DR)+∥ω∥Hs+7/2(D2)≤C[∥f∥Hs+3/2(DR)+∥f∥Hs+3/2(D2)]≤C′∥φ∥Hs(Γ).

In other words,

 SLκ1−SLκ2:Hs(Γ)→Hs+7/2(D2),SLκ1−SLκ2:Hs(Γ)→Hs+7/2(DR)

are continuous for all and any sufficiently large .

Analogously, one can show and are continuous for all and any sufficiently large .

By (7), and the continuity of the trace and normal derivative in the appropriate Sobolev spaces, we can easily conclude that

 Sκ1−Sκ2 :Hs(Γ)→Hs+3(Γ), Nκ1−Nκ2 :Hs+1(Γ)→Hs+2(Γ) K⊤κ1−K⊤κ2 :Hs(Γ)→Hs+2(Γ), Kκ1−Kκ2 :Hs+1(Γ)→Hs+3(Γ), (10)

are continuous for any . A transposition argument in proves that, again for ,

 Sκ1−Sκ2 :H−s−3(Γ)→H−s(Γ), Nκ1−Nκ2 :H−s−2(Γ)→H−s−1(Γ). Kκ1−Kκ2 :H−s−2(Γ)→H−s(Γ), K⊤κ1−K⊤κ2 :H−s−3(Γ)→H−s−1(Γ), (11)

which with (10), proves (8) for almost all . The remaining values are covered using the theory of interpolation in Sobolev spaces cf [17, App. B]. Hence, for instance, we have already proved that is continuous for . For , we have, in the usual notation of interpolation spaces, and, since , we can make use of (10)-(2.2) and the theory of interpolation spaces to conclude

 ∥Sκ1−Sκ2∥Hs(Γ)→Hs+3(Γ)≤∥Sκ1−Sκ2∥1/2Hs+2(Γ)→Hs+5(Γ)∥Sκ1−Sκ2∥1/2Hs−2(Γ)→Hs+1(Γ).

The continuity for the excluded of the remaining boundary layer operators are dealt in a similar way.

## 3 Regularized combined source integral equations for transmission problems

We derive in this section regularized combined field integral equations for transmission problems that rely on the use of adequate approximations of the Dirichlet to Neumann (DtN) operators. We start by explaining the main idea of our strategy. To this end, we introduce several notations. For a field such that and are solutions of the equations (2) and (1) respectively, where in addition is radiative, we define by the operator that maps the field to the boundary data of the transmission problems given in equations (2.1), that is

 γT(u)=(γ1Du1−γ2Du2γ1Nu1−νγ2Nu2).

At the heart of our approach there are two operators and that for a given defined as above map the given transmission boundary conditions (e.g. the operator ) to the Cauchy data on given by the exterior and interior Dirichlet and Neumann traces of . More specifically, we denote by the matrix operator which maps the difference of exterior and interior Dirichlet and Neumann traces of the field to the exterior Dirichlet and Neumann traces on of the component on . We write this as

 R1γT(u)=γ1C(u),% whereγ1C(u)=(γ1Du1γ1Nu1). (12)

Similarly, we denote by the operator that maps the difference of exterior and interior Dirichlet and Neumann traces of the field to the interior traces on of the component , that is

 R2γT(u):=γ2C(u)% whereγ2C(u)=(γ2Du2γ2Nu2). (13)

On account of the boundary conditions in equations (2.1), it follows that , where denotes the identity matrix. The field itself can be retrieved through the Green’s formulas, i.e. equations (6) from the Cauchy data on . We write this in operator form as

 u=(DL1−SL100)(γ1Du1γ1Nu1)+(00−DL2SL2)(γ2Du2γ2Nu2)

or in short form as . Obviously, if we denote by , the following identity holds

 γTCR=I. (14)

The operators can be expressed in terms of Dirichlet to Neumann operators. However, their evaluation is numerically cumbersome, if at all possible. Our idea is to use instead certain suitable approximations of the operator defined in equation (21). Once operators are constructed, their counterparts can be taken to be . We look then for a field in the form , where is a vector density defined on and . More precisely, if we denote , we look for fields defined as

 u1(z)=DL1(˜R11a+˜R12b)(z)−SL1(˜R21a+˜R22b)(z),z∈Rd∖Γ (15)

and defined as

 u2(z)=−DL2(˜R11a+˜R12b−a)(z)+ν−1SL2(˜R21a+˜R22b−b)(z),z∈Rd∖Γ. (16)

Using the jump conditions of the boundary layer potentials we are led to the integral equation

 ˜D(ab)=γTC(˜Rw)=−(γ1Duincγ1Nuinc) (17)

which we refer to as Generalized Combined Source Integral Equation (GCSIE) and which takes on the following explicit form:

 (12I−K2+(K1+K2)˜R11−(S1+ν−1S2)˜R21)a +(ν−1S2+(K1+K2)˜R12−(S1+ν−1S2)˜R22)b = −γ1Duinc (−νN2+(N1+νN2)˜R11−(K⊤1+K⊤2)˜R21)a +(12I+K⊤2+(N1+νN2)˜R12−(K⊤1+K⊤2)˜R22)b = −γ1Nuinc. (18)

In what follows we compute the operator in terms of Dirichlet to Neumann operators and we establish in what sense should the operators approximate the operator so that, in the light of formulas (14) and (17), the matrix operators in the left-hand side of GCSIE equations (3) are close to the identity matrix. Notice that is precisely what we obtain when using in (3) instead.

We present next a formal calculation of the operator based on the use of Dirichlet-to-Neumann operators for the domains . The latter operators are defined such that maps the Dirichlet trace on the boundary of a radiative solution of the Helmholtz equation with wavenumber in the domain to its Neumann trace on the boundary , i.e. and maps the Dirichlet trace on the boundary of a solution of the Helmholtz equation with wavenumber in the domain to its Neumann trace on the boundary , i.e. . By (12),

 R11(γ1Du1−γ2Du2)+R12(γ1Nu1−νγ2Nu2)=γ1Du1. (19)

Equation (19) can be further expressed in term of the admittance operators as the following (operator) linear system

 R11+R12Y1 = I R11+νR12Y2 = 0. (20)

We immediately obtain that the solution of the linear system in equation (3) is given by and . Using the admittance operator we see that the second row of the matrix operator can be obtained by composing on the left the first row of by . Hence we obtain

 R1=(−ν(Y1−νY2)−1Y2(Y1−νY2)−1−νY1(Y1−νY2)−1Y2Y1(Y1−νY2)−1). (21)

We note that the calculations that led to equation (21) are entirely formal, as the operators or may not be well defined. Furthermore, the computation of the Dirichlet to Neumann operators is expensive for general domains . Nevertheless, operators can be constructed such that the difference operators are smoother operators than (these type of operators are referred to as parametrices). We discuss in what follows what degree of smoothing must the operators have in order to lead to GCSIE operators (17) that are Fredholm of the second kind in appropriate Sobolev spaces.

## 4 Approximations of the admittance operators

Our goal is to produce appropriate approximations of the exact admittance operator so that the matrix operators that enter GCSIE formulations (17) are (i) compact perturbations of the identity (matrix) operator in appropriate Sobolev spaces and (ii) invertible in the same spaces. We establish in this section sufficient conditions on the regularity properties of the difference matrix operators that ensure the aforementioned property (i). We distinguish two cases with regards to Sobolev spaces the matrix operators act upon:

Case I we consider and which implies that the solution of the GCSIE formulations (17) has the same regularity, that is ;

Case II we consider and which implies .

Given that , we expect that once we construct operators with the desired properties (i) and (ii), the eigenvalues of the operators in the left-hand side of equation (17) will accumulate at . In addition, we strive to construct operators that are (iii) as simple as possible so that their evaluation is as numerically inexpensive as possible.

We present first a result that establishes in what sense should approximate in order for the first property (i) to hold in Case I. In order to make a more striking distinction between the two cases, we denote by the approximating operators of the operators in the spaces . Given the mapping properties of the Dirichlet to Neumann operators , and assuming that the operators are well defined, we have then the following mapping properties of the components of the matrix operator defined in equation (21): , , , and . Our first important result is given in

###### Theorem 4.1

Assume that the operators and are well defined. Let be operators such that for all we have

Then the matrix operator defined in equation (17) and corresponding to the regularizing operator has the following mapping property . Furthermore, the operator can be written in the form

 ˜Ds,s=I+(˜Dr11˜Dr12˜Dr21˜Dr22)

where , , , and . In particular, and given the compact embeddings of into for all , it follows that the matrix operator is a compact perturbation of identity in the space .

Proof. Given the assumptions about the operators , it follows that they have the same mapping properties as the operators for , and thus the mapping property of the operators follow immediately. Let us consider first the identity (see (14))

 Ds,s−I=γTCRs,s−γTCR.

The result of the Theorem follows once we show that are regularizing operators of one order for all and . A simple calculation gives

 ˜Dr11=˜Ds,s11−I=˜Ds,s11−D11=(K1+K2)(˜Rs,s11−R11)−(S1+ν−1S2)(˜Rs,s21−R21).

We have that and , together with and , and thus . We also have

 ˜Dr12=˜Ds,s12=˜Ds,s12−D12=(K1+K2)(˜Rs,s12−R12)−(S1+ν−1S2)(˜Rs,s22−R22)

and hence since , , , and . Furthermore, given that

 ˜Dr21=˜Ds,s21=˜Ds,s21−D21=(N1+νN2)(˜Rs,s11−R11)−(K⊤1+K⊤2)(˜Rs,s21−R21)

we obtain that since , , , and . Finally, we have that

 ˜Dr22=˜Ds,s22−I=˜Ds,s22