Fast rapidly convergent penetrable scattering computations

We present a fast high-order scheme for the numerical solution of a volume-surface integro-diﬀerential equation. Such equations arise in problems of scattering of time-harmonic acoustic and electromagnetic waves by inhomogeneous media with variable density wherein the material properties jump across the medium interface. The method uses a partition of unity to segregate the interior and the boundary regions of the scattering obstacle, enabling us to make use of specially designed quadratures to deal with the material discontinuities in a high-order manner. In particular, the method uses suitable changes of variables to resolve the singularities present in the integrals in conjunction with a decomposition of Green’s function via the addition theorem. To achieve a reduced computational cost, the method employs a Fast Fourier Transform (FFT) based acceleration strategy to compute the integrals over the boundary region. Moreover, the necessary oﬀgrid evaluation of the density and the inter-grid transfer of data is achieved by applying an FFT-based reﬁned-grid interpolation strategy. We validate the performance of the method through multiple scattering simulations. In particular, the numerical experiments demonstrate that the proposed method can handle high-contrast material properties without any adverse eﬀect on the number of GMRES iterations.


Introduction
In the last few decades, the scattering problem, particularly the computational aspects, has garnered a lot of interest from engineers, scientists, and mathematicians owing to its vast application base.Producing a meaningful numerical solution to this problem still remains computationally challenging, particularly where there is a discontinuity or high contrast in the material properties, and when a high-order accuracy is desired.
In this text, we consider the scatterer to be an open bounded set in R 2 with a smooth boundary , e ⊂ R 2 being the homogeneous exterior.The scalar scattering problem is modeled in terms of the Bergmann's equation [1], where (for the acoustic case) ψ denotes the time-harmonic wave field, κ is the wave number in the media, and ρ is the material density. 1he exterior field ψ e = ψ inc + ψ s satisfies the Helmholtz's equation ψ e (x) + κ 2 e ψ e (x) = 0, for x ∈ e , (2) where ψ inc is the incident wave, κ e is the wave number in e , and ψ s is the scattered wave satisfying the Sommerfeld radiation condition [2].The transmission conditions 1 are satisfied across the interface , where n denotes the outward normal and ρ e is the density in e .
For the two dimensional electromagnetic scattering problem where the electric field is given by E = (E x , E y , E z ) and the magnetic field is given by H = (H x , H y , H z ), the problem can be modeled through transverse magnetic (TM) and transverse electric (TE) modes [3] as below: where and μ are the electric permittivity and magnetic permeability, respectively.Note that for the two-dimensional electromagnetic problem, the material properties κ, and μ depend only on the x and y variables.
The main difficulty in numerically solving the wave scattering problem remains the fact that to describe highly oscillatory functions in large domains accurately, a certain fixed number of points are necessary to resolve a wavelength.The central challenge in directly solving the differential equation or its variational formulations lies in constructing an approximating boundary condition to replace the radiation condition, which in itself is a highly active area of research [4][5][6].Generally, the error associated with such boundary conditions typically dominates the error in the computed solution.Another approach is to solve an equivalent integral equation [2,7], which is favored for this problem because the solution of the integral equation inherently satisfies the radiation condition.There are two main challenges in obtaining the numerical solution of the equivalent integral equation of the scattering problem, first, singularity present in the kernel of the integral equation, and, second, quadratic cost per iteration of the linear solver.
A significant amount of work addressing the acoustic scattering problem (based on various formulations of the problem) can be found in the literature, to name a few [9][10][11][12][13][14], but most of these studies address the scattering problem where the density remains constant (ρ = ρ e ).To the best of our knowledge, only a limited number of works directly address the problem with a variable density ρ.In [15], Kriegsmann and Reiss have given a solution to the problem with variable ρ, assuming ρ/ρ e 1.They also assume that ρ is constant near , which, in general, need not be true.Although the authors consider variable material density in [16][17][18], the solution methodologies work only for spherical geometry.In [19], Bleszynski et al. address this problem in a three-dimensional setup and provide a fast algorithm.They also point out the difficulty of the high-contrast problem (ρ/ρ e >> 1), which they address in a later article [20].Though fast, the method in [20] is a low-order method.Recently in [21], the authors presented an optimized weak coupling of boundary element and finite element methods to solve acoustic scattering problems in three dimensions for homogeneous and inhomogeneous scatterers.In [22], Colton and Monk have reviewed the progress for the inverse problem with the variable density.This article proposes a fast high-order method to solve the problem with variable and highcontrast material properties.In particular, we investigate the performance of our solver in the high-contrast scenario considered in [20].
The rest of the paper is organized as follows."Numerical method" presents the proposed solution methodology; in particular, we discuss the numerical schemes used for high-order approximation of the integral operators and the density derivatives.Next, in "Computation examples", we present a variety of computational examples to showcase the efficacy of the proposed method.A brief discussion on the extension of the methodology for the corresponding three-dimensional acoustic scattering problem follows in "Extension to three dimensions", wherein we demonstrate its effectiveness via some numerical experiments.Finally, we summarize our conclusions in "Conclusions".

