## 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 [1–3]. 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 *e*^{−iω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 *u ^{s}* must satisfy the equation

*k*

_{1},

*k*

_{2}and

*k*

_{3}the wavenumber for the three layers, and by

*k*the wavenumber for the particles. The scattered field also has to satisfy the Sommerfeld radiation condition at infinity [5]:

_{p}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 **x _{0}** = (

*x*

_{0},

*y*

_{0}) in the top layer, with wavenumber

*k*

_{1}, the corresponding field is given by the (two-dimensional) free space Green’s function: ${G}_{{k}_{1}}(\mathbf{x},{\mathbf{x}}_{\mathbf{0}})=\frac{i}{4}{H}_{0}^{1}\left({k}_{1}\left|\mathbf{x}-{\mathbf{x}}_{\mathbf{0}}\right|\right)$, where ${H}_{0}^{1}(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:

*y*≠

*y*

_{0}.

In the Sommerfeld approach [6], we assume the upward scattered field ${u}_{1}^{s}$ in the top layer can be expressed as

*σ*

_{1}(

*λ*) is an unknown density on the upper interface

*y*= 0. It is straghtforward to verify that ${u}_{1}^{s}$ satisfies the Helmholtz equation with Helmholtz parameter

*k*

_{1}.

In the second layer, the scattered field
${u}_{2}^{s}$ can be written in terms of contributions from both the upper (*y* = 0) and lower (*y* = −*d*) interfaces:
${u}_{2}^{t}$ and
${u}_{2}^{b}$. These are given by

*spectral density*functions on the upper and lower interfaces.

Similarly, we can represent the scattered field
${u}_{3}^{s}$ in the third layer with an unknown density *σ*_{3}(*λ*) on the lower interface as

*Remark*1. The signs of the terms ${e}^{\pm \sqrt{{\lambda}^{2}-{{k}_{i}}^{2}y}}$ and ${e}^{\pm \sqrt{{\lambda}^{2}-{{k}_{i}}^{2}}(y+d)}$ in Eqs. (4)–(7) ensure that evanescent modes (when |

*λ*| > |

*k*|) decay away from each layer. (Physically, this is related to causality and is required in the derivation of formula [7] by contour integration).

_{i}It is worth noting that the four unknown densities *σ*_{1},
${\sigma}_{2}^{+}$,
${\sigma}_{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},
${\sigma}_{2}^{+}$,
${\sigma}_{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

*∂/∂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:

**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
${k}_{p}=\omega \sqrt{{\epsilon}_{p}\mu}$ embedded in **R**^{2}, assumed to consist of a dielectric with
${k}_{2}=\omega \sqrt{{\epsilon}_{2}\mu}$. For transverse magnetic(TM) polarization, the total electrical field *u* in the exterior of the inclusion satisfies the Helmholtz equation:

*u*can be written as the sum of the incident field

*u*and the scattered field

^{inc}*u*, where

^{s}*u*satisfies (11) and the Sommerfeld radiation condition,

^{s}*u*satisfies the Helmholtz equation with wavenumber

*k*, On the boundary of the inclusion, we must enforce the continuity conditions given by Eqs. (8) and (9).

_{p}#### 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

*r*,

*θ*) are the polar coordinates of a point in the plane,

*H*(

_{n}*r*) is the Hankel function of the first kind of order

*n*and

*J*(

_{n}*r*) is the Bessel function of order

*n*[9, 10].

We now expand the incident wave *u ^{inc}* and its normal derivative in the form:

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*:

*n*∈ ℕ.

Solving Eq. (17) determines the coefficients *β _{n}*,

*γ*:

_{n}*k*and

*k*have positive real part and non-negative imaginary part [5, 11].

_{p}**Definition 3.1.** The mapping between the incoming coefficients {*α _{n}*} and outgoing coefficients {

*β*} is referred as the scattering matrix for the disk and denoted by

_{n}*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

*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 *k _{p}* and the background medium again has wavenumer

*k*

_{2}. For each individual particle, the analysis can be carried out as above. We will denote by

*α⃗*the incoming coefficients and by

^{m}*β⃗*the outgoing coefficients for the

^{m}*m*-th particle. We have

*S*denotes the truncated (2

_{p}*p*+ 1) × (2

*p*+ 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 *u ^{inc}* and the contribution to the scattered field

*u*from all of the other particles. In order to formulate the problem concisely, given the multipole expansion for disk

^{s}*j*, we need some additional notation.

**Lemma 3.1.** *[12] Let disk m be centered at* **x _{m}**

*and let disk l be centered at*

**x**

_{l}*. Then the multipole expansion*

*induces a field on disk l of the form*

*where*

*Here,*(

*r*,

_{m}*θ*)

_{m}*and*(

*r*,

_{l}*θ*)

_{l}*denote the polar coordinates of a target point with respect to disk centers*

**x**

_{m}*and*

**x**

_{l}*, respectively and*

*θ*

_{lm}*denotes the angle between*(

**x**−

_{m}**x**)

_{l}*and the x-axis.*

*Remark* 5. We denote by *T ^{jm}* the translation operator that maps the outgoing coefficients

*β⃗*from particle

^{m}*m*to the local expansion

*α⃗*centered at particle

^{l}*l*. With this operator in place, the incoming coefficients

*α⃗*for the

^{m}*m*-th particle is

*a⃗*is the (truncated) local expansion (16) of the incident wave

^{m}*u*on particle

^{inc}*m*.

*T*is referred to as the multipole-to-local (M2L) translation operator [12].

^{jm}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:

The system (25) can be solved iteratively, using GMRES [13]. Since each translation operator *T ^{nm}* is dense, a naive matrix-vector product requires

*O*((

*M*(2

*p*+ 1))

^{2}) operations, where

*p*is the order of the truncated expansion. FMM acceleration reduces the cost to

*O*(

*M*(2

*p*+ 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 [15–17] 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 *k _{p}*, as above. Given the incident wave

*u*and the boundary conditions (8), (9), the exterior scattered field

^{inc}*u*and the field

^{s}*u*within Ω have the following representation [5]:

**S**

*and*

^{k}**D**

*are the usual single layer and double layer potentials on*

^{k}*∂*Ω,

*σ*(

*y*) and

*μ*(

*y*) are unknown charge and dipole densities that lie on the boundary

*∂*Ω. We will need the normal derivatives of

**S**

*and*

^{k}**D**

*as well:*

^{k}By construction, the representations (26) and (27) satisfy the relevant Helmholtz equation in each domain. The single layer potential **S*** ^{k}* is weakly singular and the value is well-defined for

**x**∈

*∂*Ω. The operators

**D**

*and*

^{k}**N**

*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*

^{k}**T**

*is*

^{k}*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:

*Remark*7. It is worth noting that, while

**T**

*is hypersingular, the difference kernel*

^{k}**T**

^{k2}−

**T**

^{kp}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 Ω

*lies within a disk*

_{i}*D*of radius

_{i}*R*so that the disks are not overlapping. (see Fig. 2).

In that case, we can sample the incoming field on the disk *D _{j}* rather than Ω

*as*

_{j}*D*.

_{j}Let *σ _{n}* and

*μ*denote the solution to the integral Eq. (32) with right-hand side

_{n}*u*=

^{inc}*J*(

_{n}*kr*)

*e*, $\frac{\partial {u}^{\mathit{inc}}}{\partial n}=k{J}_{n}^{\prime}(kr){e}^{\mathit{in}\theta}$. We may then precompute the multipole expansion from these source distributions

^{inθ}**y**is the location of a point on

*∂*Ω

*with respect to the center of disk*

_{j}*D*and

_{j}*θ*(

_{j}**y**) is the polar angle subtended with respect to the center of disk

*D*. The formula for

_{j}*β*is standard [12, 14] and derived from the Graf addition theorem [10].

_{l}**Definition 4.1.** As before, the mapping between the incoming coefficients {*α _{n}*} and outgoing coefficients {

*β*} is referred as the scattering matrix for the inclusion Ω

_{n}*and denoted by*

_{j}*S*.

_{j}*Remark* 8. There exist several other methods for evaluating the scattering matrix *S _{j}*. 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 2

*p*+ 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

*σ*and

_{n}*μ*, since this requires only

_{n}*O*(2

*N*(2

*p*+ 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*((2

*p*+ 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*(

*N*

^{3}) work. Each right-hand side corresponding to

*u*=

^{inc}*J*(

_{n}*k*

_{2}

*R*)

*e*and $\frac{\partial {u}^{\mathit{inc}}}{\partial n}={k}_{2}{J}_{n}^{\prime}({k}_{2}R){e}^{\mathit{in}\theta}$ can then be solved for

^{inθ}*n*= −

*p*,...,

*p*at a total cost of

*O*(

*N*

^{2}(2

*p*+ 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 2*p* + 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

*u*

_{1}and

*u*

_{3}denote the fields in the top and bottom half spaces and

*u*

_{2}denotes the field in the central layer exterior to the scattering disks

*D*. ${u}_{1}^{s}$, ${u}_{2}^{t}$, ${u}_{2}^{b}$, and ${u}_{3}^{s}$ are the Sommerfeld integrals from Section 2. Once

_{j}*u*

_{2}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},
${\sigma}_{2}^{+}$,
${\sigma}_{2}^{-}$, *σ*_{3}, and {*β⃗ _{m}*,

*m*= 1...,

*M*}.

#### 5.1. Evaluation of the Sommerfeld integral

Let us consider the function
${u}_{2}^{t}$ 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

*k*and down at −

*k*, whether

*k*is real or complex, as shown in Fig. 3).

We truncate Γ_{1} and Γ_{3} at a point *t _{max}* > 0, where the integrand of
${u}_{2}^{t}$ has decayed to a user-specified tolerance. Fortunately, the decay in the integrand is exponential once

*λ*exceeds

*k*

_{2}. (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

*N*denote the number of points used in the quadrature for the Sommerfeld contour and note that each discretization point

_{S}*λ*on the contour corresponds to a plane wave. We use the same contour and the same

_{j}*N*values {

_{S}*λ*} for each of ${u}_{1}^{s}$, ${u}_{2}^{t}$, ${u}_{2}^{b}$, and ${u}_{3}^{s}$.

_{j}#### 5.2. The full linear system

Let us denote by *σ⃗* the discretized densities on the dielectric layers,
$\overrightarrow{\sigma}={\left[{\overrightarrow{\sigma}}_{1},{\overrightarrow{\sigma}}_{2}^{+},{\overrightarrow{\sigma}}_{2}^{-},{\overrightarrow{\sigma}}_{3}\right]}^{T}$, and by *β⃗* the multipole coefficients for all *M* particles in the central layer. Each of *σ⃗*_{1},
${\overrightarrow{\sigma}}_{2}^{+}$,
${\overrightarrow{\sigma}}_{2}^{-}$, and *σ⃗*_{3} is of length *N _{S}* and the full linear system for multiple scattering in the layered medium takes the form of a block 2 × 2 linear system:

*A* itself is block diagonal 4*N _{S}* × 4

*N*matrix with 4 × 4 blocks of the form

_{S}*A*in (10), each such block corresponding to a distinct

_{λ}*λ*in the contour integral discretization. The right-hand side component

_{j}*b*is simply the right-hand side of (10) for each such

*λ*. The matrix

_{j}*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*

Suppose now that we want to compute the contribution from
${\sigma}_{2}^{+}$ to a local expansion on a disk centered at (*x*_{1}, *y*_{1}). Using Lemma 5.1, it is easy to see that

*ϕ*= arccos(

*λ*

_{j}/k_{2}),

*θ*= arccos((

*x*−

*x*

_{1})/

*r*) and $r=\sqrt{{(x-{x}_{1})}^{2}+{(y-{y}_{1})}^{2}}$. The analogous formula can be obtained for the contribution from ${\sigma}_{2}^{-}$.

The cost of using formula (40) to compute the action of the *C* block in the system matrix above is clearly *O*(*MN _{S}*(2

*p*+ 1)), where

*M*denotes the number of particles and

*N*the number of discretization points

_{S}*λ*in the Sommmerfeld contour and

_{j}*p*is the order of the expansions used in the multiple scattering representation. This is quite acceptable when either

*N*or

_{S}*M*is small. For high frequency problems with many inclusions, where

*k*

_{2}is large and

*N*=

_{S}*O*(

*k*

_{2}), 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
${u}_{2}^{t}$ for a fixed value of *y*, we have

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 *n*_{1} × *n*_{2} boxes that contain all of the *M* scattering disks. In each such box, we construct a tensor product *m*_{1} × *m*_{2} Chebyshev mesh, which will permit *q*th order local interpolation by barycentric interpolation [27]. The cost for evaluation at all grid points is *O*((*n*_{2}*m*_{2})(*n*_{1}*m*_{1} + *N _{S}*) log(

*n*

_{1}

*m*

_{1}+

*N*)) operations, using the NUFFT for each of the distinct

_{S}*n*

_{2}

*m*

_{2}locations in

*y*.

Consider now one of the scattering disks *D _{j}* of radius

*R*. If we discretize the boundary of the disk using 2

*p*+ 1 equispaced points, evaluation of the induced field at each of the points requires

*O*(

*m*

_{1}

*m*

_{2}) operations, for a net cost of

*O*(

*m*

_{1}

*m*

_{2}(2

*p*+ 1)) work. An FFT of order (2

*p*+ 1) converts these values into their Fourier transforms, which we denote by {

*u*}, for

_{n}*n*= −

*p*,···,

*p*. From this, the

*n*-th term

*a*in the incoming

_{n}*J*-expansion is simply

*Remark*10. The formula (42) will fail if the value

*k*

_{2}

*R*is a zero of the function

*J*for any

_{n}*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′*} the Fourier coefficient of the normal derivative, it is easy to see that

_{n}In summary, it requires *O*(*Mm*_{1}*m*_{2}(2*p* + 1)) operations to interpolate the field values on each of the *M* scattering disks and *O*(*M*(2*p* + 1)log(2*p* + 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 *D _{j}* with center (

*x*,

_{j}*y*) has a spectral representation on the upper layer

_{j}*y*= 0 and the lower layer

*y*= −

*d*of the form:

The formulae for ${\sigma}_{mp}^{+}(\lambda )$ and ${\sigma}_{mp}^{-}(\lambda )$ follow directly from the following theorem.

**Theorem 5.2.** *[28] Let* (*x _{j}*,

*y*)

_{j}*denote the center of a multipole expansion in the central layer, with*−

*d*<

*y*< 0

_{j}*and let*(

*r*,

*θ*)

*denote the polar coordinates of a target point with respect to that center. Then, on the upper interface,*

*and on the lower interface,*

Each multipole coefficient in the expansion about disk *D _{j}* contributes to each of the

*N*discretization points in the Sommerfeld integrals, requiring a total of

_{S}*O*((2

*p*+ 1)

*N*) work. This, then, is the cost of applying the

_{S}M*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 *N _{S}* 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

*n*th mode at each such grid point (

*x*,

_{l}*y*) for a fixed horizontal line

_{j}*y*=

*y*. If there are

_{j}*n*

_{1}such expansion centers, with

*x*coordinates {

*x*},

_{l}*l*= 1,···,

*n*

_{1}, and we denote by ${a}_{n}^{l}$ the coefficient for the

*n*th mode of the

*H*-expansion at location (

*x*,

_{l}*y*), then the total contribution to the induced spectral coefficient ${\sigma}_{mp}^{+}({\lambda}_{j})$ on the top layer is given by

_{j}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*(2*p*
+ 1)^{2}) for shifting all the *H*-expansions
and *O*(*n*_{2} (2*p*
+ 1)(*n*_{1} +
*N _{S}*) log(

*n*

_{1}+

*N*)) for the NUFFT-based work (see Table 1). The merits of the NUFFT-based schemes would become more apparent for larger

_{S}*N*.

_{S}#### 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

*β⃗*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

Given a fixed *a*_{1}, *a*_{2} and *a*_{3}, 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 *N _{S}* 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

*k*

_{1}= 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

*t*= max{|

_{max}*k*

_{1}|, |

*k*

_{2}|, |

*k*

_{3}|} + 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 *a*_{1} = 0.12, *a*_{2} = 0.04, and *a*_{3} = 3 in Eq. (50) with wavenumber *k _{p}* = 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

*k*

_{1}= 1.0,

*k*

_{2}= 3.0,

*k*

_{3}= 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.

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 (*k*_{1} = *k*_{2} = *k*_{3}). 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 *k _{p}* = 2.0 + 1.0

*i*, 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

*k*

_{2}= 1.0 + 0.01

*i*up to

*k*

_{2}= 20.0 + 0.01

*i*. Results are shown in Figs. 6 and 7 for 6 digits of accuracy.

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 *a*_{1} = 0.3, *a*_{2} = 0.1, *a*_{3} = 5 in Eq. (50) with *k _{p}* = 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

*k*

_{1}= 1.0,

*k*

_{2}= 3.0,

*k*

_{3}= 2.0. Results are shown in Figs. 8 and 9.

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 *M*^{2}.

## 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]