Abstract

In this paper, we consider acoustic or electromagnetic scattering in two dimensions from an infinite three-layer medium with thousands of wavelength-size dielectric particles embedded in the middle layer. Such geometries are typical of microstructured composite materials, and the evaluation of the scattered field requires a suitable fast solver for either a single configuration or for a sequence of configurations as part of a design or optimization process. We have developed an algorithm for problems of this type by combining the Sommerfeld integral representation, high order integral equation discretization, the fast multipole method and classical multiple scattering theory. The efficiency of the solver is illustrated with several numerical experiments.

© 2014 Optical Society of America

1. Introduction

The problem of designing composite materials that exhibit a specific acoustic or electromagnetic response is an area of active research [13]. Examples include the design of random media with a well-defined macroscopic refraction (coherent scattering) [2] and the fabrication of metamaterials [3] for cloaking, near field imaging, etc. In many experiments, the materials are designed by incorporating large numbers of identical inclusions (particles) in a layered material. When the size of each particle is comparable to the wavelength of the incoming field and the distribution of particles is reasonably dense, then the interaction of the particles involves non-negligible multiple scattering effects and methods based on homogenization [2] are not applicable. Instead, the full Helmholtz or Maxwell equations should be solved at each iteration of the design process. Numerical simulation, in the absence of suitable fast algorithms, are impractical when thousands of particles are involved.

In this paper, we develop an algorithm that accelerates the computation of electromagnetic scattering when a large number of particles are embedded in the middle of a three-layer dielectric medium. Numerical experiments show that our solver takes 1–2 minutes to evaluate the scattered field for up to 5, 000 particles on a 2.3GHz laptop. Our method combines the Sommerfeld integral representation, a well-posed integral formulation, high-order discretization, multiple scattering theory and the fast multipole method. We focus on the two dimensional setting by assuming the material is invariant in the z direction. A related three-dimensional solver was considered in [4], but the particles were assumed to be distributed in free space. A principal contribution of this paper is the development of a mathematical framework that permits them to be embedded in a layer material (which is closer to being manufacturable). While we restrict our attention here to the three-layer case, the extension to an arbitrary layered medium is straightforward.

More precisely, we consider time-harmonic scattering (with time dependence eiωt) from a three-layered medium as depicted in Fig. 1. The incident field is assumed to be driven by a point source located in the first (top) layer. The thickness of the middle layer is denoted by d. We assume the magnetic permeability μ is identical in each layer, while the electric permittivity ε is piecewise constant. There are two fundamental polarizations in the two dimensional setting to consider: the transverse magnetic (TM) polarization and the transverse electric (TE) polarization. In both cases, the Maxwell equations reduce to a scalar Helmholtz equation. For simplicity, we consider the TM polarization here, in which case the scattered field us must satisfy the equation

Δus+k2us=0,
where k=ωεμ is the wavenumber. We denote by k1, k2 and k3 the wavenumber for the three layers, and by kp the wavenumber for the particles. The scattered field also has to satisfy the Sommerfeld radiation condition at infinity [5]:
limrr(usnikus)=0,
where r=x2+y2.

 

Fig. 1 Geometry of the three-layered medium, with a large number of dielectric particles embedded in the middle layer.

Download Full Size | PPT Slide | PDF

In order to develop an especially fast solution method, we make two further assumptions. First, as in the fast multi-particle scattering (FMPS) method of [4], we assume that the particles are well separated from each other - that is, the separation between particles is at least 10% of the particle size. Second, we assume that only a finite number of distinct particle shapes are included in the simulation. The first condition ensures that a multiple scattering formalism will be accurate and the second condition ensures that precomputation of single particle scattering matrices permits a dramatic reduction in the number of degrees of freedom necessary for the solver. The particles are not assumed to be symmetric and may be placed with arbitrary orientation. Both hypotheses are common in materials design (although there are exceptions).

An outline of the paper follows. In Section 2, we introduce the Sommerfeld integral and its application to layered materials (in the absence of inclusions). In Section 3, we review classical multiple scattering theory for circular particles. Section 4 extends the scattering formalism to non-circular particles and Section 5 develops analytical tools needed to go back and forth between the Sommerfeld integral formalism and multiple scattering theory. Section 5 also combines the techniques in the preceding sections and extends the FMPS method to layered media. Numerical examples are provided in Section 6 to illustrate the efficiency of the method, followed by some concluding remarks in Section 7.

2. The Sommerfeld integral for layered media

Wave propagation in a layered medium is a well-studied problem in acoustic and electromagnetic scattering theory. Nearly a century ago, Sommerfeld developed a spectral representation involving a Fourier integral in the “transverse” variable (the x-coordinate in Fig. 1) [6]. Assuming a point source is located at x0 = (x0, y0) in the top layer, with wavenumber k1, the corresponding field is given by the (two-dimensional) free space Green’s function: Gk1(x,x0)=i4H01(k1|xx0|), where H01(x) is the first kind Hankel function of order zero. Combing the Fourier transform and contour integration [7], the Green’s function can also be written in the form:

Gk1(x,x0)=14πeλ2k12|yy0|λ2k12eiλ(xx0)dλ.
It is important to note that the Sommerfeld integral (3) is conditionally convergent and as stated, requires that yy0.

In the Sommerfeld approach [6], we assume the upward scattered field u1s in the top layer can be expressed as

u1s=14πeλ2k12yλ2k12eiλ(xx0)σ1(λ)dλ,
where σ1(λ) is an unknown density on the upper interface y = 0. It is straghtforward to verify that u1s satisfies the Helmholtz equation with Helmholtz parameter k1.

In the second layer, the scattered field u2s can be written in terms of contributions from both the upper (y = 0) and lower (y = −d) interfaces: u2t and u2b. These are given by

u2t=14πeλ2k22yλ2k22eiλ(xx0)σ2+(λ)dλ,
u2b=14πeλ2k22(y+d)λ2k22eiλ(xx0)σ2(λ)dλ,
where σ2+(λ) and σ2(λ) are used to denote spectral density functions on the upper and lower interfaces.

Similarly, we can represent the scattered field u3s in the third layer with an unknown density σ3(λ) on the lower interface as

u3s=14πeλ2k32(y+d)λ2k32eiλ(xx0)σ3(λ)dλ.
Remark 1. The signs of the terms e±λ2ki2y and e±λ2ki2(y+d) in Eqs. (4)(7) ensure that evanescent modes (when |λ| > |ki|) decay away from each layer. (Physically, this is related to causality and is required in the derivation of formula [7] by contour integration).

It is worth noting that the four unknown densities σ1, σ2+, σ2 and σ3 can be interpreted in two ways. First, they can simply be considered the spectral densities in the Fourier domain of a consistent representation for the Helmholtz equation. For those more familiar with potential theory, they can be viewed as the Fourier transforms of charge densities of four single layer potentials lying on the corresponding interfaces [8].

In the absence of any inclusions, the Sommerfeld representation for the field in each sub-domain is derived from a “mode by mode” analysis. That is, the unknown functions σ1, σ2+, σ2, and σ3 are found by enforcing the continuity conditions at the interface for each value of the argument λ. For the case of electromagnetic scattering in TM polarization, when the permeability μ is constant in each layer, this requires that

[u]=0,
[un]=0,
where [·] denotes the jump of a function along the interface, ∂/∂n is the normal derivative and u is the total field in each layer [5].

It is straightforward to check that the linear system to be solved for each λ takes the form:

(1λ2k121λ2k22eλ2k22dλ2k2200eλ2k22dλ2k221λ2k221λ2k3211eλ2k22d00eλ2k22d11)(σ1(λ)σ2+(λ)σ2(λ)σ3(λ))=(eλ2k12y0λ2k120eλ2k12y00)
Definition 2.1. We will denote the 4 × 4 matrix above by Aλ.

For the problem we consider here, the Sommerfeld integrals must be coupled to a representation of the field induced by the many particles present in the central layer. Before discussing the coupled system, however, we first summarize some well-known facts about scattering from a finite collection of inclusions in a homogeneous infinite medium.