Numerical method
We base our numerical approximation of the solution of Eq. ( 5) on the Nyström methodology.Many of the principal components of the proposed method are adaption or modification of the works thereof [9,13,23,24].The method presented in this paper can be seen as an extension of the work presented in [13] for the solution of the Lippmann-Schwinger equation, which is a special case of the integral equation (5).The main difficulties in the fast and accurate computation of the integrals in (5) are the singular (when x = x) and the near-singular (when x and x are in close proximity) behavior of the kernel.A specialized integration scheme is designed to resolve these issues in conjunction with appropriate changes of variables.The required offgrid calculations are handled using Fast Fourier Transform (FFT) based high-order interpolation schemes designed in accordance with the underlying grid as described in this section below.
The resulting linear system, arising from the Nyström discretization of (5), is solved using the matrix free GMRES solver [25].In the following, we describe the numerical scheme for approximating the integral operators and the density derivatives.

Decomposition of integral operator
We split the closure of the inhomogeneity into P overlapping curvilinear patches, say, { p } P p=1 , where the p th -patch is homeomorphic to the set (0, 1) × [0, 1), via a smooth, invertible parametrization ξ p = ξ p (u, t).By constructing a partition of unity ω p (x) = ω p (ξ p (u, t)), p = 1, . . ., P subordinated to this covering, we rewrite the volume integral (6) as a sum of integrals over P patches Note that when the closure of a patch p does not intersect the the boundary of the domain, that is, when p ∩ = ∅, the corresponding function ω p and its derivatives vanish to high-order at both ends in every direction.On the other hand, for the case when p ∩ = ∅, the corresponding function ω p vanishes at both ends only in the udirection.This distinction allows us to employ different approximation methodologies for integration over these two sets of patches.In particular, for integration over the patches, where the POU functions vanish in every direction, we make use of the addition theorembased integration methodology [9,26] for an accurate O(N log N ) approximation (see "Interior integration").
To facilitate this, we classify the patches into the following two sets of patches and the corresponding indices, namely, ( I , I) = {( p , p) : p ∩ = ∅} and ( B , B) = {( p , p) : p ∩ = ∅}.Here onward, the patches in I are called the interior patches and the patches in B are called the boundary patches.For p ∈ B, we take the parameter t to describe the transverse parameter, where the boundary coincides with t = 0.In view of this, the set of functions {ω p (u, 0), p ∈ B} serves as a partition of unity for the surface , allowing us to decompose the surface integral (7) as where p := p ∩ .Thus, using the above patch classification and the surface integral decomposition, we can write the integral operator (Lψ)(x) as where and E(x) : = p∈I ω p (x) vanishes to the boundary with high-order.Further, the integrals (K p V ψ)(x) and (S p F ψ)(x) for p ∈ B can be rewritten, in the parametric space, as and respectively, where J p in (11) denotes the Jacobian of the change of variable ξ p .Clearly, in order to achieve high-order approximation of the operator (Lψ), the method requires high-order numerical integration schemes to approximate the integrals K I V ψ, K p V ψ, and S p F ψ given in (10), (11), and (12), respectively.In addition, the method requires an accurate numerical differentiation scheme to evaluate the expressions V ψ and F ψ.
Toward this end, we start by placing a uniform parametric grid and for each overlapping patch as illustrated in Figure 2.
We denote the set of all the Nyström nodes (that is, the points where the method computes the unknown) on the whole domain by G := G B ∪ G I , where G B := ∪ p∈B G B p , and G I := ∪ p∈I G I p .Thus, the number of unknowns for the iterative solver is N = Note that a separate grid placement for the surface calculation is not required; instead, the u-direction grid at t = 0 for p ∈ B does the job, restricting the total number of unknowns to N .Before we describe the numerical strategy in detail, we point out that the method also uses an "auxiliary" grid in polar coordinates to approximate the integral in Eq. ( 10) as described in detail in "Interior integration".Note that these auxiliary polar grid points do not contribute to the number of unknowns N in the iterative solver.
In what follows, we first present the approximation strategy for the surface integral (12), then extend this strategy to approximate the boundary-patch integral (11).In both cases, a local floating POU around the target point is used to isolate the singularity, which localizes the otherwise expensive singular integration calculation and facilitates the use of the FFT-based equivalent source acceleration strategy [13,27].For approximation of the integral (10) over the interior patches, the method relies on a suitable decomposition of the Green's function G e via the addition theorem of the Hankel function [9] to speed up the otherwise bulky calculation.To efficiently transfer of the data E(x)V ψ(x) back and forth between the polar grid and the uniform parametric mesh on the interior patches, we set up two fast high-order FFT-based interpolation schemes discussed later in the text.

