A Coupling between Integral Equations and On-Surface Radiation Conditions for Diffraction Problems by Non Convex Scatterers

: The aim of this paper is to introduce an orignal coupling procedure between surface integral equation formulations and on-surface radiation condition (OSRC) methods for solving two-dimensional scattering problems for non convex structures. The key point is that the use of the OSRC introduces a sparse block in the surface operator representation of the wave ﬁeld while the integral part leads to an improved accuracy of the OSRC method in the non convex part of the scattering structure. The procedure is given for both the Dirichlet and Neumann scattering problems. Some numerical simulations show the improvement induced by the coupling method.


Introduction
During the last decades, time-harmonic wave propagation has proved to be central in many engineering and technological key developments, based, e.g., on acoustics, electromagnetism or elastic mechanisms. When one wants to simulate the associated boundary-value problem, one of the difficulties is related to the fact that the solution which has to be computed is set in an exterior unbounded domain Ω + defined as the complementary of a finite scatterer Ω − . Therefore, to use a standard numerical method, it is necessary to bound the computational domain. One well-known possibility is to use an absorbing boundary condition [1][2][3][4][5][6][7] or a Perfectly Matched Layer [8][9][10][11][12][13] to bound the domain and then solve the resulting problem by the finite element method [14][15][16]. Another widely used alternative is to rewrite equivalently the initial exterior PDE problem as an integral equation over the finite surface Γ of the scatterer Ω − based on the Green's function [17][18][19][20][21][22][23][24][25][26][27][28]. Then, this has the advantage of reducing from one the dimension of the problem. One of the main drawbacks is that, unlike the initial problem which involves partial differential operators, the integral equation is defined by construction as a nonlocal pseudodifferential operator. When a discretization technique is then applied, as the boundary element method, then the corresponding discrete version of the integral equation leads to the numerical solution of a highly indefinite complex-valued dense linear system which is particularly difficult to tackle in the high-frequency regime. Many technical aspects are then necessary to make it working correctly for some applications, for example to reduce the storage and to accelerate the solution of the linear system [17,18,29]. Various numerical methods were further developed to propose some other ways to solve, at least approximately, the initial scattering problem. One of them is the On-Surface Radiation Condition (OSRC) method introduced in [30] and further developed by many authors (see, e.g., [31]). Without giving too much details now, the OSRC approach also leads to solve an equation given over the surface of the scatterer, but defined through local surface partial differential operators. Therefore, after the application of the boundary element method, the linear system is highly sparse, yielding an efficient way to solve the scattering problem. The price to pay is that the method can be considered as a numerical asymptotic method, and therefore can lose some accuracy in some cases, in particular when the scatterer is not convex and includes some concave parts [31] where multiple scattering effects are present. Considering the problem of getting better accuracy for OSRC based techniques, a coupled formulation between the OSRC method and a two-dimensional finite element technique was proposed in [32] for scattering by single non convex obstacles. For the numerical simulation of multiple scattering problems (then leading to non convex global structures), coupled OSRC/integral equation formulations were presented in [33,34], where the shapes of the single obstacles are convex [33] or slightly non convex [34]. The aim of the present paper is to contribute to the improvement of the OSRC method for non convex single obstacles by proposing a direct simple coupling between the OSRC and the surface integral equation method. Unlike [32], the presented method is naturally only set on the boundary of the scatterer and does not suffer from the pollution error related to the finite element method.
The plan of the paper is the following: In Section 2, we introduce the two-dimensional Dirichlet scattering problem and the basic informations about the integral equation representations and their numerical approximation. Section 3 presents the notion of OSRC and its numerical discretization after writing the variational formulation. Section 4 develops the original coupling procedure for the Dirichlet problem. In particular, we explain how to formulate the problem and to improve its convergence properties if it is used in a Krylov solver, based on operator preconditioning. The coupling is validated in a simple two-dimensional example. The extension to the Neumann problem is shortly given in Section 5. Finally, Section 6 is a conclusion.