3. Wave scattering for disks

Suppose now that we have an inclusion of dielectric material with kp=ωεpμ embedded in R2, assumed to consist of a dielectric with k2=ωε2μ. For transverse magnetic(TM) polarization, the total electrical field u in the exterior of the inclusion satisfies the Helmholtz equation:

Δu+k22u=0.
Further, the total field u can be written as the sum of the incident field uinc and the scattered field us, where us satisfies (11) and the Sommerfeld radiation condition,
limrr(usnik2us)=0,
where r=x2+y2. Within the inclusion, the field u satisfies the Helmholtz equation with wavenumber kp,
Δu+kp2u=0.
On the boundary of the inclusion, we must enforce the continuity conditions given by Eqs. (8) and (9).

3.1. A single disk

When the inclusion is a disk of radius R centered at the origin, it is straightforward to represent the solution using separation of variables, with

us=n=βnHn(k2r)einθ
in the exterior and
u=n=γnJn(kpr)einθ
in the interior. Here, (r, θ) are the polar coordinates of a point in the plane, Hn(r) is the Hankel function of the first kind of order n and Jn(r) is the Bessel function of order n [9, 10].

We now expand the incident wave uinc and its normal derivative in the form:

uinc=n=αnJn(k2r)einθ,uincr=n=αnk2Jn(k2r)einθ.

Enforcing the continuity conditions (8), (9) on the boundary of the disk for each Fourier mode, we easily obtain the following linear equation for mode n:

[Hn(k2R)Jn(kpR)k2Hn(k2R)kpJn(kpR)][βnγn]=[αnJn(k2R)αnk2Jn(k2R)],
where n ∈ ℕ.

Solving Eq. (17) determines the coefficients βn, γn:

βn=[kpJn(k2R)Jn(kpR)k2Jn(k2R)Jn(kpR)k2Hn(k2R)Jn(kR)kpJn(kpR)Hn(k2R)]αn,
γn=[k2Jn(k2R)Hn(k2R)k2Jn(k2R)Hn(k2R)k2Hn(k2R)Jn(kpR)kpJn(kpR)Hn(k2R)]αn.
It is straightforward to verify that the denominator in the preceding expressions cannot vanish if k and kp have positive real part and non-negative imaginary part [5, 11].

Definition 3.1. The mapping between the incoming coefficients {αn} and outgoing coefficients {βn} is referred as the scattering matrix for the disk and denoted by S.

Remark 2. While we restrict our attention here to dielectric particles, the method can easily be extended to perfectly conducting disks. Since the interior field u in (15) is zero for perfect conductors, from (17), we have

βn=Jn(k2R)Hn(k2R)αn.
Remark 3. In the remainder of this paper, we will refer expansions based on Hankel functions, such as (14), as multipole expansions or H-expansions, and expansions based on Bessel functions, such as (15), as local expansions or J-expansions.

Remark 4. In practice, we will truncate the expansions after, say, p terms with the value of p to be determined later. We then define α⃗ ≡ (αp, αp+1,..., α0, α1,..., αp) and β⃗ ≡ (βp, βp+1,..., β0, β1,..., βp).

3.2. Multiple disks

Suppose now that we have M well separated, identical dielectric disks randomly distributed in a homogeneous medium. Each disk is assumed to have radius R and wavenumber kp and the background medium again has wavenumer k2. For each individual particle, the analysis can be carried out as above. We will denote by α⃗m the incoming coefficients and by β⃗m the outgoing coefficients for the m-th particle. We have

βm=Sp[αm],form=1,,M.
where Sp denotes the truncated (2p + 1) × (2p + 1) scattering matrix acting on the truncated expansion.

The principle difference between the single particle and multi-particle scattering problem is that, in the latter case, the incoming field experienced by each particle consists of two parts: the (applied) incident field uinc and the contribution to the scattered field us from all of the other particles. In order to formulate the problem concisely, given the multipole expansion for disk j, we need some additional notation.

Lemma 3.1. [12] Let disk m be centered at xm and let disk l be centered at xl. Then the multipole expansion

n=βnmHn(k2rm)einθm
induces a field on disk l of the form
u=n=αnlJn(k2rl)einθl
where
αnl=n=ein(θlmπ)βnnmHn(k2xmxl).
Here, (rm, θm) and (rl, θl) denote the polar coordinates of a target point with respect to disk centers xm and xl, respectively and θlm denotes the angle between (xmxl) and the x-axis.

Remark 5. We denote by Tjm the translation operator that maps the outgoing coefficients β⃗m from particle m to the local expansion α⃗l centered at particle l. With this operator in place, the incoming coefficients α⃗m for the m-th particle is

αm=am+j=1jmMTjmβj,
where a⃗m is the (truncated) local expansion (16) of the incident wave uinc on particle m. Tjm is referred to as the multipole-to-local (M2L) translation operator [12].

Combining Eqs. (21) and (24), one can easily eliminate the incoming coefficients α⃗m and obtain the following linear system that only involves the outgoing coefficients:

(𝒮1𝒯)[β1β2βM]=[a1a2aM],
where
𝒮=[SpSpSp],𝒯=[0T21TM1T120TM2T1MT2M0].

The system (25) can be solved iteratively, using GMRES [13]. Since each translation operator Tnm is dense, a naive matrix-vector product requires O((M(2p + 1))2) operations, where p is the order of the truncated expansion. FMM acceleration reduces the cost to O(M(2p + 1)2) work, for which we refer the reader to [12, 14]. Further, (25) has a simple diagonal preconditioner. Multiplying through by the block diagonal matrix 𝒮, results in the preconditioned system matrix I𝒮𝒯. This significantly reduces the number of iterations.

We now extend the multiple scattering approach to arbitrarily shaped particles.

Remark 6. It is worth emphasizing that the multiple scattering theory as discussed here is hardly new. We refer the reader to [1517] and the references therein.

4. Wave scattering for arbitrarily shaped particles

When the dielectric inclusions are of arbitrary shape, multiple scattering theory cannot be used quite so easily. Suppose, however, that an inclusion Ω is compactly supported with boundary Ω and that it is composed of a homogeneous material with wavenumber kp, as above. Given the incident wave uinc and the boundary conditions (8), (9), the exterior scattered field us and the field u within Ω have the following representation [5]:

us=Sk2σ+Dk2μ,forxΩc,
u=Skpσ+Dkpμ,forxΩ,
where Sk and Dk are the usual single layer and double layer potentials on Ω,
Skσ=ΩGk(x,y)σ(y)dsy,
Dkμ=ΩGk(x,y)n(y)μ(y)dsy.
σ(y) and μ(y) are unknown charge and dipole densities that lie on the boundary Ω. We will need the normal derivatives of Sk and Dk as well:
Nkσ=ΩGk(x,y)n(x)σ(y)dsy,Tkμ=Ω2Gk(x,y)n(x)n(y)μ(y)dsy.

By construction, the representations (26) and (27) satisfy the relevant Helmholtz equation in each domain. The single layer potential Sk is weakly singular and the value is well-defined for xΩ. The operators Dk and Nk are define on the boundary in the principal value sense (and have different limits when approaching the boundary from the interior and the exterior). The operator Tk is hypersingular with its value on the boundary defined in the Hadamard finite part sense. For further details, we refer the reader to [5].

Enforcing the interface conditions (8), (9) and taking appropriate limits [5] yields the following system of Fredholm integral equations of the second kind:

μ+[Sk2Skp]σ+[Dk2Dkp]μ=uinc,
σ+[Nk2Nkp]σ+[Tk2Tkp]μ=uincn.
Remark 7. It is worth noting that, while Tk is hypersingular, the difference kernel Tk2Tkp is only logarithmically singular and compact as are all the other difference operators in (32), at least for smooth boundaries. We use Nyström discretization for the system of equations based on the high order hybrid Gauss-trapezoidal rule of Alpert [18]. In this paper, we restrict our attention to smooth inclusions that are about one wavelength in size, so that 12 digits of accuracy are easily achieved with modest values of N using the Gauss-trapezoidal rule for logarithmic singularities of order 16. We refer the reader to [8] and the references therein for further details.