Approximation of the surface integral
The difficulties in numerical approximation of the surface integral (12) are similar for each p ∈ B; hence it is sufficient to discuss the method for a particular p.First, when the target point x is sufficiently away from the surface p , the integrand is smooth and periodic (owing to the vanishing properties at both end).Thus, an application of the trapezoidal rule can effect a high-order approximation.On the other hand, when x ∈ p or x is close to the surface patch p , the kernel G e has a singular or a near-singular (due to the rapid change in the kernel value) behavior, respectively.To achieve a high-order accuracy without compromising efficiency, as mentioned above, the method localizes the singularity using a smooth cutoff function η u centered at u; in our computation we used for γ 1 > γ 0 ≥ 0 and d = |u − u |.Using η u we rewrite Eq. ( 12) as a sum of two integrals where ).The first integral on the right hand side (r.h.s.) of ( 16) has a smooth periodic integrand, and can be effected to a high-order approximation using the trapezoidal rule.On the other hand, the integrand of the second integral on r.h.s. of ( 16) contains a singularity.In order to resolve the singularity, we utilize a change of variable where υ(τ ) is an odd function that vanishes at τ = 0 along with M (M ≥ 1) of its derivatives.One such "change of variable" is for an even, non-negative integer M.This change of variable renders the integrand in ( 19) below smooth with uniformly bounded M derivatives, which vanishes to high-order at the endpoints τ = −τ 1 and τ = τ 1 , where τ 1 = υ −1 (γ 1 ).Thus, a high-order integration is possible using the trapezoidal rule (cf.Remark 1).Moreover, the change of variable given by ( 17) also resolves the near singularity [11].Notably, this change of variable necessitates density values at some offgrid points, which we obtain via the one-dimensional version of the high-order TrigPoly interpolation introduced in [23].

Approximation of the boundary-patch integral
The present section provides a concise description of the numerical strategy employed to approximate the boundary-patch integrals.Similar to the surface calculation, it is sufficient to discuss the methodology for a particular p ∈ B. The difficulties in numerical approximation of these integrals significantly depend on the position of the target point, leading to the following three different instances.

Non-singular integral
In case, the target point x resides outside the integrating patch p , the integrand is smooth and vanishes at two ends in the u-direction; see Figure 1.A combination of the trapezoidal rule in the u-variable, and a high-order Newton-Cotes quadrature in the transverse t-variable, effects a rapid convergence in this case.

Singular integral
Second, we consider the case where the target point x resides within the integrating patch p and it coincides with a quadrature point x p ij in G B p (see (13)).Here, we encounter a singularity at ξ p (u , t ) = x and a near singularity when these points are in close proximity.In order to treat the singularity, the method utilizes a one-dimensional floating POU η x (u ; t ) in the u-variable centered at the target point x, leading to where V p ψ(u , t ) = ω p (ξ p (u , t ))V ψ(ξ p (u , t ))J p (u , t ), and The first integral on the r.h.s. of ( 20) is free of singularity and can be treated to a highorder approximation following the non-singular case described above in "Non-singular integral".On the other hand, I p (t ; x, u) has a corner type singularity at t = t, which can be resolved by splitting the second integral in the r.h.s. of (20) at t = t [13], where each integral piece is approximated to high-order using Newton-Cotes quadrature.The values of I p (t ; x, u) as required at the transverse Newton-Cotes quadrature points are obtained in a high-order manner using the trapezoidal rule in conjunction with the change of variable u = u (τ ) given by Eq. ( 17).