The Two-Dimensional Scattering Problem
Let us consider Ω − as a scatterer with polygonal boundary Γ := Ω − . The homogeneous isotropic exterior domain of propagation is denoted by Ω + = R 2 \ Ω − . For the sake of conciseness in the presentation, we first assume that the scatterer is acoustically sound-soft (i.e., Dirichlet boundary condition). Nevertheless, the case of a sound-hard scattering problem (Neumann boundary condition) is also treated shortly in Section 5. We now consider a time-harmonic incident plane wave u inc (x) = e ikθ inc ·x (with x = (x 1 , x 2 ) ∈ R 2 ) illuminating Ω − , with an incidence direction θ inc = (cos(θ inc ), sin(θ inc )) for a time dependence e −iωt , setting ω as the wave pulsation and k as the wavenumber. The sound-soft scattering problem of u inc by Ω − leads to the computation of the scattered wavefield u ∈ C 2 (Ω + ) ∩ C 0 (Ω + ) as the solution to the boundary-value problem We designate by (∆ + k 2 ) the Helmholtz operator, where ∆ = ∂ 2 x 1 + ∂ 2 x 2 is the laplacian. The gradient operator is ∇ and ||x|| = √ x · x, where x · y is the scalar product of two vectors x and y of R 2 . The last equation of (1) is the well-known Sommerfeld's radiation condition at infinity that ensures the uniqueness of the scattered wave field u. The outwardly directed unit normal vector to Ω − is n. A schematic representation of the problem is given in Figure 1. The existence and uniqueness of the solution to this BVP is well-known and detailed in [20], Theorem 2.12. From standard arguments connecting formulations in classical and Sobolev spaces ( [20], p. 107), then u ∈ H 1 loc (Ω + ) and u is C ∞ up to Γ, excluding the corners of the polygonal curves ( [20], Lemma 2.35). Following [20], Theorem 2.12, then it can also be proved that ∂u ∂n ∈ L 2 (Γ).

Integral Operators for Scattering
Let us define G as the two-dimensional free-space Green's kernel where H

Proposition 1.
If v ∈ C 2 (Ω + ) ∩ H 1 loc (Ω + ) is a solution to the Helmholtz equation in an unbounded domain Ω + which also satisfies the Sommerfeld radiation condition, then the following relation holds In addition, let us assume that v − ∈ C 2 (Ω − ) ∩ H 1 (Ω − ) is solution to the Helmholtz equation in a bounded domain Ω − . One can write Let us now introduce the volume single-and double-layer integral operators, respectively, denoted by L and M , that are defined for ρ ∈ L 2 (Γ) by The wave fields v and v − (see (3) and (4)) can be expressed as Furthermore, the single-and double-layer integral operators provide some outgoing solutions to the Helmholtz equation [20].

Direct Boundary Integral Equations for the Dirichlet Problem
The aim of this section is to provide without any detail the standard integral equation formulations for solving the two-dimensional scattering problem with Dirichlet boundary condition that will be used later. More details can be found in [17,20,28] for the derivation and properties of these integral equations (well-posedness, existence of resonant modes,. . . ). Let us introduce the following boundary integral operators Then, based on the expressions of the trace and normal derivative trace of the volume single-and double-layer potentials (see, e.g., [20] for further details), a first formulation is based on the trace of the single-layer operator where the unknown density ρ is an element of L 2 (Γ). The equation is well-posed and equivalent to the exterior scattering problem (1) as soon as k is not an irregular interior frequency of the associated Dirichlet boundary-value problem [17,20,28]. This integral equation is called Electric Field Integral Equation (EFIE) in electromagnetism. In the sequel, this formulation will be denoted by Single-Layer Integral Equation (SLIE). In the case of a closed boundary Γ, which is the situation in the paper, it can be proved that the spurious internal modes do not radiate. Therefore, there is no pollution in the far-field computation, which justifies that the SLIE can be considered as a reference solution. In addition, for an open surface Γ, the SLIE is the only possible integral equation that can be written. Another surface integral formulation is given by (see [20] for the functional framework associated to the integral equations) It is also well-posed and equivalent to the exterior scattering problem (1) if k is not an interior Neumann resonance [17,28]. This formulation is often designated by Magnetic Field Integral Equation (MFIE). Nevertheless, this integral equation is not recommended in practice since the spurious modes radiate and introduce some errors when computing the far-field pattern. To avoid the interior resonance problem, Burton and Miller [17,19,28] proposed to rather use a linear combination between the EFIE and MFIE. If α is a realvalued parameter such that 0 < δ < 1 and if η is a complex parameter with (η) = 0, where (η) is the imaginary part of η, then one gets the Combined Field Integral Equation (CFIE) [17,23,28] This integral equation is well-posed for any wavenumber k but can only be applied to closed surfaces.
An important point is that all these integral equations are based on the single-layer representation where the unknown surface field ρ is the physical quantity defined by the jump (set as the difference between the interior and exterior traces) To compute the far-field pattern, let us recall that we have: u = L ρ + M λ, where ρ and λ are two unknown densities. In the polar coordinates system (r, θ), the use of asymptotic expansions when r → +∞ leads to the following relation [22] ∀θ where a L and a M are the radiated far-fields for the single-and double-layer potentials, respectively, defined for any angle θ of [0, 2π] by with θ := (cos(θ), sin(θ)). In addition, the Sonar Cross Section (SCS) (in dB) is such that When using the single-layer representation (9), only a L is needed while a M is set to zero. In the paper, we will focus on the SLIE.