The integral Eq. (32) was introduced in electromagnetics by Müller [19], and in the scalar case by Kress, Rokhlin, Haider, Shipman and Venakides [11, 20, 21].

4.1. The scattering matrix

Suppose now that we have M inclusions Ω1,..., ΩM that are identical up to rotation, and well separated in the sense that each inclusion Ωi lies within a disk Di of radius R so that the disks are not overlapping. (see Fig. 2).

 

Fig. 2 Two inclusions and their enclosing disks. The scattering matrix Si for each inclusion Ωi with wavenumber kp is defined as the map from an incoming field on Di to the corresponding outgoing field. It is computed by solving a sequence of boundary value problems on the inclusion itself in a precomputation phase (see text). In this paper, we assume that all the inclusions are identical but may be rotated, as drawn here.

Download Full Size | PPT Slide | PDF

In that case, we can sample the incoming field on the disk Dj rather than Ωj as

u=n=ppαnJn(k2r)einθ,
using a polar coordinate system centered on the disk Dj.

Let σn and μn denote the solution to the integral Eq. (32) with right-hand side uinc = Jn(kr)einθ, uincn=kJn(kr)einθ. We may then precompute the multipole expansion from these source distributions

u=l=ppβlnHl(k2r)eilθ,
where
βln=Ωj[Jl(k2|y|)eilθj(y)σn(y)]+n[Jl(w|y|)eilθj(y)μn(y)]dsy.
Here, y is the location of a point on Ωj with respect to the center of disk Dj and θj(y) is the polar angle subtended with respect to the center of disk Dj. The formula for βl is standard [12, 14] and derived from the Graf addition theorem [10].

Definition 4.1. As before, the mapping between the incoming coefficients {αn} and outgoing coefficients {βn} is referred as the scattering matrix for the inclusion Ωj and denoted by Sj.

Remark 8. There exist several other methods for evaluating the scattering matrix Sj. A popular approach is the Extended Boundary Condition Method (EBCM) [22]. The EBCM does not involve singular integrals, which makes the implementation easier. However, it can suffer from numerical instabilities when the particle has a large aspect ratio. We have chosen to use the integral equation approach described above because it results in a well conditioned linear system for any geometry.

The reason for permitting a different scattering matrix for each inclusion is that the Ωj may be distinct in terms of geometry or dielectric properties. For the sake of simplicity, we assume here that the wavenumbers are the same in each inclusion and that the shapes are the same up to rotation. This permits us to solve only 2p + 1 integral equations on a single prototype inclusion in the enclosing disk. The scattering matrix for each rotated copy is then trivial to construct. Moreover, we can easily store the densities σn and μn, since this requires only O(2N(2p + 1)) storage, where N is the number of points used to discretize the boundary Ω. The amount of memory required to store the scattering matrix is O((2p + 1)2). For modest values of N, as is the case in the present paper, we compute the LU factors of the integral equation system matrix corresponding to (31), (32) only once, at a cost of O(N3) work. Each right-hand side corresponding to uinc = Jn(k2R)einθ and uincn=k2Jn(k2R)einθ can then be solved for n = −p,..., p at a total cost of O(N2(2p + 1)) work.

4.2. Multiple scattering

If we were interested in solving the multiple scattering problem in an infinite medium, we could now proceed as in the previous section. The number of degrees of freedom is only 2p + 1 per inclusion rather than N points per inclusions (the number needed to discretize the domain boundaries Ωj). For complicated inclusions, this permits a vast reduction in the number of degrees of freedom required and forms the basis for the FMPS method [4]. Moreover, the block-diagonal preconditioned multiple scattering equations are much better conditioned than the integral Eqs. (31) and (32) itself and FMM acceleration is particularly fast in this setting.

Remark 9. Extending the method to more than one type of substructure is straightforward as long as the assumption that the enclosed circles are well separated still holds. The additional cost is the bookkeeping for different scattering matrices of these substructures.

5. Multi-particle scattering in a layered medium

To this point, we have discussed the layered medium and multiple scattering problem separately. For the full problem, we now assume that multiple inclusions have been placed in the middle of a three-layered medium. We assume that the inclusions are well separated, so that the multiple scattering formalism applies within the layer. Then, we may write

u1(x)=Gk1(x,x0)+u1su2(x)=u2t+u2b+j=1Mn=ppβnmHn(k2rm)einθmu3(x)=u3s
where u1 and u3 denote the fields in the top and bottom half spaces and u2 denotes the field in the central layer exterior to the scattering disks Dj. u1s, u2t, u2b, and u3s are the Sommerfeld integrals from Section 2. Once u2 is known, the field within the scattering disks and the inclusions themselves is easily obtained.

It remains to discuss the discretization of the Sommerfeld integral, and the setup of the global linear system for the unknowns σ1, σ2+, σ2, σ3, and {β⃗m, m = 1...,M}.

5.1. Evaluation of the Sommerfeld integral

Let us consider the function u2t defined by (6). Its computation is a standard problem in acoustic and electromagnetic scattering and often handled by contour deformation. It is typical to deform the integration contour by pushing it from the real line into the second and fourth quadrants of the complex λ-plane in order to avoid the square root singularities in the integrand. One option is to use a hyperbolic tangent contour [8], which yields spectral accuracy with the trapezoidal rule and is extremely efficient. In our numerical simulation, we have chosen to use the piecewise smooth contour shown in Fig. 3 instead. This is slightly less efficient, but will permit us to evaluate the Sommerfeld integral using the non-uniform FFT, as explained further below. The contour consists of three segments: Γ1, Γ2 and Γ3, where