Offgrid boundary-patch evaluation
Finally we consider the case, when a target point resides in the integrating patch, but not a quadrature point, that is, x ∈ p and x / ∈ G B p ; which only happens when the target point comes from the overlapping region.Such offgrid evaluations are dealt by interpolation of the precomputed integral values Using the properties of a partitions of unity {W p : p ∈ B} set over the boundary patches (where the functions W p vanishes only in the u-direction), we can express and the precomputed values of K B V ψ can be used to set up an FFT-refined polynomial interpolation [14] for every boundary patch making accurate and efficient evaluations of K B V ψ at offgrid points.

FFT-based equivalent source acceleration of the non-singular surface and boundary-patch integrals
It is straightforward to see that while the method discussed above for the non-singular integrals converges with high-order, it has a O(N B N ) + O(NN S ) computational complexity, where N B and N S are the number of the boundary-patch and the surface unknowns, respectively.To improve the complexity to O(N log N ), the method utilizes the two-face equivalent source acceleration strategy introduced in [27] for three dimensional surface computation and later adapted for volumetric calculation in [11].Below we provide a concise description of this acceleration strategy.
We consider a square C of side A containing the domain , which is then further partitioned into L 2 identical non-overlapping squares C ij (i, j = 1, . . ., L) of side H = A/L such that the cells C ij do not admit inner resonance.To accelerate the non-singular integrals, we split the boundary-patch and the surface integral in (9) into two, namely, "adjacent" and "non-adjacent" interactions.To precisely define the adjacency, for a grid point x ∈ C ij which we refer as a "true source" point, we first define the neighbourhood set N (x) which contain at most eight neighbouring cells.Two source points x and x are considered to be adjacent if x ∈ N (x), and to be non-adjacent if otherwise.Utilising the definition of neighbourhood set and by ensuring that the support of the cut-off functions η x (u ; t ) and η u (u ) (see (15)) are contained within the neighbourhood set N (x), the boundarypatch and the surface integrals can be rewritten as sum of the adjacent and non-adjacent interaction as follows: The adjacent interactions above in Eqs. ( 22) and ( 23) can be approximated to high-order by employing the quadrature schemes discussed above for the non-singular integrals in O(N ) complexity.The non-adjacent interactions, free of the cutoff functions, take a simpler form of the discrete convolution that reads w l G e (x, x l )ϕ(x l ), (24) where ϕ(x l ) = F ψ(x l ) + V ψ(x l ) if x l ∈ G B ∩ , and ϕ(x l ) = V ψ(x l ) otherwise.We point out that the true source values in (24) includes contributions from the non-adjacent surface calculations.Note that the convolution in (24) can not be computed directly using FFT as the trues source points x l ∈ G B are irregularly distributed.We also note that while the contributing sources range only in x l ∈ G B , but the convolution (24) needs to be computed at additional points which will not lie in the boundary region.We call these extra points as "virtual source point", and include them in summation (24) with zero weight w l , which facilitates an efficient evaluation of the non-adjacent interactions at all these virtual source locations using FFT-based acceleration strategy that we briefly describe next for completeness.
The key idea of the acceleration technique is to compute (24) by replacing the true sources with a certain set of "equivalent sources" distributed uniformly on the two parallel faces of the cell C ij (see Fig. 3).Let σ m ij,l and σ d ij,l (l = 1, • • • , N eq ) denote the acoustic monopoles and dipoles at the Cartesian grid x ij,l , respectively.Then, the field induced by these equivalent sources is Let K reg,true ij,B (x) denotes the field induced by true sources, say, given by The unknowns σ m ij,l and σ d ij,l in (25) are obtained by solving a over determinant linear system Aσ = b, where the matrix A is obtained by evaluating (25) at n coll = 4N eq points located on the boundary of the adjacent set N (x), and the vector b is obtained by evaluating (26) at n coll points.It is well known (see [27]) that the quantity K reg, eq ij,B (x) provides an accurate approximation to the field K reg,true ij,B (x) generated by true sources in cell C ij .Owing to the identical geometry of each cell C ij , the QR decomposition of matrix A is obtained once and stored for repeated use.This whole process requires O N 3/2 /L 3 + O N 3/2 /L 2 operations in total, see [13].
For any x ∈ C ij , the non-adjacent interaction induced by all Cartesian grids can be obtained by subtracting the adjacent interaction from the whole convolution.To be more precise, for any x ∈ C ij , the quantity provides an accurate approximation of the non-adjacent interactions defined in (24).The convolution on the right-hand side is described over an equidistant Cartesian grid, hence can be computed using FFT.Finally, to obtain the field values at the true source and the virtual source locations within each cell C ij , a Dirichlet problem for the free space Helmholtz equation is solved.Owing to the non-resonance of each cell this problem is uniquely solvable.The solution to all these well-posed problems can be obtained efficiently by discretized the plane wave expansion in a coordinate system local to cell C ij of the form where d l is the unit vector on the surface of the unit disc centred at x ij,c , the centre of the cell C ij .The coefficients β l are computed by solving an overdetermined linear system constructed by matching (28) on the Cartesian grid points of the parallel faces of cell C ij .It was shown in [13] that the acceleration strategy has a computational complexity of O(N log N ).