Numerical Approximation
To numerically solve Equation (6) (or (8)), we first introduce a polygonal interpolating surface Γ h which approximates Γ. The triangularization of Γ h is built by using n K linear segments K j of size h. Therefore, we have Γ h = ∪ n K j=1 K j . In the numerical examples, we take about 10 elements per wavelength which is enough to get a suitable accuracy for the Boundary Element Method (BEM). The notation P designates the space of complex-valued polynomials of order . In (6) (or (8)), no derivative operator is applied to the unknown ρ. However, as seen below, the OSRC method, and the integral equation technique for the Neumann problems introduce some tangential derivatives that apply to the surface fields. For this reason, we propose to use the linear BEM all along the paper, where the conformal finite element space V h of P 1 piecewise continuous functions is defined by For the Dirichlet problem, a P 0 (constant per segment) BEM could also be applied. For the numerical approximation, we naturally use the weak form of the integral Equation (6) which leads to [27] ∀q ∈ L 2 (Γ), Let us consider now that the BEM is applied. To compute the elementary interaction for two segments K 1 and K 2 of Γ h , a standard semi-analytical integration formula is used. This means that there is a first outer integration in (13) according to x on K 1 and based on a two-points Gauss-Legendre quadrature. Next, for the inner integral with respect to y on K 2 , if K 1 = K 2 , then again a two-points Gauss-Legendre quadrature is applied, an exact integration formula is used [35] when K 1 = K 2 . By using the interpolating surface Γ h and the linear approximations of both ρ and q in V h , the discrete form of (12) yields the linear system where L is the complex-valued matrix associated with the single-layer operator and M is the surface mass matrix. If n P is the number of points of the curve Γ h , then all the matrices are elements of the space M n P (C) of complex-valued matrices of size n P × n P . Indeed, assuming that Γ is a closed boundary, then the number of degrees of freedom of the linear boundary element method, i.e., the number of points n P := dim(V h ), is equal to the number of segments: n P = n K . In addition, the unknown complex-valued nodal vector ρ and the nodal incident vector u inc are in C n P . For computing the coefficients of L, some semi-analytical quadrature formula are used to integrate the kernel singularity. If one rather uses the BWIE (8), the discrete form leads to where ∂ n u inc is the nodal complex-valued vector related to the normal trace of the incident field. Finally, we have ρ = −(∂ n u + ∂ n u inc ).