{Γ1:tib,t(0,),Γ2:it,t[b,b],Γ3:t+ib,t(,0).
The branch cuts for the square root in the integrand are chosen to ensure that waves are decaying away from the interface (up at k and down at −k, whether k is real or complex, as shown in Fig. 3).

 

Fig. 3 The Sommerfeld contour in the complex λ plane: Each segment in the contour is discretized using Gauss-Legendre quadrature. The branch cut (shown in red) points upward from k and downward from −k.

Download Full Size | PPT Slide | PDF

We truncate Γ1 and Γ3 at a point tmax > 0, where the integrand of u2t has decayed to a user-specified tolerance. Fortunately, the decay in the integrand is exponential once λ exceeds k2. (The precise rate of decay depends on the distance from the interface of the scattering disks and the point source generating the incoming field.) We let NS denote the number of points used in the quadrature for the Sommerfeld contour and note that each discretization point λj on the contour corresponds to a plane wave. We use the same contour and the same NS values {λj} for each of u1s, u2t, u2b, and u3s.

5.2. The full linear system

Let us denote by σ⃗ the discretized densities on the dielectric layers, σ=[σ1,σ2+,σ2,σ3]T, and by β⃗ the multipole coefficients for all M particles in the central layer. Each of σ⃗1, σ2+, σ2, and σ⃗3 is of length NS and the full linear system for multiple scattering in the layered medium takes the form of a block 2 × 2 linear system:

[ABCD][σβ]=[b0].

A itself is block diagonal 4NS × 4NS matrix with 4 × 4 blocks of the form Aλ in (10), each such block corresponding to a distinct λj in the contour integral discretization. The right-hand side component b is simply the right-hand side of (10) for each such λj. The matrix D = 𝒮𝒯I is simply the multiple scattering system for the particles from (25). The off-diagonal blocks B and C are more complicated. B is a matrix that translates the multipole expansion coefficients to a Sommerfeld representation on the upper and lower interfaces of the layered medium, while C requires the evaluation of the Sommerfield integral contributions from the interfaces in terms of incoming local expansions on the scattering disks themselves. We turn now to the efficient application of the matrices B and C.

5.2.1. The Sommerfeld-to-local operator

A straightforward mechanism to map from the σ⃗ variables to local expansions on the M disks is to use the Jacobi-Anger formula [10].

Lemma 5.1. Given r ∈ ℝ, k ∈ ℂ, we have

eikrcosθ=n=inJn(kr)einθ.

Suppose now that we want to compute the contribution from σ2+ to a local expansion on a disk centered at (x1, y1). Using Lemma 5.1, it is easy to see that

eλj2k22y+iλj(xx0)=eλj2k22y1+iλj(x1x0)n=inJn(k2r)ein(ϕ+θ),
where ϕ = arccos(λj/k2), θ = arccos((xx1)/r) and r=(xx1)2+(yy1)2. The analogous formula can be obtained for the contribution from σ2.

The cost of using formula (40) to compute the action of the C block in the system matrix above is clearly O(MNS(2p + 1)), where M denotes the number of particles and NS the number of discretization points λj in the Sommmerfeld contour and p is the order of the expansions used in the multiple scattering representation. This is quite acceptable when either NS or M is small. For high frequency problems with many inclusions, where k2 is large and NS = O(k2), we have developed a more efficient scheme, based on the nonuniform FFT (NUFFT).

5.2.2. The Sommerfeld-to-local operator using the NUFFT

Instead of mapping the contribution from the Sommerfeld integral to each disk separately, we seek a fast algorithm for evaluating the integral on a grid of points in the central layer, after which we can use high order interpolation to get the desired local expansion.

Restricting our attention to u2t for a fixed value of y, we have

14πΓ1eλ2k22yλ2k22eiλ(xx0)σ2+(λ)dλ=14πeb(xx0)0tmaxg(t)eitxdt,
where
g(t)=e(tib)2k22y(tib)2k22eitx0σ2+(tib).

Note now that the integral on the right-hand side of (41) is a finite Fourier transform. If we could compute it rapidly, we would have an efficient method for evaluating the Sommerfeld integral at a fine grid in the x variable for a fixed y. The discretization points in t, however, lie at Gauss-Legendre nodes, so the FFT itself does not apply. Fortunately, the nonuniform FFT (NUFFT) of Dutt and Rokhlin [23, 24] permits this to be done in nearly linear time. In our numerical simulations, we use the version discussed in [25,26]. The analogous method permits the rapid evaluation of the Sommerfeld integral on the contour Γ3. For the integral on Γ2, the NUFFT cannot be applied, but only a few discretization points are required, so we evaluate that contribution directly.

To provide rapid access to the field induced by the Sommerfeld integral at any location in the central layer, we superimpose on it a grid of n1 × n2 boxes that contain all of the M scattering disks. In each such box, we construct a tensor product m1 × m2 Chebyshev mesh, which will permit qth order local interpolation by barycentric interpolation [27]. The cost for evaluation at all grid points is O((n2m2)(n1m1 + NS) log(n1m1 + NS)) operations, using the NUFFT for each of the distinct n2m2 locations in y.

Consider now one of the scattering disks Dj of radius R. If we discretize the boundary of the disk using 2p + 1 equispaced points, evaluation of the induced field at each of the points requires O(m1m2) operations, for a net cost of O(m1m2(2p + 1)) work. An FFT of order (2p + 1) converts these values into their Fourier transforms, which we denote by {un}, for n = −p,···, p. From this, the n-th term an in the incoming J-expansion is simply

an=unJn(k2R).
Remark 10. The formula (42) will fail if the value k2R is a zero of the function Jn for any n from −p ..., p. This can be avoided if we also compute the normal derivative of the Sommerfeld integral on the boundary of each scattering disk. If we denote by {u′n} the Fourier coefficient of the normal derivative, it is easy to see that
an=unJn(k2R)+unkJn(k2R)Jn2(k2R)+(kJn(k2R))2,forj=p,,p.
The evaluation of the gradient of the Sommerfeld integral can be computed by an obvious modification of the formula (41) or (with a reduction in order) by computing the gradient of the tensor product Chebyshev series discussed above.

In summary, it requires O(Mm1m2(2p + 1)) operations to interpolate the field values on each of the M scattering disks and O(M(2p + 1)log(2p + 1)) operations to obtain the coefficients of the J-expansions. This completes the computation of the C block in the system matrix.

5.3. The multipole-to-Sommerfeld operator

The off-diagonal B block in (38) requires a formula for recasting the multipole expansion to the corresponding Sommerfeld representation on either the upper or lower interface of the layered medium. More precisely, each H-expansion in the central layer, centered on disk Dj with center (xj, yj) has a spectral representation on the upper layer y = 0 and the lower layer y = −d of the form:

ujt=14π1λ2k22eiλ(xx0)σmp+(λ)dλ,
ujb=14π1λ2k22eiλ(xx0)σmp(λ)dλ,
respectively.

The formulae for σmp+(λ) and σmp(λ) follow directly from the following theorem.

Theorem 5.2. [28] Let (xj, yj) denote the center of a multipole expansion in the central layer, withd < yj < 0 and let (r, θ) denote the polar coordinates of a target point with respect to that center. Then, on the upper interface,

Hn(kr)einθ=(1)n4πeλ2k2yjλ2k2eiλ(xxj)(λ2k2+k2k2)ndλ,
and on the lower interface,
Hn(kr)einθ=(1)n4πeλ2k2(d+yj)λ2k2eiλ(xxj)(λ2k2k2k2)ndλ.

Each multipole coefficient in the expansion about disk Dj contributes to each of the NS discretization points in the Sommerfeld integrals, requiring a total of O((2p + 1)NSM) work. This, then, is the cost of applying the B block of the system matrix directly.

5.3.1. The multipole-to-Sommerfeld operator using the NUFFT

Because of the computational complexity of applying the B block in the manner described above, it is important to develop a fast algorithm for the case where M and NS are large. We do so by essentially inverting the method of section 5.2.2. Assume first that all the centers of the H-expansions lie at the nodes of a uniform grid in the central layer and let us consider the contributions from the nth mode at each such grid point (xl, yj) for a fixed horizontal line y = yj. If there are n1 such expansion centers, with x coordinates {xl}, l = 1,···, n1, and we denote by anl the coefficient for the nth mode of the H-expansion at location (xl, yj), then the total contribution to the induced spectral coefficient σmp+(λj) on the top layer is given by

{σmp+(λj)}neλj2k22yj(λj2k22+k22k22)nl=1n1anleiλj(xjx0){σmp(λj)}neλj2k22(d+yj)(λj2k22k22k22)nl=1n1anleiλj(xjx0)

The formulae (48) imply that for each row, one can use the NUFFT to compute the induced coefficients for each discrete quadrature node λj on Γ1 or Γ3. As above, we use direct computation for the contributions to discretization nodes on Γ2. In the general case, the centers of the H-expansions are not aligned on a grid, but we can first shift the center of each H-expansion to the nearest grid point, using the multipole-to-multipole translation operator [12, 14] based on the Graf addition theorem [10]. After M such shifts, we may apply the transformation of (48).

The total computational cost is O(M(2p + 1)2) for shifting all the H-expansions and O(n2 (2p + 1)(n1 + NS) log(n1 + NS)) for the NUFFT-based work (see Table 1). The merits of the NUFFT-based schemes would become more apparent for larger NS.

Tables Icon

Table 1. Comparison of CPU time in seconds for the Sommerfeld-to-local and multipole-to-Sommerfeld operators, using both the direct and NUFFT-based schemes (see text). The Sommerfeld contour is discretized with 500 Gauss-Legendre points (240 points for Γ1 and Γ3, with 20 points for Γ2).

5.4. Iterative solution of the system matrix

We will solve Eq. (38) using the iterative method GMRES [13]. However, the unknowns σ⃗ and β⃗ may be poorly scaled with respect to each other. However, A is block diagonal, as noted above, with simple 4 × 4 blocks. Thus, we first invert A directly and use GMRES on the Schur complement of (38). In other words, we solve the system

[DCA1B][β]=CA1b
instead. This is much better conditioned and involves only the β⃗ unknowns. The Schur complement formalism has a simple physical interpretation: it is, in essence, a reformulation of the scattering problem using the layered medium Green’s function.

6. Numerical experiments

In this section, we illustrate the performance of our algorithm with three examples. For simplicity, we use a single class of inclusions, parametrized by

{x=(a1+a2cos(a3t))cos(t),y=(a1+a2cos(a3t))sin(t),for0t<2π.
As discussed in section 4, inclusions with more complicated boundaries do not introduce any essential difficulty in our scheme except that the precomputation of the scattering matrix is a little more involved, particulalry if corners are present [29, 30].

Given a fixed a1, a2 and a3, multiple copies of the inclusion are randomly distributed in the central layer of the medium with random orientations. To ensure the inclusions are well separated but confined in a fixed region, we use a bin sorting algorithm to construct the random distribution. We begin with inclusions located on a regular grid and then perturb their positions randomly, accepting the random move if the inclusion remains inside the region and stays well separated from the others. Several such sweeps are carried out to randomize the positions further.

We have not, as yet, specified the parameter NS used to discretize the Sommerfeld integral in (37). While special techniques have been developed by many authors to handle sources near the interface (see [7, 8, 31]), we simply assume that the source defining the incoming field is at least 0.2 wavelengths from the top interface. More precisely, in our examples, the source point in the top layer is placed at (1, 1) (which is roughly 0.2 wavelengths away for wavenumber k1 = 1). We also assume that the nearest inclusions get to either one of the interfaces in the layered medium is at least 0.5 wavelengths. Under these assumptions, we let tmax = max{|k1|, |k2|, |k3|} + 20, b = 0.2, and discretize Γ1 and Γ3 using 240 Gauss-Legendre points and Γ2 by 20 Gauss-Legendre points. This is sufficient to achieve about 10 digits of accuracy.

All computations are carried out using a 2.3GHz Intel Core i5 laptop, with 4GB RAM.

6.1. Example 1: scattering from large numbers of inclusions

In our first example, we consider the scattering of inclusions defined by parameters a1 = 0.12, a2 = 0.04, and a3 = 3 in Eq. (50) with wavenumber kp = 2.0. To obtain the scattering matrix with p = 10, we solve the integral Eqs. (31) and (32) by discretizing the boundary of the particle using N = 300 equispaced points. We assume the wavenumbers of the layered medium are given by k1 = 1.0, k2 = 3.0, k3 = 1.0. The thickness of the central layer is determined by the parameter d = 32. We consider distributions of M = 100, 500, 1, 000, 5, 000 inclusions and solve the mulitple scattering problem using GMRES with FMM acceleration. We terminate the iteration once the residual is less than 10−6. Results are presented in Fig. 4 and 5.

 

Fig. 4 Real part of the total field with 5, 000 dielectric inclusions randomly distributed in a three-layered medium. The wavenumber for each particle is kp = 2.0 and the wavenumbers for the three layers are k1 = 1.0, k2 = 3.0, k3 = 1.0. The size of each particle is approximately 0.1 wavelength for the wavenumber k2.

Download Full Size | PPT Slide | PDF

 

Fig. 5 Convergence behavior of GMRES and the CPU time required for various numbers of inclusions embedded in either (a) free space or (b)a three layered medium. For (a), we set k1 = k2 = k3 = 3.0 and for (b), we set k1 = 1.0, k2 = 3.0, k3 = 1.0.

Download Full Size | PPT Slide | PDF

Figure 4 shows the total field in the case M = 5, 000. The field distortion due to the inclusions is apparent. It requires 127s to achieve 6 digits of accuracy. Figure 5 shows the convergence behavior of GMRES as the number of inclusions is increased as well as the total CPU time. Clearly, more iterations are required for larger numbers of particles. Nevertheless, the time scales roughly linearly with the number of particles. In Fig. 5(a), we study the convergence rate when the background is homogeneous, by setting the material parameters to be the same for the three layers (k1 = k2 = k3). As expected, convergence is more rapid than when the inclusions are embedded in a true layered medium, because of the multiple reflections from the interfaces themselves.

6.2. Example 2: scattering in high contrast materials

In our second example, we consider the same inclusion shape as above, with kp = 2.0 + 1.0i, but with higher contrast materials. We fix the number of particles to be 200 and the thickness of the middle layer to be d = 12. We allow the wavenumber in the middle layer to vary from k2 = 1.0 + 0.01i up to k2 = 20.0 + 0.01i. Results are shown in Figs. 6 and 7 for 6 digits of accuracy.

 

Fig. 6 Real part of the total field for 200 dielectric inclusions distributed in a three layer medium with wavenumbers k1 = 1.0, k2 = 10.0 + 0.01i, k3 = 1.0. For each inclusion, the wavenumber is kp = 2.0 + 1.0i. The inclusions are approximately 0.3 wavelength in size for the wavenumber k2.

Download Full Size | PPT Slide | PDF

 

Fig. 7 Convergence behavior of GMRES iteration and the CPU time required for 200 inclusions embedded in the central layer, where k2 is allowed to vary from 1.0 + 0.01i to 20.0 + 0.01i. In (a), we create a homogeneous background by setting k1 = k2 = k3, while in (b), k1 and k3 are fixed at 1.0, and k2 varies.

Download Full Size | PPT Slide | PDF

In Fig. 7, we compare the convergence behavior in an infinite medium (a) vs. a layer medium (b). Note that the convergence is slower at high constrast and that this effect is more pronounced in the layerd medium case, where the central layer involves strong scattering and reflection.

6.3. Example 3: scattering from smoothed five-pointed stars

In our last example, we consider the scattering from a different inclusion shape, setting a1 = 0.3, a2 = 0.1, a3 = 5 in Eq. (50) with kp = 2.0. The inclusions are smoothed pentagons, as shown in Fig. 8. We discretize the boundary of the inclusion using N = 300 equispaced points and solve Eq. (31) and (32) to obtain the scattering matrix with p = 10. We consider M = 100, 200, 500 and 1, 000 inclusions. For the three-layered medium, we set k1 = 1.0, k2 = 3.0, k3 = 2.0. Results are shown in Figs. 8 and 9.

 

Fig. 8 Real part of the total field when 1, 000 inclusions are embedded in a three-layered medium with k1 = 1.0, k2 = 3.0, k3 = 2.0. Each inclusion is a smoothed five-pointed star, approximately 0.2 wavelengths in size.

Download Full Size | PPT Slide | PDF

 

Fig. 9 Convergence behavior of GMRES for a tolerance of 10−6 and the CPU time required as the number of inclusions embedded in the central layer varies. In (a), we create a homogeneous background by setting k1 = k2 = k3 = 3.0, while in (b), k1 = 1.0, k2 = 3.0, k3 = 2.0

Download Full Size | PPT Slide | PDF

Note that in order to obtain 6 digits of accuracy, 69.6 secs. are required for 1, 000 inclusions in a homogeneous background, while 140 seconds are required for the three layered medium. The time for convergence increases more or less in proportion to M2.

7. Conclusions

We have developed a fast algorithm to simulate electromagnetic scattering from a microstructured, three-layered material. Our methodology permits inclusions of arbitrary shape using a scattering matrix formalism combined with the use of Sommerfeld integrals to account for the influence of the layered material. We have designed efficient procedures to evaluate the Sommerfeld integral at arbitrary locations in the layered material using the non-uniform FFT and an effective preconditioner that allows the multiple scattering problem to be solved using GMRES with a modest number of iterations. As one would expect from physical considerations, the performance of the method degrades when the packing of inclusions is dense and when the contrast is high. While the method is suitable for parallel implementation, we are also investigating the possibility of replacing GMRES iteration with a fast direct solver [32].

Extension of the present method to the quasi-periodic case, where the incoming field impinges on a periodic microstructure will be reported at a later date.

Acknowledgments

This work was supported by the U.S. Department of Energy Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DEFGO288ER25053 and by the Office of the Assistant Secretary of Defense for Research and Engineering and AFOSR under NSSEFF Program Award FA9550-10-1-0180.

References and links

1. G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).

2. W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010). [CrossRef]  

3. Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009). [CrossRef]  

4. Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013). [CrossRef]  

5. D. Colton and R. Kress, Integral Equation Method in Scattering Theory (Wiley-Interscience, 1983).

6. W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995).

7. M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014). [CrossRef]  

8. A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011). [CrossRef]  

9. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93 (Springer-Verlag, 1998). [CrossRef]  

10. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

11. R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978). [CrossRef]  

12. V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” J. Comput. Phys. 86, 414–439 (1990). [CrossRef]  

13. Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).

14. H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006). [CrossRef]  

15. L. L. Foldy, “The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers,” Phys. Rev. 67, 107–119 (1945). [CrossRef]  

16. N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007). [CrossRef]  

17. K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013). [CrossRef]   [PubMed]  

18. B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. 20, 1551–1584 (1999). [CrossRef]  

19. C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer-Verlag, 1969). [CrossRef]  

20. M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002). [CrossRef]  

21. V. Rokhlin, “Solution of acoustic scattering problems by means of second kind integral equations,” Wave Motion 5, 257–272 (1983). [CrossRef]  

22. A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

23. A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993). [CrossRef]  

24. A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995). [CrossRef]  

25. L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004). [CrossRef]  

26. J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005). [CrossRef]  

27. J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004). [CrossRef]  

28. H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006). [CrossRef]  