Interior integration
We first note that the integration over the interior patches can be re-written as where B R ⊃ is a circular region of radius R containing and E(x) : = p∈I ω p (x) vanishes to the boundary .As mentioned earlier in "Numerical method" to approximate the integral over B R in Eq. ( 29), we employ a quadrature that relies upon decomposition of the kernel.In particular, using the addition theorem for the Hankel function H 1 0 (κ e |x − x |) [2], in polar coordinates x = (r cos θ, r sin θ) and x = (r cos θ , r sin θ ), the Eq. ( 29) yields where, denoting the -th order Bessel and Hankel function of the first kind by J and H 1 , respectively, and J κ e , (r, r ) = i 4 H 1 (κ e max(r, r ))J (κ e min(r, r )).Clearly, we have a logarithmic singularity as r → 0 and a jump discontinuity at r = r due to the presence of J κ e , (r, r ) in the integrand.Note that, owing to the presence of the factor E(x), the integrand EV ψ(r, θ) is smooth in R 2 .
For any fixed r ∈ [0, R], the evaluation of the angular integral amounts to computing the th Fourier coefficients of a smooth and 2π periodic function.An accurate and highorder approximation can be effected using the trapezoidal quadrature rule, which can be accomplished by one time application of the FFT.Thus, we place a uniform grid along the angular direction whereas we place a piecewise Chebyshev grid (see Fig. 2) for the radial integration as described in below.
Denoting the th Fourier coefficient of g by (g) , we have where  31) and (32) as where the integral moments need to be computed only once at the beginning of each run and stored for future use as the integrands are known analytically.Storing these moments requires O(N c N r FN i ) memory, where N r = N c +1 (see [26]).One can use a high-order quadrature, for example, Clenshaw-Curtis quadrature, to approximate the moments in Eqs. ( 37) and (38).Once the moments are available, an accurate approximation of the radial integral can be obtained by adding and subtracting scaled values of these moments.The integral I 2,log can be approximated efficiently by first, resolving the logarithmic singularity present in Y log using integration by parts and then, utilizing the approximated values of I 1 (r j ) [9].
Transfer of data between the polar grid and the uniform parametric grid.The method uses two FFT-based interpolation strategies to evaluate the density V ψ on the polar grid as required by the integration strategy described above in this section and to transfer the computed integral values to the solver on the uniform parametric grid.First, (for interpolation of V ψ onto the polar grid)-since the POU is smoothly vanishing in every direction for each interior patch, we use an FFT-refined local polynomial interpolation strategy in both u-and v-directions.And secondly, (for interpolating the computed integral values K I V ψ onto the parametric grid)-we set a piecewise Chebyshev interpolation along the radial direction and a refined FFT-based polynomial interpolation along the angular direction.Since the local polynomial interpolation scheme has an O(1) cost, the interpolation strategies described above have the same computational complexity of O(N log N ) as the FFT procedure.

