A fast solver for multi-particle scattering in a layered medium
AA fast solver for multi-particle scattering in a layered medium
Jun Lai ∗ Motoki Kobayashi ∗ Leslie Greengard ∗† Abstract
In this paper, we consider acoustic or electromagnetic scattering in two dimensionsfrom an infinite three-layer medium with thousands of wavelength-size dielectric particlesembedded in the middle layer. Such geometries are typical of microstructured compositematerials, and the evaluation of the scattered field requires a suitable fast solver for eithera single configuration or for a sequence of configurations as part of a design or optimizationprocess. We have developed an algorithm for problems of this type by combining theSommerfeld integral representation, high order integral equation discretization, the fastmultipole method and classical multiple scattering theory. The efficiency of the solver isillustrated with several numerical experiments.
Keywords : Helmholtz equation, multiple scattering, layered medium, Sommerfeld inte-gral, composite material design
The problem of designing composite materials that exhibit a specific acoustic or electro-magnetic response is an area of active research [2, 26, 30]. Examples include the design ofrandom media with a well-defined macroscopic refraction (coherent scattering) [26] and thefabrication of metamaterials [30] 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 ofthe incoming field and the distribution of particles is reasonably dense, then the interactionof the particles involves non-negligible multiple scattering effects and methods based onhomogenization [26] are not applicable. Instead, the full Helmholtz or Maxwell equationsshould be solved at each iteration of the design process. Numerical simulation, in theabsence of suitable fast algorithms, are impractical when thousands of particles are involved.In this paper, we develop an algorithm that accelerates the computation of electromag-netic 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 minutesto evaluate the scattered field for up to 5 ,
000 particles on a 2 . z direction. Arelated three-dimensional solver was considered in [15], but the particles were assumed tobe distributed in free space. A principal contribution of this paper is the development ofa mathematical framework that permits them to be embedded in a layer material (whichis closer to being manufacturable). While we restrict our attention here to the three-layercase, the extension to an arbitrary layered medium is straightforward. ∗ Courant Institute of Mathematical Sciences, New York University, NY 10012 (Email: [email protected],[email protected], [email protected]) † Simons Center for Data Analysis, Simons Foundation, New York, NY 10010 a r X i v : . [ m a t h . NA ] J u l igure 1: Geometry of the three-layered medium, with a large number of dielectric particlesembedded in the middle layer. More precisely, we consider time-harmonic scattering (with time dependence e iωt ) froma three-layered medium as depicted in Fig. 1. The incident field is assumed to be drivenby a point source located in the first (top) layer. The thickness of the middle layer isdenoted by d . We assume the magnetic permeability µ is identical in each layer, while theelectric permittivity (cid:15) is piecewise constant. There are two fundamental polarizations inthe two dimensional setting to consider: the transverse magnetic (TM) polarization andthe transverse electric (TE) polarization. In both cases, the Maxwell equations reduce to ascalar Helmholtz equation. For simplicity, we consider the TM polarization here, in whichcase the scattered field u s must satisfy the equation∆ u s + k u s = 0 , (1)where k = ω √ µ(cid:15) is the wavenumber. We denote by k , k and k the wavenumber for thethree layers, and by k p the wavenumber for the particles. The scattered field also has tosatisfy the Sommerfeld radiation condition at infinity [10]:lim r →∞ √ r (cid:18) ∂u s ∂n − iku s (cid:19) = 0 , (2)where r = (cid:112) x + y .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 [15], we assume that theparticles are well separated from each other - that is, the separation between particles isat least 10% of the particle size. Second, we assume that only a finite number of distinctparticle shapes are included in the simulation. The first condition ensures that a multiplescattering formalism will be accurate and the second condition ensures that precomputationof single particle scattering matrices permits a dramatic reduction in the number of degreesof freedom necessary for the solver. The particles are not assumed to be symmetric andmay 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 andits application to layered materials (in the absence of inclusions). In Section 3, we reviewclassical multiple scattering theory for circular particles. Section 4 extends the scatteringformalism to non-circular particles and Section 5 develops analytical tools needed to goback and forth between the Sommerfeld integral formalism and multiple scattering theory.Section 5 also combines the techniques in the preceding sections and extends the FMPSmethod to layered media. Numerical examples are provided in Section 6 to illustrate theefficiency of the method, followed by some concluding remarks in Section 7. The Sommerfeld integral for layered media
Wave propagation in a layered medium is a well-studied problem in acoustic and elec-tromagnetic scattering theory. Nearly a century ago, Sommerfeld developed a spectralrepresentation involving a Fourier integral in the “transverse” variable (the x -coordinatein Fig. 1) [8]. Assuming a point source is located at x = ( x , y ) in the top layer, withwavenumber k , the corresponding field is given by the (two-dimensional) free space Green’sfunction: G k ( x , x ) = i H (1)0 ( k | x − x | ), where H (1)0 ( x ) is the first kind Hankel functionof order zero. Combing the Fourier transform and contour integration [25], the Green’sfunction can also be written in the form: G k ( x , x ) = 14 π (cid:90) ∞−∞ e − √ λ − k | y − y | (cid:112) λ − k e iλ ( x − x ) dλ. (3)It is important to note that the Sommerfeld integral (3) is conditionally convergent and asstated, requires that y (cid:54) = y .In the Sommerfeld approach ([8]), we assume the upward scattered field u s in the toplayer can be expressed as u s = 14 π (cid:90) ∞−∞ e − √ λ − k y (cid:112) λ − k e iλ ( x − x ) σ ( λ ) dλ, (4)where σ ( λ ) is an unknown density on the upper interface y = 0. It is straghtforward toverify that u s satisfies the Helmholtz equation with Helmholtz parameter k .In the second layer, the scattered field u s can be written in terms of contributions fromboth the upper ( y = 0) and lower ( y = − d ) interfaces: u t and u b . These are given by u t = 14 π (cid:90) ∞−∞ e √ λ − k y (cid:112) λ − k e iλ ( x − x ) σ +2 ( λ ) dλ, (5) u b = 14 π (cid:90) ∞−∞ e − √ λ − k ( y + d ) (cid:112) λ − k e iλ ( x − x ) σ − ( λ ) dλ, (6)where σ +2 ( λ ) and σ − ( λ ) are used to denote spectral density functions on the upper andlower interfaces.Similarly, we can represent the scattered field u s in the third layer with an unknowndensity σ ( λ ) on the lower interface as u s = 14 π (cid:90) ∞−∞ e √ λ − k ( y + d ) (cid:112) λ − k e iλ ( x − x ) σ ( λ ) dλ. (7) Remark . The signs of the terms e ± √ λ − k i y and e ± √ λ − k i ( y + d ) in eqs. (4)-(7) ensurethat evanescent modes (when | λ | > | k i | ) decay away from each layer. (Physically, this isrelated to causality and is required in the derivation of formula [25] by countour integration).It is worth noting that the four unknown densities σ , σ +2 , σ − and σ can be interpretedin two ways. First, they can simply be considered the spectral densities in the Fourierdomain of a consistent representation for the Helmholtz equation. For those more familiarwith potential theory, they can be viewed as the Fourier transforms of charge densities offour single layer potentials lying on the corresponding interfaces [3].In the absence of any inclusions, the Sommerfeld representation for the field in eachsubdomain is derived from a “mode by mode” analysis. That is, the unknown functions σ , σ +2 , σ − , and σ are found by enforcing the continuity conditions at the interface for each alue 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 , (8) (cid:20) ∂u∂n (cid:21) = 0 , (9)where [ · ] denotes the jump of a function along the interface, ∂/∂n is the normal derivativeand u is the total field in each layer [10].It is straightforward to check that the linear system to be solved for each λ takes theform: √ λ − k − √ λ − k − e − √ λ − k d √ λ − k e − √ λ − k d √ λ − k √ λ − k − √ λ − k − e − √ λ − k d e − √ λ − k d − − σ ( λ ) σ +2 ( λ ) σ − ( λ ) σ ( λ ) = − e − √ λ − k y √ λ − k e − √ λ − k y (10) Definition 2.1.
We will denote the 4 × A λ .For the problem we consider here, the Sommerfeld integrals must be coupled to arepresentation of the field induced by the many particles present in the central layer. Beforediscussing the coupled system, however, we first summarize some well-known facts aboutscattering from a finite collection of inclusions in a homogeneous infinite medium. Suppose now that we have an inclusion of dielectric material with k = ω √ (cid:15) p µ embeddedin R , assumed to consist of a dielectric with k = ω √ (cid:15) µ . For transverse magnetic(TM)polarization, the total electrical field u in the exterior of the inclusion satisfies the Helmholtzequation: ∆ u + k u = 0 . (11)Further, the total field u can be written as the sum of the incident field u inc and thescattered field u s , where u s satisfies (11) and the Sommerfeld radiation condition,lim r →∞ √ r (cid:18) ∂u s ∂n − ik u s (cid:19) = 0 , (12)where r = (cid:112) x + y . Within the inclusion, the field u satisfies the Helmholtz equationwith wavenumber k p , ∆ u + k p u = 0 . (13)On the boundary of the inclusion, we must enforce the continuity conditions given by Eq(8) and (9). .1 A single disk When the inclusion is a disk of radius R centered at the origin, it is straightforward torepresent the solution using separation of variables, with u s = ∞ (cid:88) n = −∞ β n H n ( k r ) e inθ (14)in the exterior and u = ∞ (cid:88) n = −∞ γ n J n ( k p r ) e inθ (15)in the interior. Here, ( r, θ ) are the polar coordinates of a point in the plane, H n ( r ) is theHankel function of the first kind of order n and J n ( r ) is the Bessel function of order n [11, 24].We now expand the incident wave u inc and its normal derivative in the form: u inc = ∞ (cid:88) n = −∞ α n J n ( k r ) e inθ , ∂u inc ∂r = ∞ (cid:88) n = −∞ α n k J (cid:48) n ( k r ) e inθ . (16)Enforcing the continuity conditions (8), (9) on the boundary of the disk for each Fouriermode, we easily obtain the following linear equation for mode n : (cid:20) − H n ( k R ) J n ( k p R ) − k H n (cid:48) ( k R ) k p J (cid:48) n ( k p R ) (cid:21) (cid:20) β n γ n (cid:21) = (cid:20) α n J n ( k R ) α n k J (cid:48) n ( k R ) (cid:21) , (17)where n ∈ N .Solving Eq. (17) determines the coefficients β n , γ n : β n = (cid:20) k p J n ( k R ) J (cid:48) n ( k p R ) − k J (cid:48) n ( k R ) J n ( k p R ) k H n (cid:48) ( k R ) J n ( kR ) − k p J (cid:48) n ( k p R ) H n ( k R ) (cid:21) α n , (18) γ n = (cid:20) k J n ( k R ) H n (cid:48) ( k R ) − k J (cid:48) n ( k R ) H n ( k R ) k H n (cid:48) ( k R ) J n ( k p R ) − k p J (cid:48) n ( k p R ) H n ( k R ) (cid:21) α n . (19)It is straightforward to verify that the denominator in the preceding expressions cannotvanish if k and k p have positive real part and non-negative imaginary part [10, 21]. Definition 3.1.
The mapping between the incoming coefficients { α n } and outgoing coef-ficients { β n } is referred as the scattering matrix for the disk and denoted by S . Remark . While we restrict our attention here to dielectric particles, the method caneasily be extended to perfectly conducting disks. Since the interior field u in (15) is zerofor perfect conductors, from (17), we have β n = − J n ( k R ) H n ( k R ) α n . (20) Remark . In the remainder of this paper, we will refer expansions based on Hankelfunctions, such as (14), as multipole expansions or H -expansions, and expansions basedon Bessel functions, such as (15), as local expansions or J -expansions. Remark . In practice, we will truncate the expansions after, say, p terms with the valueof p to be determined later. We then define (cid:126)α ≡ ( α − p , α − p +1 , . . . , α , α , . . . , α p ) and (cid:126)β ≡ ( β − p , β − p +1 , . . . , β , β , . . . , β p ). .2 Multiple disks Suppose now that we have M well separated, identical dielectric disks randomly distributedin a homogeneous medium. Each disk is assumed to have radius R and wavenumber k p andthe background medium again has wavenumer k . For each individual particle, the analysiscan be carried out as above. We will denote by (cid:126)α m the incoming coefficients and by (cid:126)β m the outgoing coefficients for the m -th particle. We have (cid:126)β m = S p [ (cid:126)α m ] , for m = 1 , · · · , M. (21)where S p denotes the truncated (2 p +1) × (2 p +1) scattering matrix acting on the truncatedexpansion.The principle difference between the single particle and multi-particle scattering problemis that, in the latter case, the incoming field experienced by each particle consists of twoparts: the (applied) incident field u inc and the contribution to the scattered field u s fromall of the other particles. In order to formulate the problem concisely, given the multipoleexpansion for disk j , we need some additional notation. Lemma 3.1. [28] Let disk m be centered at x m and let disk l be centered at x l . Then themultipole expansion ∞ (cid:88) n = −∞ β mn H n ( k r m ) e inθ m (22) induces a field on disk l of the form u = ∞ (cid:88) n (cid:48) = −∞ α ln (cid:48) J n (cid:48) ( k r l ) e in (cid:48) θ l (23) where α ln (cid:48) = ∞ (cid:88) n = −∞ e − in ( θ lm − π ) β mn (cid:48) − n H n ( k (cid:107) x m − x l (cid:107) ) . Here, ( r m , θ m ) and ( r l , θ l ) denote the polar coordinates of a target point with respect to diskcenters x m and x l , respectively and θ lm denotes the angle between ( x m − x l ) and the x -axis.Remark . We denote by T jm the translation operator that maps the outgoing coefficients (cid:126)β m from particle m to the local expansion (cid:126)α l centered at particle l . With this operator inplace, the incoming coefficients (cid:126)α m for the m -th particle is (cid:126)α m = (cid:126)a m + M (cid:88) j =1 j (cid:54) = m T jm (cid:126)β j , (24)where (cid:126)a m is the (truncated) local expansion (16) of the incident wave u inc on particle m . T jm is referred to as the multipole-to-local (M2L) translation operator [28].Combining eqs. (21) and (24), one can easily eliminate the incoming coefficients (cid:126)α m and obtain the following linear system that only involves the outgoing coefficients: (cid:0) S − − T (cid:1) (cid:126)β (cid:126)β ... (cid:126)β M = (cid:126)a (cid:126)a ... (cid:126)a M , (25)where S = S p S p . . . S p , T = T · · · T M T · · · T M ... ... . . . ... T M T M · · · . he system (25) can be solved iteratively, using GMRES [29]. Since each translationoperator T nm is dense, a naive matrix-vector product requires O (( M (2 p + 1)) ) operations,where p is the order of the truncated expansion. FMM acceleration reduces the cost to O ( M (2 p + 1) ) work, for which we refer the reader to [28, 7]. Further, (25) has a simplediagonal preconditioner. Multiplying through by the block diagonal matrix S , results in thepreconditioned system matrix I − ST . This significantly reduces the number of iterations.We now extend the multiple scattering approach to arbitrarily shaped particles. Remark . It is worth emphasizing that the multiple scattering theory as discussed here ishardly new. We refer the reader to [14, 18, 31] and the references therein.
When the dielectric inclusions are of arbitrary shape, multiple scattering theory cannot beused quite so easily. Suppose, however, that an inclusion Ω is compactly supported withboundary ∂ Ω and that it is composed of a homogeneous material with wavenumber k p , asabove. Given the incident wave u inc and the boundary conditions (8), (9), the exteriorscattered field u s and the field u within Ω have the following representation [10]: u s = S k σ + D k µ, for x ∈ Ω c , (26) u = S k p σ + D k p µ, for x ∈ Ω , (27)where S k and D k are the usual single layer and double layer potentials on ∂ Ω, S k σ = (cid:90) ∂ Ω G k ( x , y ) σ ( y ) ds y , (28) D k µ = (cid:90) ∂ Ω ∂G k ( x , y ) ∂n ( y ) µ ( y ) ds y . (29) σ ( y ) and µ ( y ) are unknown charge and dipole densities that lie on the boundary ∂ Ω. Wewill need the normal derivatives of S k and D k as well: N k σ = (cid:90) ∂ Ω ∂G k ( x , y ) ∂n ( x ) σ ( y ) ds y , T k µ = (cid:90) ∂ Ω ∂ G k ( x , y ) ∂n ( x ) ∂n ( y ) µ ( y ) ds y . (30)By construction, the representations (26) and (27) satisfy the relevant Helmholtz equa-tion in each domain. The single layer potential S k is weakly singular and the value iswell-defined for x ∈ ∂ Ω. The operators D k and N k are define on the boundary in theprincipal value sense (and have different limits when approaching the boundary from theinterior and the exterior). The operator T k is hypersingular with its value on the boundarydefined in the Hadamard finite part sense. For further details, we refer the reader to [10].Enforcing the interface conditions (8), (9) and taking appropriate limits [10] yields thefollowing system of Fredholm integral equations of the second kind: µ + [ S k − S k p ] σ + [ D k − D k p ] µ = − u inc , (31) − σ + [ N k − N k p ] σ + [ T k − T k p ] µ = − ∂u inc ∂n . (32) Remark . It is worth noting that, while T k is hypersingular, the difference kernel T k − T k p is only logarithmically singular and compact as are all the other difference operators in (32),at least for smooth boundaries. We use Nystr¨om discretization for the system of equationsbased on the high order hybrid Gauss-trapezoidal rule of Alpert [1]. In this paper, werestrict our attention to smooth inclusions that are about one wavelength in size, so that 12digits of accuracy are easily achieved with modest values of N using the Gauss-trapezoidal k k D m D l ⌦ l ⌦ m Wednesday, July 9, 14
Figure 2: Two inclusions and their enclosing disks. The scattering matrix S i for each inclusionΩ i with wavenumber k p is defined as the map from an incoming field on D i to the correspondingoutgoing field. It is computed by solving a sequence of boundary value problems on the inclusionitself in a precomputation phase (see text). In this paper, we assume that all the inclusions areidentical but may be rotated, as drawn here. rule for logarithmic singularities of order 16. We refer the reader to [3] and the referencestherein for further details.The integral equation (32) was introduced in electromagnetics by M¨uller [23], and inthe scalar case by Kress, Rokhlin, Haider, Shipman and Venakides [19, 21, 27]. Suppose now that we have M inclusions Ω , . . . , Ω M that are identical up to rotation, andwell separated in the sense that each inclusion Ω i lies within a disk D i of radius R so thatthe disks are not overlapping. (see Fig. 2).In that case, we can sample the incoming field on the disk D j rather than Ω j as u = p (cid:88) n = − p α n J n ( k r ) e inθ , (33)using a polar coordinate system centered on the disk D j .Let σ n and µ n denote the solution to the integral equation (32) with right-hand side u inc = J n ( kr ) e inθ , ∂u inc ∂n = kJ (cid:48) n ( kr ) e inθ . We may then precompute the multipole expansionfrom these source distributions u = p (cid:88) l = − p β nl H l ( k r ) e ilθ , (34)where β nl = (cid:90) ∂ Ω j [ J l ( k | y | ) e − ilθ j ( y ) σ n ( y )] + n · ∇ [ J l ( w | y | ) e − ilθ j ( y ) µ n ( y )] ds y . (35)Here, y is the location of a point on ∂ Ω j with respect to the center of disk D j and θ j ( y )is the polar angle subtended with respect to the center of disk D j . The formula for β l isstandard [28, 7] and derived from the Graf addition theorem [24]. Definition 4.1.
As before, the mapping between the incoming coefficients { α n } andoutgoing coefficients { β n } is referred as the scattering matrix for the inclusion Ω j anddenoted by S j . he 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 shapesare the same up to rotation. This permits us to solve only 2 p + 1 integral equations ona single prototype inclusion in the enclosing disk. The scattering matrix for each rotatedcopy is then trivial to construct. Moreover, we can easily store the densities σ n and µ n ,since this requires only O (2 N (2 p + 1)) storage, where N is the number of points used todiscretize the boundary ∂ Ω. The amount of memory required to store the scattering matrixis O ((2 p + 1) ). 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 ) work. Each right-hand side corresponding to u inc = J n ( k R ) e inθ and ∂u inc ∂n = k J (cid:48) n ( k R ) e inθ can then be solved for n = − p, . . . , p at a total cost of O ( N (2 p +1))work. If we were interested in solving the multiple scattering problem in an infinite medium, wecould now proceed as in the previous section. The number of degrees of freedom is only2 p + 1 per inclusion rather than N points per inclusions (the number needed to discretizethe domain boundaries ∂ Ω j ). For complicated inclusions, this permits a vast reduction inthe number of degrees of freedom required and forms the basis for the FMPS method [15].Moreover, the block-diagonal preconditioned multiple scattering equations are much betterconditioned than the integral equation (31), (32) itself and FMM acceleration is particularlyfast in this setting. Remark . Extending the method to more than one type of substructure is straightforwardas long as the assumption that the enclosed circles are well separated still holds. Theadditional cost is the bookkeeping for different scattering matrices of these substructures.
To this point, we have discussed the layered medium and multiple scattering problemspearately. For the full problem, we now assume that multiple inclusions have been placedin 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 ( x ) = G k ( x , x ) + u s u ( x ) = u t + u b + M (cid:88) j =1 p (cid:88) n = − p β mn H n ( k r m ) e inθ m (36) u ( x ) = u s where u and u denote the fields in the top and bottom half spaces and u denotes thefield in the central layer exterior to the scattering disks D j . u s , u t , u b , and u s are theSommerfeld integrals from Section 2. Once u is known, the field within the scatteringdisks and the inclusions themselves is easily obtained.It remains to discuss the discretization of the Sommerfeld integral, and the setup of theglobal linear system for the unknowns σ , σ +2 , σ − , σ , and { (cid:126)β m , m = 1 . . . , M } . Let us consider the function u t defined by (6). Its computation is a standard problem inacoustic and electromagnetic scattering and often handled by contour deformation. It istypical to deform the integration contour by pushing it from the real line into the second λ plane: Each segment in the contour isdiscretized using Gauss-Legendre quadrature. The branch cut (shown in red) points upwardfrom k and downward from − k . and fourth quadrants of the complex λ -plane in order to avoid the square root singularitiesin the integrand. One option is to use a hyperbolic tangent contour [3], which yields spectralaccuracy 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 slightlyless efficient, but will permit us to evaluate the Sommerfeld integral using the non-uniformFFT, as explained further below. The contour consists of three segments: Γ , Γ and Γ ,where Γ : t − ib, t ∈ (0 , ∞ ) , Γ : it, t ∈ [ − b, b ] , Γ : t + ib, t ∈ ( −∞ , . (37)The branch cuts for the square root in the integrand are chosen to ensure that waves aredecaying away from the interface. (up at k and down at − k as shown in Fig. 3).We truncate Γ and Γ at a point t max >
0, where the integrand of u t has decayed toa user-specified tolerance. Fortunately, the decay in the integrand is exponential once λ exceeds k . (The precise rate of decay depends on the distance from the interface of thescattering disks and the point source generating the incoming field.) We let N S denotethe number of points used in the quadrature for the Sommerfeld contour and note thateach discretization point λ j on the contour corresponds to a plane wave. We use the samecontour and the same N S values { λ j } for each of u s , u t , u b , and u s . Let us denote by (cid:126)σ the discretized densities on the dielectric layers, (cid:126)σ = [ (cid:126)σ , (cid:126)σ +2 , (cid:126)σ − , (cid:126)σ ] T ,and by (cid:126)β the multipole coefficients for all M particles in the central layer. Each of (cid:126)σ , (cid:126)σ +2 , (cid:126)σ − , and (cid:126)σ is of length N S and the full linear system for multiple scattering in the layeredmedium takes the form of a block 2 × (cid:20) A BC D (cid:21) (cid:20) (cid:126)σ(cid:126)β (cid:21) = (cid:20) b (cid:21) . (38) A itself is block diagonal 4 N S × N S matrix with 4 × A λ in (10),each such block corresponding to a distinct λ j in the contour integral discretization. Theright-hand side component b is simply the right-hand side of (10) for each such λ j . Thematrix D = ST − 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 translatesthe multipole expansion coefficients to a Sommerfeld representation on the upper and lower nterfaces of the layered medium, while C requires the evaluation of the Sommerfield integralcontributions from the interfaces in terms of incoming local expansions on the scatteringdisks themselves. We turn now to the efficient application of the matrices B and C . A straightforward mechanism to map from the (cid:126)σ variables to local expansions on the M disks is to use the Jacobi-Anger formula [24]. Lemma 5.1.
Given r ∈ R , k ∈ C , we have e ikr cos θ = ∞ (cid:88) n = −∞ i n J n ( kr ) e inθ . (39)Suppose now that we want to compute the contribution from σ +2 to a local expansionon a disk centered at ( x , y ). Using Lemma 5.1, it is easy to see that e √ λ j − k y + iλ j ( x − x ) = e √ λ j − k y + iλ j ( x − x ) ∞ (cid:88) n = −∞ i n J n ( k r ) e in ( φ + θ ) , (40)where φ = arccos( λ j /k ), θ = arccos(( x − x ) /r ) and r = (cid:112) ( x − x ) + ( y − y ) . Theanalogous formula can be obtained for the contribution from σ − .The cost of using formula (40) to compute the action of the C block in the systemmatrix above is clearly O ( M N S (2 p + 1)), where M denotes the number of particles and N S the number of discretization points λ j in the Sommmerfeld contour and p is the order of theexpansions used in the multiple scattering representation. This is quite acceptable wheneither N S or M is small. For high frequency problems with many inclusions, where k islarge and N S = O ( k ), we have developed a more efficient scheme, based on the nonuniformFFT (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 t for a fixed value of y , we have14 π (cid:90) Γ e √ λ − k y (cid:112) λ − k e iλ ( x − x ) σ +2 ( λ ) dλ = 14 π e b ( x − x ) (cid:90) t max g ( t ) e itx dt, (41)where g ( t ) = e √ ( t − ib ) − k y (cid:112) ( t − ib ) − k e − itx σ +2 ( t − ib ) . 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 theSommerfeld 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, thenonuniform FFT (NUFFT) of Dutt and Rokhlin [12, 13] permits this to be done in nearlylinear time. In our numerical simulations, we use the version discussed in [17, 22]. Theanalogous method permits the rapid evaluation of the Sommerfeld integral on the contourΓ . For the integral on Γ , the NUFFT cannot be applied, but only a few discretizationpoints are required, so we evaluate that contribution directly.To provide rapid access to the field induced by the Sommerfeld integral at any locationin the central layer, we superimpose on it a grid of n × n boxes that contain all of the M cattering disks. In each such box, we construct a tensor product m × m Chebyshev mesh,which will permit q th order local interpolation by barycentric interpolation [4]. The cost forevaluation at all grid points is O (( n m ) ( n m + N S ) log( n m + N S )) operations, usingthe NUFFT for each of the distinct n m locations in y .Consider now one of the scattering disks D j of radius R . If we discretize the boundaryof the disk using 2 p + 1 equispaced points, evaluation of the induced field at each of thepoints requires O ( m m ) operations, for a net cost of O ( m m (2 p + 1)) work. An FFT oforder (2 p + 1) converts these values into their Fourier transforms, which we denote by { u n } ,for n = − p, · · · , p . From this, the n -th term a n in the incoming J -expansion is simply a n = u n J n ( k R ) . (42) Remark . The formula (42) will fail if the value k R is a zero of the function J n for any n from − p . . . , p . This can be avoided if we also compute the normal derivative of theSommerfeld integral on the boundary of each scattering disk. If we denote by { u (cid:48) n } theFourier coefficient of the normal derivative, it is easy to see that a n = u n J n ( k R ) + u (cid:48) n kJ (cid:48) n ( k R ) J n ( k R ) + ( kJ (cid:48) n ( k R )) , for j = − p, · · · , p. (43)The evaluation of the gradient of the Sommerfeld integral can be computed by an obviousmodification of the formula (41) or (with a reduction in order) by computing the gradientof the tensor product Chebyshev series discussed above.In summary, it requires O ( M m m (2 p + 1)) operations to interpolate the field valueson each of the M scattering disks and O ( M (2 p + 1) log(2 p + 1)) operations to obtain thecoefficients of the J -expansions. This completes the computation of the C block in thesystem matrix. The off-diagonal B block in (38) requires a formula for recasting the multipole expansionto the corresponding Sommerfeld representation on either the upper or lower interface ofthe layered medium. More precisely, each H -expansion in the central layer, centered ondisk D j with center ( x j , y j ) has a spectral representation on the upper layer y = 0 and thelower layer y = − d of the form: u tj = 14 π (cid:90) ∞−∞ (cid:112) λ − k e iλ ( x − x ) σ + mp ( λ ) dλ, (44) u bj = 14 π (cid:90) ∞−∞ (cid:112) λ − k e iλ ( x − x ) σ − mp ( λ ) dλ, (45)respectively.The formulae for σ + mp ( λ ) and σ − mp ( λ ) follow directly from the following theorem. Theorem 5.2. [6] Let ( x j , y j ) denote the center of a multipole expansion in the centrallayer, with − d < y j < and let ( r, θ ) denote the polar coordinates of a target point withrespect to that center. Then, on the upper interface, H n ( kr ) e inθ = ( − n π (cid:90) ∞−∞ e √ λ − k y j √ λ − k e iλ ( x − x j ) (cid:18) √ λ − k + k k (cid:19) n dλ, (46) and on the lower interface, H n ( kr ) e inθ = ( − n π (cid:90) ∞−∞ e −√ λ − k ( d + y j ) √ λ − k e iλ ( x − x j ) (cid:18) √ λ − k − k k (cid:19) n dλ. (47) ach multipole coefficient in the expansion about disk D j contributes to each of the N S discretization points in the Sommerfeld integrals, requiring a total of O ((2 p + 1) N S M )work. This, then, is the cost of applying the B block of the system matrix directly. Because of the computational complexity of applying the B block in the manner describedabove, 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 thecenters of the H -expansions lie at the nodes of a uniform grid in the central layer and letus consider the contributions from the n th mode at each such grid point ( x l , y j ) for a fixedhorizontal line y = y j . If there are n such expansion centers, with x coordinates { x l } , l = 1 , · · · , n , and we denote by a ln the coefficient for the n th mode of the H -expansion atlocation ( x l , y j ), then the total contribution to the induced spectral coefficient σ + mp ( λ j ) onthe top layer is given by { σ + mp ( λ j ) } n := e √ λ j − k y j (cid:18) (cid:113) λ j − k + k k (cid:19) n n (cid:88) l =1 a ln e − iλ j ( x j − x ) (48) { σ − mp ( λ j ) } n := e − √ λ j − k ( d + y j ) (cid:18) (cid:113) λ j − k − k k (cid:19) n n (cid:88) l =1 a ln e − iλ j ( x j − x ) The formulae (48) imply that for each row, one can use the NUFFT to compute theinduced coefficients for each discrete quadrature node λ j on Γ or Γ . As above, we usedirect computation for the contributions to discretization nodes on Γ . In the general case,the centers of the H -expansions are not aligned on a grid, but we can first shift the centerof each H -expansion to the nearest grid point, using the multipole-to-multipole translationoperator [28, 7] based on the Graf addition theorem [24]. After M such shifts, we mayapply the transformation of (48).The total computational cost is O ( M (2 p + 1) ) for shifting all the H -expansions and O ( n (2 p + 1) ( n + N S ) log( n + N S )) for the NUFFT-based work (see Table 1). Themerits of the NUFFT-based schemes would become more apparent for larger N S . Table 1: Comparison of CPU time in seconds for the Sommerfeld-to-local and multipole-to-Sommerfeld operators, unsing both the direct and NUFFT-based schemes (see text). TheSommerfeld contour is discretized with 500 Gauss-Legendre points (240 points for Γ and Γ ,with 20 points for Γ ). (a) Computation of the Sommerfeld-to-local operator Number of scatterers 100 500 1,000 5,000Direct method 2.21e-2 9.90e-2 2.06e-1 9.86e-1NUFFT 1.31e-1 1.77e-1 2.28e-1 4.93e-1 (b) Computation of the multipole-to-Sommerfeld operator
Number of scatterers 100 500 1,000 5,000Direct method 3.49e-2 1.82e-1 3.60e-1 1.80NUFFT 6.05e-2 1.38e-1 1.61e-1 2.77e-113 .4 Iterative solution of the system matrix
We will solve equation (38) using the iterative method GMRES [29]. However, the un-knowns (cid:126)σ and (cid:126)β may be poorly scaled with respect to each other. However, A is blockdiagonal, as noted above, with simple 4 × A directly and useGMRES on the Schur complement of (38). In other words, we solve the system[ D − CA − B ][ (cid:126)β ] = − CA − b (49)instead. This is much better conditioned and involves only the (cid:126)β unknowns. The Schurcomplement formalism has a simple physical interpretation: it is, in essence, a reformulationof the scattering problem using the layered medium Green’s function. In this section, we illustrate the performance of our algorithm with three examples. Forsimplicity, we use a single class of inclusions, parametrized by (cid:26) x = ( a + a cos( a t )) cos( t ) ,y = ( a + a cos( a t )) sin( t ) , for 0 ≤ t < π. (50)As discussed in section 4, inclusions with more complicated boundaries do not introduce anyessential difficulty in our scheme except that the precomputation of the scattering matrixis a little more involved, particulalry if corners are present [5, 20].Given a fixed a , a and a , multiple copies of the inclusion are randomly distributed inthe central layer of the medium with random orientations. To ensure the inclusions are wellseparated but confined in a fixed region, we use a bin sorting algorithm to construct therandom distribution. We begin with inclusions located on a regular grid and then perturbtheir positions randomly, accepting the random move if the inclusion remains inside theregion and stays well separated from the others. Several such sweeps are carried out torandomize the positions further.We have not, as yet, specified the parameter N S used to discretize the Sommerfeldintegral in (37). While special techniques have been developed by many authors to handlesources near the interface (see [3, 9, 25]), we simply assume that the source definingthe incoming field is at least 0.2 wavelengths from the top interface. More precisely, inour examples, the source point in the top layer is placed at (1 ,
1) (which is roughly 0 . k = 1). We also assume that the nearest the inclusionsget to either one of the interfaces in the layered medium is at least 0 . t max = max {| k | , | k | , | k |} + 20, b = 0 .
2, and discretize Γ andΓ using 240 Gauss-Legendre points and Γ by 20 Gauss-Legendre points. This is sufficientto achieve about 10 digits of accuracy.All computations are carried out using a 2.3GHz Intel Core i5 laptop, with 4GB RAM. In our first example, we consider the scattering of inclusions defined by parameters a =0 . , a = 0 .
04, and a = 3 in eq. (50) with wavenumber k p = 2 .
0. To obtain thescattering matrix with p = 10, we solve the integral equation (31) and (32) by discretizingthe boundary of the particle using N = 300 equispaced points. We assume the wavenumbersof the layered medium are given by k = 1 . , k = 3 . , k = 1 .
0. The thickness ofthe central layer is determined by the parameter d = 32. We consider distributions of M = 100 , , , , ,
000 inclusions and solve the mulitple scattering problem usingGMRES with FMM acceleration. We terminate the iteration once the residual is less than10 − . Results are presented in Fig. 4 and 5. ,
000 dielectric inclusions randomly distributed ina three-layered medium. The wavenumber for each particle is k p = 2 . k = 1 . k = 3 . k = 1 .
0. The size of each particle is approximately0 . k . (a) (b) Figure 5: Convergence behavior of GMRES and the CPU time required for various numbersof inclusions embedded in either (a) free space or (b)a three layered medium. For (a), we set k = k = k = 3 . k = 1 . k = 3 . k = 1 . ig. 4 shows the total field in the case M = 5 , k = k = k ). As expected, convergence ismore rapid than when the inclusions are embedded in a true layered medium, because ofthe multiple reflections from the interfaces themselves. In our second example, we consider the same inclusion shape as above, with k p = 2, butwith higher contrast materials. We fix the number of particles to be 200 and the thicknessof the middle layer to be d = 12. We allow the wavenumber in the middle layer to varyfrom k = 1 up to k = 20. 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 layermedium (b). Note that the convergence is slower at high constrast and that this effectis more pronounced in the layerd medium case, where the central layer involves strongscattering and reflection. In our last example, we consider the scattering from a different inclusion shape, setting a = 0 . , a = 0 . , a = 5 in eq. (50) with k p = 2 .
0. The inclusions are smoothedpentagons, as shown in Fig. 8. We discretize the boundary of the inclusion using N = 300equispaced 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, weset k = 1 . , k = 3 . , k = 2 .
0. Results are shown in Figs. 8 and 9.Note that in order to obtain 6 digits of accuracy, 69 . , M . We have developed a fast algorithm to simulate electromagnetic scattering from a mi-crostructured, three-layered material. Our methodology permits inclusions of arbitraryshape using a scattering matrix formalism combined with the use of Sommerfeld integralsto account for the influence of the layered material. We have designed efficient proceduresto evaluate the Sommerfeld integral at arbitrary locations in the layered material usingthe non-uniform FFT and an effective preconditioner that allows the multiple scatteringproblem to be solved using GMRES with a modest number of iterations. As one wouldexpect from physical considerations, the performance of the method degrades when thepacking of inclusions is dense and when the contrast is high. While the method is suitablefor parallel implementation, we are also investigating the possibility of replacing GMRESiteration with a fast direct solver [16].Extension of the present method to the quasi-periodic case, where the incoming fieldimpinges on a periodic microstructure will be reported at a later date.
Acknowledgements
This work was supported in part by the Applied Mathematical Sciences Program of theU.S. Department of Energy under Contract DEFGO288ER25053 and by the Office of the k = 1 . k = 10 . k = 1 .
0. For each inclusion, the wavenumberis k p = 2 .
0. The inclusions are approximately 0 . k . (a) (b) Figure 7: Convergence behavior of GMRES iteration and the CPU time required for 200inclusions embedded in the central layer, where k is allowed to vary from 1 to 20. In (a), wecreate a homogeneous background by setting k = k = k , while in (b), k and k are fixed at1, and k varies. 17igure 8: Real part of the total field when 1 ,
000 inclusions are embedded in a three-layeredmedium with k = 1 . k = 3 . k = 2 .
0. Each inclusion is a smoothed pentagon,approximately 0 . (a) (b) Figure 9: Convergence behavior of GMRES for a tolerance of 10 − and the CPU timerequired as the number of inclusions embedded in the central layer varies. In (a), we createa homogeneous background by setting k = k = k = 3 .
0, while in (b), k = 1 . k = 3 . k = 2 . ssistant Secretary of Defense for Research and Engineering and AFOSR under NSSEFFProgram Award FA9550-10-1-0180. References [1] B. K. Alpert. Hybrid Gauss-trapezoidal quadrature rules.
SIAM J. Sci. Comput. ,20:1551–1584, 1999.[2] G. Bao and J. Lai. Radar cross section reduction of a cavity in the ground plane.
Commun. Comput. Phys. , 15:895–910, 2014.[3] A. Barnett and L. Greengard. A new integral representation for quasi-periodicscattering problems in two dimensions.
BIT Numerical Mathematics , 51(1):67–90,2011.[4] J.-P. Berrut and L. N. Trefethen. Barycentric lagrange interpolation.
SIAM Rev ,46:501–517, 2004.[5] J. Bremer, V. Rokhlin, and I. Sammis. Universal quadratures for boundary integralequations on two-dimensional domains with corners.
J. Comput. Phys. , 229:8259–8280,2010.[6] H. Cheng, J. Huang, and T. J. Leiterman. An adaptive fast solver for the modifiedhelmholtz equation in two dimensions.
Journal of Computational Physics , 211(2):616– 637, 2006.[7] H. Cheng, C. W. Y., Z. Gimbutas, L. Greengard, J. Huang, V. Rokhlin, N. Yarvin,and J. Zhao. Remarks on the implementation of the wideband fmm for the helmholtzequation in two dimensions.
Contemporary Mathematics , 408:99–110, 2006.[8] W. C. Chew.
Waves and fields in inhomogeneous media . IEEE Press, New York, 1995.[9] M. H. Cho and W. Cai. A parallel fast algorithm for computing the Helmholtz integraloperator in 3-d layered media.
Journal of Computational Physics , 231(17):5910–5925,2012.[10] D. Colton and R. Kress.
Integral Equation Method in Scattering Theory . Wiley-Interscience, New York, 1983.[11] D. Colton and R. Kress.
Inverse Acoustic and Electromagnetic Scattering Theory,Applied Mathematical Sciences 93 . Springer-Verlag, Berlin, 1998.[12] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data.
SIAM Journalon Scientific Computing , 14(6):1368–1393, 1993.[13] A. Dutt and V. Rokhlin. Fast fourier transforms for nonequispaced data, II.
Appliedand Computational Harmonic Analysis , 2(1):85 – 100, 1995.[14] L. L. Foldy. The multiple scattering of waves. i. general theory of isotropic scatteringby randomly distributed scatterers.
Phys. Rev. , 67:107–119, Feb 1945.[15] Z. Gimbutas and L. Greengard. Fast multi-particle scattering: A hybrid solver for theMaxwell equations in microstructured materials.
J. Comput. Phys. , 232:22–32, 2013.[16] K. L. Greengard, L. Ho and J.-Y. Lee. A fast direct solver for scattering from periodicstructures with multiple material interfaces in two dimensions.
J. Comput. Phys. ,258:738–751, 2014.[17] L. Greengard and J. Lee. Accelerating the nonuniform fast fourier transform.
SIAMReview , 46(3):443–454, 2004.[18] N. A. Gumerov and R. Duraiswami. A scalar potential formulation and translationtheory for the time-harmonic maxwell equations.
Journal of Computational Physics ,225(1):206 – 236, 2007.
19] M. Haider, S. Shipman, and S. Venakides. Boundary-integral calculations of two-dimensional electromagnetic scattering in infinite photonic crystal slabs: Channeldefects and resonances.
SIAM J. Appl. Math. , 62:2129–2148, 2002.[20] J. Helsing and R. Ojala. Corner singularities for elliptic problems: Integral equations,graded meshes, quadrature, and compressed inverse preconditioning.
Journal ofComputational Physics , 227(20):8820 – 8840, 2008.[21] R. Kress and G. F. Roach. Transmission problems for the helmholtz equation.
Journalof Mathematical Physics , 19(6):1433–1437, 1978.[22] J. Lee and L. Greengard. The type 3 nonuniform FFT and its applications.
Journalof Computational Physics , 206(1):1 – 5, 2005.[23] C. M¨uller.
Foundations of the Mathematical Theory of Electromagnetic Waves .SPringer Verlag, 1969.[24] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark.
NIST Handbook ofMathematical Functions . Cambridge University Press, New York, 2010.[25] 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):1 –13, 2014.[26] W. J. Parnell, I. D. Abrahams, and P. R. Brazier-Smith. Effective properties ofa composite half-space: Exploring the relationship between homogenization andmultiple-scattering theories.
The Quarterly Journal of Mechanics and AppliedMathematics , 63(2):145–175, 2010.[27] V. Rokhlin. Solution of acoustic scattering problems by means of second kind integralequations.
Wave Motion , 5:257–272, 1983.[28] V. Rokhlin. Rapid solution of integral equations of scattering theory in two dimensions.
J. Comput. Phys. , 86:414–439, 1990.[29] Y. Saad and M. Schultz. Gmres: A generalized minimal residual algorithm for solvingnonsymmetric linear-systems.
IAM J. Sci. Stat. Comput. , 7:856–869, 1986.[30] Y. Wu and Z.-Q. Zhang. Dispersion relations and their symmetry properties ofelectromagnetic and elastic metamaterials in two dimensions.
Phys. Rev. B , 79:195111,May 2009.[31] Y. Xu. Electromagnetic scattering by an aggregate of spheres.
Appl. Opt. , 34(21):4573–4588, Jul 1995., 34(21):4573–4588, Jul 1995.