29. J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010). [CrossRef]  

30. J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008). [CrossRef]  

31. M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012). [CrossRef]  

32. K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014). [CrossRef]  

References

  • View by:
  • |
  • |
  • |

  1. G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).
  2. W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
    [Crossref]
  3. Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009).
    [Crossref]
  4. Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013).
    [Crossref]
  5. D. Colton and R. Kress, Integral Equation Method in Scattering Theory (Wiley-Interscience, 1983).
  6. W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995).
  7. M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
    [Crossref]
  8. A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011).
    [Crossref]
  9. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93 (Springer-Verlag, 1998).
    [Crossref]
  10. F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).
  11. R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978).
    [Crossref]
  12. V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” J. Comput. Phys. 86, 414–439 (1990).
    [Crossref]
  13. Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).
  14. H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
    [Crossref]
  15. L. L. Foldy, “The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers,” Phys. Rev. 67, 107–119 (1945).
    [Crossref]
  16. N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007).
    [Crossref]
  17. K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
    [Crossref] [PubMed]
  18. B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. 20, 1551–1584 (1999).
    [Crossref]
  19. C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer-Verlag, 1969).
    [Crossref]
  20. M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
    [Crossref]
  21. V. Rokhlin, “Solution of acoustic scattering problems by means of second kind integral equations,” Wave Motion 5, 257–272 (1983).
    [Crossref]
  22. A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).
  23. A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993).
    [Crossref]
  24. A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995).
    [Crossref]
  25. L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004).
    [Crossref]
  26. J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005).
    [Crossref]
  27. J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004).
    [Crossref]
  28. H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
    [Crossref]
  29. J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
    [Crossref]
  30. J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008).
    [Crossref]
  31. M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012).
    [Crossref]
  32. K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
    [Crossref]

2014 (3)

G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).

M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
[Crossref]

K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
[Crossref]

2013 (2)

Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013).
[Crossref]

K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
[Crossref] [PubMed]

2012 (1)

M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012).
[Crossref]

2011 (1)

A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011).
[Crossref]

2010 (2)

W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
[Crossref]

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
[Crossref]

2009 (1)

Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009).
[Crossref]

2008 (1)

J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008).
[Crossref]

2007 (1)

N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007).
[Crossref]

2006 (2)

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
[Crossref]

2005 (1)

J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005).
[Crossref]

2004 (2)

J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004).
[Crossref]

L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004).
[Crossref]

2002 (1)

M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
[Crossref]

1999 (1)

B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. 20, 1551–1584 (1999).
[Crossref]

1995 (1)

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995).
[Crossref]

1993 (1)

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993).
[Crossref]

1990 (1)

V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” J. Comput. Phys. 86, 414–439 (1990).
[Crossref]

1986 (1)

Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).

1983 (1)

V. Rokhlin, “Solution of acoustic scattering problems by means of second kind integral equations,” Wave Motion 5, 257–272 (1983).
[Crossref]

1978 (1)

R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978).
[Crossref]

1945 (1)

L. L. Foldy, “The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers,” Phys. Rev. 67, 107–119 (1945).
[Crossref]

Abrahams, I. D.

W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
[Crossref]

Alpert, B. K.

B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. 20, 1551–1584 (1999).
[Crossref]

Bao, G.

G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).

Barnett, A.

A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011).
[Crossref]

Berrut, J.-P.

J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004).
[Crossref]

Boisvert, R. F.

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

Brazier-Smith, P. R.

W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
[Crossref]

Bremer, J.

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
[Crossref]

Cai, W.

M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012).
[Crossref]

Cheng, H.

H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
[Crossref]

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Chew, W. C.

W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995).

Cho, M. H.

M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012).
[Crossref]

Clark, C. W.

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

Colton, D.

D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93 (Springer-Verlag, 1998).
[Crossref]

D. Colton and R. Kress, Integral Equation Method in Scattering Theory (Wiley-Interscience, 1983).

Crutchfield, W.

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Duraiswami, R.

N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007).
[Crossref]

Dutt, A.

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995).
[Crossref]

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993).
[Crossref]

Foldy, L. L.

L. L. Foldy, “The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers,” Phys. Rev. 67, 107–119 (1945).
[Crossref]

Gimbutas, Z.

Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013).
[Crossref]

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Greengard, K. L.

K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
[Crossref]

Greengard, L.

M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
[Crossref]

Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013).
[Crossref]

A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011).
[Crossref]

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005).
[Crossref]

L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004).
[Crossref]

Gumerov, N. A.

N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007).
[Crossref]

Haider, M.

M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
[Crossref]

Helsing, J.

J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008).
[Crossref]

Ho, L.

K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
[Crossref]

Huang, J.

H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
[Crossref]

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Huang, K.

K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
[Crossref] [PubMed]

Kress, R.

R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978).
[Crossref]

D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93 (Springer-Verlag, 1998).
[Crossref]

D. Colton and R. Kress, Integral Equation Method in Scattering Theory (Wiley-Interscience, 1983).

Lacis, A. A.

A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Lai, J.

G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).

Lee, J.

J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005).
[Crossref]

L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004).
[Crossref]

Lee, J.-Y.

K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
[Crossref]

Leiterman, T. J.

H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
[Crossref]

Li, P.

K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
[Crossref] [PubMed]

Lozier, D. W.

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

Mishchenko, M. I.

A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Müller, C.

C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer-Verlag, 1969).
[Crossref]

O’Neil, M.

M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
[Crossref]

Ojala, R.

J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008).
[Crossref]

Olver, F. W. J.

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

Parnell, W. J.

W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
[Crossref]

Pataki, A.

M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
[Crossref]

Roach, G. F.

R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978).
[Crossref]

Rokhlin, V.

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
[Crossref]

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995).
[Crossref]

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993).
[Crossref]

V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” J. Comput. Phys. 86, 414–439 (1990).
[Crossref]

V. Rokhlin, “Solution of acoustic scattering problems by means of second kind integral equations,” Wave Motion 5, 257–272 (1983).
[Crossref]

Saad, Y.

Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).

Sammis, I.

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
[Crossref]

Schultz, M.

Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).

Shipman, S.

M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
[Crossref]

Travis, L. D.

A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Trefethen, L. N.

J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004).
[Crossref]

Venakides, S.

M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
[Crossref]

Wu, Y.

Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009).
[Crossref]

Yarvin, N.

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Zhang, Z.-Q.

Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009).
[Crossref]

Zhao, H.

K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
[Crossref] [PubMed]

Zhao, J.

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

Appl. Comput. Harmon. Anal. (1)

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data, II,” Appl. Comput. Harmon. Anal. 2, 85–100 (1995).
[Crossref]

BIT Numer. Math. (1)

A. Barnett and L. Greengard, “A new integral representation for quasi-periodic scattering problems in two dimensions,” BIT Numer. Math. 51, 67–90 (2011).
[Crossref]