First and Second-Order OSRCs
The on-surface radiation condition (OSRC) method was introduced in the middle of the eighties by Kriegsmann, Taflove and Umashankar [30]. At that time, the main idea was to develop an approximate but efficient and low memory numerical solution for scattering problems, most particularly in the high-frequency regime. Starting from local approximations of the Dirichlet-to-Neumann (DtN) operator, they were able to propose the computation of the scattered field by two-dimensional simple obstacles. Since then, the OSRC method has received much attention from many researchers and many improvements and extensions have been proposed (see, e.g., [31]). For the sake of conciseness, we restrict our presentation to the so-called first-and second-order Bayliss-Turkel-like radiation conditions [36].
The first-order OSRC (that we denote by OSRC 1 ) is given by In the above equation, (u 1 , ∂ n u 1 ) is the OSRC Cauchy data that approximate the exact Cauchy data (u, ∂ n u) on Γ. Here, let us remark that (u 1 , ∂ n u 1 ) must be understood as a notation since the OSRC is an approximate DtN map, which means that we do not a priori know if ∂ n u 1 is the normal derivative trace of a function u 1 . The function κ := κ(s) is the curvature at a point s of the surface, where s is the curvilinear abscissa counterclockwise directed along Γ. Equation (17) can be seen as a simple impedance boundary condition for the exterior domain Ω + . Various ways of deriving (17) are available in the literature. Let us mention for example [1] were formal and rigorous approaches are developed.
The first-order OSRC (17) can be improved, leading to the so-called symmetrical second-order Bayliss-Turkel condition (denoted by OSRC 2 in the sequel) In (18), the curvilinear derivative operator is written ∂ s u = ∇u · τ, where τ is the tangent vector to Γ (n · τ = 0). In the same spirit, an Engquist-Majda-like OSRC can be derived [1,31]. However, some numerical simulations in various papers show that the boundary condition (18) provides a higher accuracy. Therefore, in the present paper, we restrict our study to (18).

Numerical Approximation
To solve (18), we introduce the weak formulation for some suitable test-functions v in H 1 (Γ). The pair of approximate Cauchy data (u 2 , ∂ n u 2 ) is discretized in the finite element space V h × V h . When the BEM is introduced, then the discretization of (19) can be rewritten at the matrix level as where the functions α and β are defined for OSRC 2 by The matrices S α and M β are, respectively, the (sparse) n P × n P generalized stiffness and mass matrices associated with the functions α and β. In addition, the vectors u 2 and ∂ n u 2 are complex-valued vector fields of C n K . Finally, the second-order OSRC method requires the following stable approximation scheme for the curvature [36]: let K = (a 1 a 2 a 3 ) be a triangle whose vertices a j , j = 1, 2, 3, are points on Γ h , then, the curvature at the vertex a 2 can be approximated by where a j , for j = 1, . . . , 3, are the lengths of the edges of K ordered with respect to the increasing size. This formula is directly applied on the triangles associated with the surface mesh and built on two adjacent segments. Using Formula (21) implies that the curvature is not equal to zero at the corner of the square since the triangle used to compute the numerical curvature is not flat. For a uniform surface mesh with size h, the numerical curvature is given as √ 2/h. This allows us to reproduce the scattering phenomenon arising at the corner (see also [36]).
Since we are solving a scattering problem with given Dirichlet boundary condition, then u 2 is given by −u inc and then This means that ∂ n u 2 is simply obtained through the solution of a complex-valued sparse linear system defined by the mass matrix. Once the unknown is obtained then one can compute an approximation of the jump of the normal derivative trace by the relation Another way of writing this relation is that ρ 2 is solution to The far-field can be directly obtained through expressions (11). Even if the boundary condition (17) can be a priori applied to any scatterer Ω − , we will see during the numerical simulations that a serious loss of accuracy is observed when the scatterer presents some concave parts (see [36] and Figures 1-3). Indeed, in the concavity, the presence of multiply bounced rays cannot be modeled by local differential operators since the nature of this phenomena is nonlocal. Therefore, the aim of the next sections is to provide an original way to directly couple the OSRC formulation in the convex part of the scatterer to the SLIE restricted to the concave part to improve the accuracy.