Computation of density derivatives
The last important component of the algorithm we discuss is the computation of the density derivatives.We have utilized the Lower Degree Chebyshev (LDC) method proposed in [24] for efficient computation of the density derivatives without significant degradation in accuracy.The basic idea of the LDC algorithm is to first, approximate the function via a low order Chebyshev polynomial expansion from the given approximated data, then compute the Chebyshev coefficients for the derivative iteratively.The data on the corresponding Chebyshev grid point are again computed by use of FFT refined polynomial interpolation.A suitable degree is used in the polynomial interpolation to obtain desired accuracy in the density values at the underlying Chebyshev points.Once these values are obtained, the Chebyshev coefficients are evaluated by one application of Discrete Cosine Transform (DCT) and the Chebyshev coefficients for the derivative are then obtained by the recurrence formula [28].For a detailed discussion and rigorous theoretical analysis of the LDC algorithm we refer the readers to [24].
Remark 1 (Expected convergence rate) The rate of convergence depends on several parameters including the degree of interpolating polynomial, the degree of approximating Chebyshev polynomial, the order of the Newton-Cotes quadrature, the value of M (cf.Eq. 18) etc.

Computational examples
This section demonstrates the performance of the method through several computational examples.The relative error reported in the tables are obtained using the following formulae and , where ψ exact (x) and ψ approx (x) denotes the analytical/refined-grid solution and the computed solution at the point x, respectively.The numerical "Order" of convergence presented in the tables in this document is obtained as-Order = log 2 ( N / 2N ), where N denote the relative error corresponding to N number of discretization points."Iter" denotes the number of iteration required for the linear solver to converge up to the given tolerance level.Unless mentioned otherwise, in all the examples presented in this section, we used -a three patch covering of the obstacle, that is, P = 3 in (8), -a 5-point Newton-Cotes quadrature for integration over the boundary-patch region, and -the plane wave incident field exp(iκ e x • d) with d = (1, 0).All of the numerical experiments reported in this paper were run on a computer with 2.6 Ghz Intel Xeon Processor with 256 Gb of memory.The two-dimensional results were run on a single thread, whereas for the three-dimensional computations multi-threading was used.
Remark 2 Note that the sizes of the polar grid are suppressed in the convergence studies as addition theorem based calculation in this case, due to vanishing density EV ψ, converges super-algebraically [29] and significantly faster than the convergence rate of the other components of the method.The size of the polar grid need not be increased with the same ratio as the parametric grid to maintain the polynomial convergence rate.
Example 1 (Forward map convergence) For our first numerical experiment, we consider the convergence study of the forward map computation (i.e.convergence study for the Table 1 Convergence study for the forward map of the operator Lψ for a unit disc with κ/κ e = √ 2.0 and ρ/ρ e = 1.0 using 5-point Newton-Cotes quadrature for integration in the t-direction over the boundary region proposed quadrature for approximation of the operator Lψ).In this example, we consider a unit disc scatterer centered at the origin with constant material properties, namely κ = √ 2, κ e = 1, and ρ = ρ e = 1.0 in which case the integration can be evaluated analytically.The convergence study presented in the Table 1 demonstrate a high-order convergence rates.
Example 2 (Scattering by a lossy medium.)For the second exercise, we consider a circular inhomogeneity of radius one centered at the origin with lossy medium (modeled using complex valued material properties), namely, κ = √ 3 + i √ 2 with κ e = 1, ρ = 2 and ρ e = 1.Table 2 presents a convergence study by comparing the approximated solution against the Mie series solution.The results presented in Table 2 demonstrate a high-order convergence of the approximated solution.
Example 3 (High contrast medium) For this experiment, we consider an example with a high-contrast between the material values across the interface of the unit disk scattering obstacle.In particular, we consider κ/κ e = 10 −04 (with κ = 10 −04 and κ e = 1) and ρ/ρ e = 10 3 ( with ρ = 10 3 and ρ e = 1); the same ratios as in the article [20].We present the convergence study in Table 3.A quick comparison between Tables 2 and 3 shows that the two scenarios presented by these tables require a similar number of GMRES iterations to converge and produce similar errors in the solution.This exemplifies the robustness of the method presented above.
Example 4 (Tempered growth in computational cost.)We showcase the efficiency of the method by demonstrating a tempered growth in the computational cost of the proposed algorithm while maintaining a fixed error level by comparing the approximated operator against the continuous operator in Eq. (5).In every step while we double our grid size, we also double the wave number κ e to keep the number of points per wavelength unchanged and avoid any deterioration in the accuracy level.The times required per iteration, as  Example 5 [(High contrast medium for (relatively) complex geometry.)]In this example, we consider a bean shaped scatterer described by the parametric representation with a = 0.65, b = 1.5, r = 1 and 0 ≤ θ ≤ 2π.In this experiment, we consider κ e = 1, κ = 10 −4 , ρ e = 1 and ρ = 10 3 .The convergence study presented in Table 5 demonstrates a high-order convergence.
Example 6 (Variable material properties) So far in this article, we have considered only constant material properties inside the scatterer, although, as mentioned previously it is not necessarily so.Toward this, we demonstrate the convergence results for variable material properties in Tables 6 and 7 for the bean shaped scatterer (given by Eq. 39).In the convergence study presented in Tables 6, only the wave number κ varies within the scatterer while the interior density ρ remains constant.In Tables 7, on the other hand, both the wave number and the density vary within the scatterer.
Example 7 (Graphical presentation) In this example, we visualize the scattering for the plane wave incident field exp(iκ e x • d) with d = (1/ √ 2, 1/ √ 2) by a bean shaped scatterer given by (39) in the setting of Example 5. We consider a lossy medium with exterior wave Table 6 Convergence study for scattering by a bean shape with κ i = sin(π x 1 ) cos(π x 2 ) and ρ i /ρ e = 5, where 5-point Newton-Cotes quadrature is used for integration in the t-direction Table 7 Convergence study for scattering by a bean shaped scatterer with κ e a = 6.0, κ i = sin(π x 1 ) cos(π x 2 ) and ρ i = 2.0 + sin(πx 1 ) cos(π x 2 ) where 5-point Newton-Cotes quadrature is used for integration in t-direction over the boundary region  Example 8 (Variable material properties) As a last example of this section, we simulate the scattering by a bean shaped inhomogeneity given by Eq. (39) with variable material properties for a plane wave incident field exp(iκ e d • x) with d = (1, 0) and acoustic size κ e a = 20.We consider the interior wave number κ = (sin(πx 1 ), cos(πx 2 )) and density ρ = 2.0 + sin(πx 1 ) cos(πx 2 ).The real and absolute value of the computed solution are plotted in Figure 5.