Commun. Comput. Phys. (1)

G. Bao and J. Lai, “Radar cross section reduction of a cavity in the ground plane,” Commun. Comput. Phys. 15, 895–910 (2014).

Contemp. Math. (1)

H. Cheng, W. Crutchfield, Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin, and J. Zhao, “Remarks on the implementation of the wideband fmm for the helmholtz equation in two dimensions,” Contemp. Math. 408, 99–110 (2006).
[Crossref]

J. Comput. Phy. (1)

J. Helsing and R. Ojala, “Corner singularities for elliptic problems: Integral equations, graded meshes, quadrature, and compressed inverse preconditioning,” J. Comput. Phy. 227, 8820–8840 (2008).
[Crossref]

J. Comput. Phys. (9)

M. H. Cho and W. Cai, “A parallel fast algorithm for computing the Helmholtz integral operator in 3-d layered media,” J. Comput. Phys. 231, 5910–5925 (2012).
[Crossref]

K. L. Greengard, L. Ho, and J.-Y. Lee, “A fast direct solver for scattering from periodic structures with multiple material interfaces in two dimensions,” J. Comput. Phys. 258, 738–751 (2014).
[Crossref]

H. Cheng, J. Huang, and T. J. Leiterman, “An adaptive fast solver for the modified helmholtz equation in two dimensions,” J. Comput. Phys. 211, 616–637 (2006).
[Crossref]

J. Bremer, V. Rokhlin, and I. Sammis, “Universal quadratures for boundary integral equations on two-dimensional domains with corners,” J. Comput. Phys. 229, 8259–8280 (2010).
[Crossref]

J. Lee and L. Greengard, “The type 3 nonuniform FFT and its applications,” J. Comput. Phys. 206, 1–5 (2005).
[Crossref]

V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” J. Comput. Phys. 86, 414–439 (1990).
[Crossref]

N. A. Gumerov and R. Duraiswami, “A scalar potential formulation and translation theory for the time-harmonic maxwell equations,” J. Comput. Phys. 225, 206–236 (2007).
[Crossref]

K. Huang, P. Li, and H. Zhao, “An efficient algorithm for the generalized Foldy–Lax formulation,” J. Comput. Phys. 234, 376–398 (2013).
[Crossref] [PubMed]

Z. Gimbutas and L. Greengard, “Fast multi-particle scattering: A hybrid solver for the Maxwell equations in microstructured materials,” J. Comput. Phys. 232, 22–32 (2013).
[Crossref]

J. Math. Phys. (1)

R. Kress and G. F. Roach, “Transmission problems for the helmholtz equation,” J. Math. Phys. 19, 1433–1437 (1978).
[Crossref]

Phys. Rev. (1)

L. L. Foldy, “The multiple scattering of waves. i. general theory of isotropic scattering by randomly distributed scatterers,” Phys. Rev. 67, 107–119 (1945).
[Crossref]

Phys. Rev. B (1)

Y. Wu and Z.-Q. Zhang, “Dispersion relations and their symmetry properties of electromagnetic and elastic metamaterials in two dimensions,” Phys. Rev. B 79, 195111 (2009).
[Crossref]

Quantum J. Mech. Appl. Math. (1)

W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith, “Effective properties of a composite half-space: Exploring the relationship between homogenization and multiple-scattering theories,” Quantum J. Mech. Appl. Math. 63, 145–175 (2010).
[Crossref]

SIAM J. Appl. Math. (1)

M. Haider, S. Shipman, and S. Venakides, “Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channel defects and resonances,” SIAM J. Appl. Math. 62, 2129–2148 (2002).
[Crossref]

SIAM J. Sci. Comput. (2)

A. Dutt and V. Rokhlin, “Fast fourier transforms for nonequispaced data,” SIAM J. Sci. Comput. 14, 1368–1393 (1993).
[Crossref]

B. K. Alpert, “Hybrid Gauss-trapezoidal quadrature rules,” SIAM J. Sci. Comput. 20, 1551–1584 (1999).
[Crossref]

SIAM J. Sci. Stat. Comput. (1)

Y. Saad and M. Schultz, “Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear-systems,” SIAM J. Sci. Stat. Comput. 7, 856–869 (1986).

SIAM Rev. (2)

J.-P. Berrut and L. N. Trefethen, “Barycentric lagrange interpolation,” SIAM Rev. 46, 501–517 (2004).
[Crossref]

L. Greengard and J. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM Rev. 46, 443–454 (2004).
[Crossref]

Wave Motion (2)

V. Rokhlin, “Solution of acoustic scattering problems by means of second kind integral equations,” Wave Motion 5, 257–272 (1983).
[Crossref]

M. O’Neil, L. Greengard, and A. Pataki, “On the efficient representation of the half-space impedance green’s function for the helmholtz equation,” Wave Motion 51, 1–13 (2014).
[Crossref]

Other (6)

D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Applied Mathematical Sciences 93 (Springer-Verlag, 1998).
[Crossref]

F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark, NIST Handbook of Mathematical Functions (Cambridge University, 2010).

D. Colton and R. Kress, Integral Equation Method in Scattering Theory (Wiley-Interscience, 1983).

W. C. Chew, Waves and Fields in Inhomogeneous Media (IEEE, 1995).

C. Müller, Foundations of the Mathematical Theory of Electromagnetic Waves (Springer-Verlag, 1969).
[Crossref]

A. A. Lacis, L. D. Travis, and M. I. Mishchenko, Scattering, Absorption, and Emission of Light by Small Particles (Cambridge University, 2002).

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (9)

Fig. 1
Fig. 1 Geometry of the three-layered medium, with a large number of dielectric particles embedded in the middle layer.
Fig. 2
Fig. 2 Two inclusions and their enclosing disks. The scattering matrix Si for each inclusion Ωi with wavenumber kp is defined as the map from an incoming field on Di to the corresponding outgoing field. It is computed by solving a sequence of boundary value problems on the inclusion itself in a precomputation phase (see text). In this paper, we assume that all the inclusions are identical but may be rotated, as drawn here.
Fig. 3
Fig. 3 The Sommerfeld contour in the complex λ plane: Each segment in the contour is discretized using Gauss-Legendre quadrature. The branch cut (shown in red) points upward from k and downward from −k.
Fig. 4
Fig. 4 Real part of the total field with 5, 000 dielectric inclusions randomly distributed in a three-layered medium. The wavenumber for each particle is kp = 2.0 and the wavenumbers for the three layers are k1 = 1.0, k2 = 3.0, k3 = 1.0. The size of each particle is approximately 0.1 wavelength for the wavenumber k2.
Fig. 5
Fig. 5 Convergence behavior of GMRES and the CPU time required for various numbers of inclusions embedded in either (a) free space or (b)a three layered medium. For (a), we set k1 = k2 = k3 = 3.0 and for (b), we set k1 = 1.0, k2 = 3.0, k3 = 1.0.
Fig. 6
Fig. 6 Real part of the total field for 200 dielectric inclusions distributed in a three layer medium with wavenumbers k1 = 1.0, k2 = 10.0 + 0.01i, k3 = 1.0. For each inclusion, the wavenumber is kp = 2.0 + 1.0i. The inclusions are approximately 0.3 wavelength in size for the wavenumber k2.
Fig. 7
Fig. 7 Convergence behavior of GMRES iteration and the CPU time required for 200 inclusions embedded in the central layer, where k2 is allowed to vary from 1.0 + 0.01i to 20.0 + 0.01i. In (a), we create a homogeneous background by setting k1 = k2 = k3, while in (b), k1 and k3 are fixed at 1.0, and k2 varies.
Fig. 8
Fig. 8 Real part of the total field when 1, 000 inclusions are embedded in a three-layered medium with k1 = 1.0, k2 = 3.0, k3 = 2.0. Each inclusion is a smoothed five-pointed star, approximately 0.2 wavelengths in size.
Fig. 9
Fig. 9 Convergence behavior of GMRES for a tolerance of 10−6 and the CPU time required as the number of inclusions embedded in the central layer varies. In (a), we create a homogeneous background by setting k1 = k2 = k3 = 3.0, while in (b), k1 = 1.0, k2 = 3.0, k3 = 2.0