Weak Coupling Procedure and Boundary Element Approximation
Since the OSRC solution is locally inaccurate in the non convex part of the domain, we propose to build a solution which is computed first by the OSRC method and improved thanks to the integral equation where the quality of the OSRC approximation deteriorates. To this end, let us assume that the geometry Γ can be decomposed into two non-overlapping parts Γ := Γ 1 ∪ Γ 2 , with Γ 1 ∩ Γ 2 = {a 1 ; a 2 }. The geometry Γ 1 is related to the boundary part of Γ which is convex, the complementary (where we will use an integral equation) is Γ 2 (see Figure 1). The proposed procedure uses 1) a global computation of an approximate surface density ρ 2 through the OSRC on Γ and 2) injects the restriction of this approximation to Γ 1 into the global integral equation formulation on Γ to obtain an approximate density on Γ 2 . This way of partitioning the computation allows us to solve a smaller size boundary element system and can be seen as a "one shot" computation.
In the following, if A is a global matrix on Γ, we denote by A j the extracted matrix related to the interaction between the part Γ j of the boundary and Γ , j, = 1, 2. We then have In a similar way, the restriction of a nodal complex-valued vector z to the part of the boundary Γ is written z , for = 1, 2. Then, we propose the following algorithm for the SLIE formulation: compute ρ 2 := (ρ 2 1 ,ρ 2 2 ) ∈ C n P such that (1) OSRC: extract ρ 2 1 = (ρ 2 ) 1 ∈ C n P1 −2 from the computation of Mρ 2 = g 2 , (2) SLIE: computeρ 2 2 ∈ C n P2 as the solution to L 22ρ Even if other coupling possibilities are available, this one has the advantage that no modification of an existing code is required. Let us denote by n Kj (respectively, n Pj ) the number of segments (respectively, points) of the boundary Γ j , j = 1, 2, based on the uniform mesh used for the full BEM on Γ h . Then, we have n K = n K1 + n K2 and n P = n P1 + n P2 − 2 (the two junction points a 1 and a 2 must be counted once). In step 1), we only retain the values of the OSRC solution that are not considered in the cavity, including the two endpoints a 1 and a 2 . This first step needs O(n K ) operations to solve the sparse complex-valued linear system (by using an LU factorization or a preconditioned iterative Krylov solver [37]) while the memory storage is O(n K1 ). In the second step, the solution to a complex-valued full matrix is necessary. If a full storage is considered, then the linear system solution needs O(n 3 K2 ) elementary operations and the memory storage is O(n 2 K2 ). Usually, in integral equation solvers, most particularly for high-frequency scattering, it is preferable to use a subspace Krylov solver [17,37,38] in conjunction with a fast matrix-vector product algorithm like for example the Multilevel Fast Multipole Method (FMM) [24], high-order solvers [18] or the recent direct Adaptive Cross Approximation (ACA) techniques [29]. In this case, the memory requirement is O(n K2 ) for a computational cost O(n K2 log n K2 ). In addition, it is well-known that preconditioning is a requirement to get a fast convergence. The second equation of (26) is defined through the single-layer integral operator. Then, the integral equation is a first-kind Fredholm equation which is known to be badly conditioned. To improve this, we propose a preconditioned version of the algorithm (26) based on the Caldèron relations [17] and well-adapted to fast iterative solvers but without modifying the solution of the initial coupling procedure. Let D be the discrete version of the normal derivative trace operator D and D 22 its restriction to Γ 2 . Then, the preconditioned version of (26) is given by: find ρ 2 := (ρ 2 1 ,ρ 2 2 ) ∈ C n P such that (1) compute ρ 2 1 = (ρ 2 ) 1 ∈ C n P1 −2 solution to Mρ 2 = g 2 , (2) obtainρ 2 2 ∈ C n P2 as the solution to D 22 L 22ρ The second equation of (27) is defined by a second-kind integral equation formulation which has better clustering properties for the convergence of Krylov subspace solvers. The price to pay is that each iteration of a Krylov solver requires to apply L 22 first and next D 22 . The preconditioned Caldèron SLIE coupling procedure (27) is called SLIE-OSRC 2 . Onceρ 2 := (ρ 2 1 ,ρ 2 2 ) ∈ C n P has been computed, all the usual quantities of interest can be obtained like for the SCS given by (11).
The SLIE-OSRC 2 formulation may suffer from the existence of interior resonances. Even if the spurious modes do not radiate for the SLIE, it can be better to have access to a stable formulation. To this end, the following BWIE-OSRC 2 coupling procedure could alternatively be used: computeρ 2 := (ρ 2 1 ,ρ 2 2 ) ∈ C n P such that Unlike the SLIE-OSRC 2 formulation, the BWIE-OSRC 2 approach is defined by a second-kind Fredholm integral equation and does not really need to be preconditioned when correctly choosing δ and η (see, e.g., [17]). An extra computational cost and memory storage is required since the evaluation of the double-layer potentials N 22 and N 21 are required.
In the following examples, we will report only the results related to the SLIE-OSRC 2 formulation since no resonance were met and we always plot the SCS. In addition, we do not use any iterative solver for the resulting linear systems because we are considering toys problems for a proof of concept. However, for three-dimensional problems, where the method directly extends, the difference between all the formulations may be important in terms of convergence rate. Here, we rather focus on the accuracy improvement. Further investigations are therefore needed, as well as improved formulations and implementation of higher order OSRCs.