Extension to three dimensions
The solution methodology presented in "Numerical method" for the two-dimensional scattering problem has an algorithmic extension that allows for an approximation of the solution to the scattering problem in three dimensions.The difficulties in high-order approximation of the integral operators and evaluating the density derivative in three dimensions are largely analogous with Spherical Harmonic Transformation replacing the Fourier to approximate the angular integral (cf.Eq. 46) in the interior computation.Below we succinctly describe some of the salient points of the extension.As in the two-dimensional case, we begin with a decomposition of the scatterer into boundary and interior patches.For instance, Figure 6 provides an illustration for a three patch decomposition of a bean shaped domain.Subsequently, the integral operator is decomposed over the boundary and the interior patches as in (9).We note that one can write the integral (K p V ψ) over a boundary patch (i.e., p ∈ B) as where ξ p = ξ p (u v , t ), V p ψ(ξ p ) = V ψ(ξ p )J p (ξ p )ω p (ξ p ) J p is the Jacobian of the transformation ξ p : p → (0, 1) 2 × [0, 1).Again, {ω p (x| t=0 ), p ∈ B} serves as a partition of unity for the surface .Thus, the planar surface integral can be written as where F p ψ(ξ p ) = F ψ(ξ p )ω p (ξ p )J s p (ξ p ), and J s p is the surface Jacobian.In order to approximate the integral (42), we again segregate the singular integral using a floating POU centered at the target point with a circular support.The method then approximates the regular part of the surface integral using the trapezoidal rule in u and v variable; whereas the singular integral over the localized circular domain is treated with first, going to the polar coordinates (r, θ) centered at (u, v), and then a polynomial change of variable r = τ M .The non-singular integral over a boundary-patch is approximated using a combination of the trapezoidal rule in the (u, v) variables, and a high-order Newton-Cotes quadrature in t variable.To approximate the singular integral, first, we change to polar coordinates (r, θ) centered at projection of ξ −1 p (x) onto the integration plane (u , v , t ), and then make a polynomial change of variable r = τ M to resolve the singularity.The corner singularity in t-direction is again treated by splitting the integral in two at t = t ; see [11] for detail.
For approximation of the interior integral K I V ψ, the method relies upon use of the addition theorem [30] (in spherical coordinates) of Green's function of the Helmholtz's equation in three dimensions [31].The truncated sum approximates the operator with high-order [31], where (K I V ψ) m (r) = iκ e R 0 h (1) (κ e max(r, r ))j (κ e min(r, r ))(V ψ) m (r )r 2 dr , (45) and Y m are the spherical harmonics of degree and order m with h (1) and j denoting the first kind spherical Hankel and Bessel functions, respectively.The integrals (V ψ) m (r ) can be seen as the spherical harmonic coefficients of V ψ, and can be computed by an application of the spherical harmonic transforms [31].The radial integral can then be  expressed in the form of Eq. ( 30) using spherical Bessel functions and the spherical harmonic coefficients (V ψ) m (r ).Finally, the radial integral can be approximated by expanding the spherical harmonic coefficients in Chebyshev series and utilizing the precomputed integral moments as in the two-dimensional case; see [31].
In order to demonstrate the high-order convergence of the three dimensional counterpart of the method, we present a convergence study in Tables 8 for the scattering of plane wave ψ inc = exp(iκ e x 3 ) by a penetrable spherical inhomogeneity.The material properties for this study are-κ = 10 −04 , κ e = 1, ρ = 10 3 and ρ e = 1.The errors presented in Tables 8 are obtained by comparing the computed solutions against the Mie series solution.The study clearly exhibits rapid convergence of the computed solutions.
As a final example, we consider simulation of scattering by a bean shaped obstacle in Figure 7 of the incident wave ψ inc = exp(iκ e x 3 ).The boundary of the scatterer is given by The largest dimension of the geometry is 2R, along the x 3 -axis.The absolute value of the computed solution is displayed in Fig. 7.