Tables (1)

Tables Icon

Table 1 Comparison of CPU time in seconds for the Sommerfeld-to-local and multipole-to-Sommerfeld operators, using both the direct and NUFFT-based schemes (see text). The Sommerfeld contour is discretized with 500 Gauss-Legendre points (240 points for Γ1 and Γ3, with 20 points for Γ2).

Equations (53)

Equations on this page are rendered with MathJax. Learn more.

Δ u s + k 2 u s = 0 ,
lim r r ( u s n i k u s ) = 0 ,
G k 1 ( x , x 0 ) = 1 4 π e λ 2 k 1 2 | y y 0 | λ 2 k 1 2 e i λ ( x x 0 ) d λ .
u 1 s = 1 4 π e λ 2 k 1 2 y λ 2 k 1 2 e i λ ( x x 0 ) σ 1 ( λ ) d λ ,
u 2 t = 1 4 π e λ 2 k 2 2 y λ 2 k 2 2 e i λ ( x x 0 ) σ 2 + ( λ ) d λ ,
u 2 b = 1 4 π e λ 2 k 2 2 ( y + d ) λ 2 k 2 2 e i λ ( x x 0 ) σ 2 ( λ ) d λ ,
u 3 s = 1 4 π e λ 2 k 3 2 ( y + d ) λ 2 k 3 2 e i λ ( x x 0 ) σ 3 ( λ ) d λ .
[ u ] = 0 ,
[ u n ] = 0 ,
( 1 λ 2 k 1 2 1 λ 2 k 2 2 e λ 2 k 2 2 d λ 2 k 2 2 0 0 e λ 2 k 2 2 d λ 2 k 2 2 1 λ 2 k 2 2 1 λ 2 k 3 2 1 1 e λ 2 k 2 2 d 0 0 e λ 2 k 2 2 d 1 1 ) ( σ 1 ( λ ) σ 2 + ( λ ) σ 2 ( λ ) σ 3 ( λ ) ) = ( e λ 2 k 1 2 y 0 λ 2 k 1 2 0 e λ 2 k 1 2 y 0 0 )
Δ u + k 2 2 u = 0 .
lim r r ( u s n i k 2 u s ) = 0 ,
Δ u + k p 2 u = 0 .
u s = n = β n H n ( k 2 r ) e i n θ
u = n = γ n J n ( k p r ) e i n θ
u inc = n = α n J n ( k 2 r ) e in θ , u inc r = n = α n k 2 J n ( k 2 r ) e in θ .
[ H n ( k 2 R ) J n ( k p R ) k 2 H n ( k 2 R ) k p J n ( k p R ) ] [ β n γ n ] = [ α n J n ( k 2 R ) α n k 2 J n ( k 2 R ) ] ,
β n = [ k p J n ( k 2 R ) J n ( k p R ) k 2 J n ( k 2 R ) J n ( k p R ) k 2 H n ( k 2 R ) J n ( k R ) k p J n ( k p R ) H n ( k 2 R ) ] α n ,
γ n = [ k 2 J n ( k 2 R ) H n ( k 2 R ) k 2 J n ( k 2 R ) H n ( k 2 R ) k 2 H n ( k 2 R ) J n ( k p R ) k p J n ( k p R ) H n ( k 2 R ) ] α n .
β n = J n ( k 2 R ) H n ( k 2 R ) α n .
β m = S p [ α m ] , for m = 1 , , M .
n = β n m H n ( k 2 r m ) e in θ m
u = n = α n l J n ( k 2 r l ) e in θ l
α n l = n = e in ( θ l m π ) β n n m H n ( k 2 x m x l ) .
α m = a m + j = 1 j m M T j m β j ,
( 𝒮 1 𝒯 ) [ β 1 β 2 β M ] = [ a 1 a 2 a M ] ,
𝒮 = [ S p S p S p ] , 𝒯 = [ 0 T 21 T M 1 T 12 0 T M 2 T 1 M T 2 M 0 ] .
u s = S k 2 σ + D k 2 μ , for x Ω c ,
u = S k p σ + D k p μ , for x Ω ,
S k σ = Ω G k ( x , y ) σ ( y ) d s y ,
D k μ = Ω G k ( x , y ) n ( y ) μ ( y ) d s y .
N k σ = Ω G k ( x , y ) n ( x ) σ ( y ) d s y , T k μ = Ω 2 G k ( x , y ) n ( x ) n ( y ) μ ( y ) d s y .
μ + [ S k 2 S k p ] σ + [ D k 2 D k p ] μ = u inc ,
σ + [ N k 2 N k p ] σ + [ T k 2 T k p ] μ = u inc n .
u = n = p p α n J n ( k 2 r ) e in θ ,
u = l = p p β l n H l ( k 2 r ) e i l θ ,
β l n = Ω j [ J l ( k 2 | y | ) e i l θ j ( y ) σ n ( y ) ] + n [ J l ( w | y | ) e i l θ j ( y ) μ n ( y ) ] d s y .
u 1 ( x ) = G k 1 ( x , x 0 ) + u 1 s u 2 ( x ) = u 2 t + u 2 b + j = 1 M n = p p β n m H n ( k 2 r m ) e in θ m u 3 ( x ) = u 3 s
{ Γ 1 : t i b , t ( 0 , ) , Γ 2 : i t , t [ b , b ] , Γ 3 : t + i b , t ( , 0 ) .
[ A B C D ] [ σ β ] = [ b 0 ] .
e i k r cos θ = n = i n J n ( k r ) e in θ .
e λ j 2 k 2 2 y + i λ j ( x x 0 ) = e λ j 2 k 2 2 y 1 + i λ j ( x 1 x 0 ) n = i n J n ( k 2 r ) e in ( ϕ + θ ) ,
1 4 π Γ 1 e λ 2 k 2 2 y λ 2 k 2 2 e i λ ( x x 0 ) σ 2 + ( λ ) d λ = 1 4 π e b ( x x 0 ) 0 t max g ( t ) e i t x d t ,
g ( t ) = e ( t i b ) 2 k 2 2 y ( t i b ) 2 k 2 2 e i t x 0 σ 2 + ( t i b ) .
a n = u n J n ( k 2 R ) .
a n = u n J n ( k 2 R ) + u n k J n ( k 2 R ) J n 2 ( k 2 R ) + ( k J n ( k 2 R ) ) 2 , for j = p , , p .
u j t = 1 4 π 1 λ 2 k 2 2 e i λ ( x x 0 ) σ m p + ( λ ) d λ ,
u j b = 1 4 π 1 λ 2 k 2 2 e i λ ( x x 0 ) σ m p ( λ ) d λ ,
H n ( k r ) e in θ = ( 1 ) n 4 π e λ 2 k 2 y j λ 2 k 2 e i λ ( x x j ) ( λ 2 k 2 + k 2 k 2 ) n d λ ,
H n ( k r ) e in θ = ( 1 ) n 4 π e λ 2 k 2 ( d + y j ) λ 2 k 2 e i λ ( x x j ) ( λ 2 k 2 k 2 k 2 ) n d λ .
{ σ m p + ( λ j ) } n e λ j 2 k 2 2 y j ( λ j 2 k 2 2 + k 2 2 k 2 2 ) n l = 1 n 1 a n l e i λ j ( x j x 0 ) { σ m p ( λ j ) } n e λ j 2 k 2 2 ( d + y j ) ( λ j 2 k 2 2 k 2 2 k 2 2 ) n l = 1 n 1 a n l e i λ j ( x j x 0 )
[ D CA 1 B ] [ β ] = CA 1 b
{ x = ( a 1 + a 2 cos ( a 3 t ) ) cos ( t ) , y = ( a 1 + a 2 cos ( a 3 t ) ) sin ( t ) , for 0 t < 2 π .

Metrics