A Numerical Example-Validation of the Procedure
As previously said, we use the SLIE-OSRC 2 formulation. The model toy problem is the following. We consider the obstacle Ω − as being composed of the square cylinder centered at the origin and with side length 2, with an inner square cavity defined by the corners (1/3, −1/3), (1, −1/3), (1, 1/3) and (1/3, 1/3). The boundary Γ 2 is then defined as the internal boundary to the cavity (blue curve on Figure 1) and Γ 1 is the complementary boundary on Γ (red curve on Figure 1). For an incident plane wave, the reference solution is given by the SLIE (14) discretized by the BEM (with 10 points per wavelength). The results are clearly improved compared to the pure OSRC approach for wave numbers k ≤ 15 and this, independently of the angle of attack θ inc . As it can be observed on the first picture of Figure 2 (Case 1, for k = 2), we almost obtain the reference solution by using the SLIE-OSRC 2 approach while a pure OSRC approximation gives a poor accuracy which deteriorates as k increases. For a higher wave number (Case 2, k = 8), the results are still accurate (see Figure 2). By increasing the frequency, we still have an acceptable solution as seen on the third example of Figure 2 (Case 3, k = 14). For both the integral equations and OSRCs, since we consider a mesh with ten elements per wavelength, the size of the matrix L 2,2 is about n K2 = 2.5k, which is more than four times less than for a pure SLIE solution.

The Neumann Scattering Problem
The Neumann scattering problem, i.e., considering (1) but with the boundary condition: ∂ n u = −∂ n u inc on Γ, can be treated similarly with the EFIE [17] based on the normal derivative trace of the double-layer potential written under its weak form (called DLIE). A hyper-singular kernel must then be integrated carefully by semi-integration techniques. For the OSRC, the method applies quite similarly and leads to the solution of a sparse linear system.
In general, the coupling procedure (called DLIE-OSRC 2 ) has proved to be less accurate than for the Dirichlet problem (i.e., for SLIE-OSRC 2 ). We think that this is due to the fact that the OSRC approach requires a higher order operator to provide a suitable fast solution during the first step of the method. We recommend to rather use the square-root OSRC developed in [2]. The resulting coupling technique will be analyzed in a forthcoming work. For low and moderate wave numbers (k ≤ 4), we still get quite acceptable results having always in mind the lower computational cost of the DLIE-OSRC 2 method. However, there is a moderate accuracy for wave numbers such that k ≥ 10 as seen on the two cases reported in Figure 3. Nevertheless, this must be counterbalanced by the lower cost of the procedure and the possibility to increase the OSRC accuracy in the convex part. In addition, we always obtain a good prediction of the main lobs where the energy is mostly radiating. As a general concluding remark, the coupling DLIE-OSRC 2 algorithm provides an interesting alternative to the DLIE. Furthermore, the method is always more accurate than a classical OSRC approach and gives a possibility of its extension to non-convex obstacles.

Conclusions
In this paper, we introduced a simple and original algorithm coupling the surface integral equation method and the OSRC technique for solving the scattering problem by nonconvex scatterers. While simple, the coupling method leads to an improved accuracy compared with the pure OSRC approach and reduces the computational cost of a direct integral equation formulation. In addition, we explain how operator preconditioning can be directly included into the formulations. The method is validated on a simple twodimensional problem as a proof of concept. In particular, the method is accurate for the Dirichlet problem, but still needs to be further investigated by using high-order OSRCs in the case of the Neumann problem. We expect that the formulation can be useful for three-dimensional high-frequency scattering problems solved iteratively by Krylov solvers with acceleration algorithms.
Author Contributions: Conceptualization, software, formal analysis, investigation, writing-original draft preparation, writing-review and editing, S.M.A., X.A. and C.C.; project administration, funding acquisition, S.M.A. and C.C. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors would like to thank the Deanship of Scientific Research at Umm Al-Qura University for supporting this work by Grant code 18-SCI-1-01-0017.
Informed Consent Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.