Conclusions
In this paper, we present a fast high-order method for solving an integro-differential equation that arise in problems of wave scattering by penetrable inhomogeneity with variable density with possible jump in the material properties across the boundary of the geometry.High-order convergence is achieved using a specialized quadrature for integration over boundary patches in conjunction with the use of addition theorem (for the Green's function) based approximation of the interior integration.The performance of the method is demonstrated through a variety of numerical examples including the scattering problem for high-contrast medium.

Fig. 1
Fig. 1 Illustrations of the curvilinear patches and the corresponding POU functions for a three patch cover for a bean shaped scatterer.a Shows the scatterer shape.b Illustrates the POU function for the interior patch whereas (c) and (d) illustrates two boundary patches.e, f Illustrate the collection of both the boundary patches and all of the patches, respectively

Fig. 2
Fig. 2 An illustration of the grid placement for a three patch decomposition of a bean shaped domain (given by, e.g., the blue curve in the top-right subfigure.)The figure shows on the top left-the grid over two boundary patches; at the top right-grid over the interior-patch region; at the bottom left-union of the grids over the boundary-and interior-patches; and at the bottom right-the auxiliary polar grid encompassing the domain

Fig. 3
Fig. 3 An illustration for the acceleration of the non-singular calculations of the boundary-patch integrals, namely, the boundary-patch non-adjacent interactions presented in "FFT-based equivalent source acceleration of the non-singular surface and boundary-patch integrals".The locations x ij, of the equivalent sources on the parallel faces of cells C ij are depicted by solid black circles.The true source points are shown as red diamonds.The virtual source points (with zero weights) are shown as green stars.The image also depicts the neighbourhood N (x) for a target point x

Fig. 7
Fig. 7 Visualization of the simulation of scattering by a bean shaped scatterer for the plane wave incident field exp(iκ e x 3 ) propagating along x 3 -axis with κ/κ e = 10 −4 with κ e = 5 and ρ/ρ e = 10 3 with ρ e = 1.Acoustic size of the scatterer is 10

Table 2
Convergence study for scattering in the setting of Example 2. The 5-point Newton-Cotes quadrature is used for integration in the t-direction over the boundary region

Table 3
Convergence study for scattering by a disc centered at origin with radius 1

Table 4
Tempered growth in computational cost in the setting ofExample 4

Table 5
Convergence study for scattering by a bean shaped scatterer of acoustic size κ e a = 3, with κ/κ e = 10 −4 and ρ/ρ e = 10 3 in the setting of Example 5

Table 8
Convergence study for scattering by a sphere of unit radius and centered at origin, with κ i /κ e = 10 −4 and ρ i /ρ e = 10 3