A fast and robust solver for the scattering from a layered periodic structure containing multiparticle inclusions
Abstract
We present a solver for plane wave scattering from a periodic dielectric grating with a large number of inclusions lying in each period of its middle layer. Such composite material geometries have a growing role in modern photonic devices and solar cells. The highorder scheme is based on boundary integral equations, and achieves many digits of accuracy with ease. The usual way to periodize the integral equation—via the quasiperiodic Green’s function—fails at Wood’s anomalies. We instead use the freespace Green’s kernel for the near field, add auxiliary basis functions for the far field, and enforce periodicity in an expanded linear system; this is robust for all parameters. Inverting the periodic and layer unknowns, we are left with a square linear system involving only the inclusion scattering coefficients. Preconditioning by the singleinclusion scattering matrix, this is solved iteratively in time using a fast matrixvector product. Numerical experiments show that a diffraction grating containing inclusions per period can be solved to 9digit accuracy in under 5 minutes on a laptop.
Keywords: fast solver, periodic scattering, multiparticle scattering, layered medium, diffraction
1 Introduction
The modeling and design of periodic dielectric structures plays a central role in modern optics. Tools such as diffraction gratings, photonic crystals, metamaterials, plasmonics, and other microscale structures, are becoming key to efficient devices, including lasers, sensors, antireflective surfaces and absorbers [21], and solar cells [3]. For instance, in thinfilm solar cell design [44, 29] the use of periodic structures, and nanoparticle inclusions, in ordered or disordered composites, enhances absorption. One then seeks a grating structure with a specific arrangement of inclusions that maximizes absorption. Other optimization problems include the design of photonic crystal lenses [35]. Related is the inverse problem of inferring a structure from measurements [37, 5]. Such tasks demand a large number of solutions of the direct (forward) scattering problem. Similar periodic and multiparticle wave scattering problems arise in acoustics and elastodynamics, and in general whenever a supercell is used to approximate the response of a random composite material (e.g. [36]). Such considerations have spurred the development of efficient methods for solving Helmholtz and Maxwell frequencydomain boundary value problems in periodic geometries [21, 7, 14, 9, 10, 23, 13, 18]. High accuracy can be challenging to achieve due to guided modes, resonances, and extreme parameter sensitivity.
Therefore, in this paper we consider the monochromatic scattering from a layered periodic structure containing a large number of inclusions (“particles”) at given locations, as in a (generalized) photonic crystal. As shown in Fig. 1, the structure is periodic in the direction, layered in the direction and invariant along the direction. Because of the twodimensional (2D) geometry, there exist two fundamental polarizations in the electromagnetic scattering: transverse magnetic (TM) where the magnetic field is transverse to the plane, and transverse electric (TE) where the electric field is transverse to the plane. We will focus on TM polarization, noting that our technique applies to TE polarization without any essential difficulty.
The grating scattering problem has been mathematically very well studied. It has been proved that for an arbitrary periodic dielectric and incident angle the problem has a unique solution for all frequencies with the possible exception of a countable set of resonances (singular frequencies [11]) at which the solution is not unique. Such physical resonances are not to be confused with Wood’s anomalies (for the definition see the next section), which are frequencies where at least one of the Bragg diffraction orders points along the grating, i.e. in the direction. A Wood’s anomaly does not prevent the solution from being unique, although it does cause arbitrarily large sensitivity with respect to the incident wave angle or frequency [33], and also causes problems with certain integral equation methods [10]. One of the advantages of our scheme is that it is applicable and accurate at or near Wood’s anomalies, without any modifications.
There exists a wide range of numerical methods for periodic diffraction, including boundary integral equations [2, 14, 10, 23, 13, 18], finite element methods [4, 8], Fourier expansion based methods [38], and continuation methods [15]. In the time domain, the finite difference scheme has been discussed in [26]. The advantages of the integral approach over finite elements and finite differences are that it reduces the dimension by one (vastly reducing the number of unknowns), and achieves highorder accuracy with appropriate surface quadratures. However, the resulting linear system is often dense, making a naive matrixvector product expensive when the number of unknowns is large. In this paper, we will reduce this cost via the fast multipole method (FMM) [24].
More specifically, we propose an integral approach based on the free space Green’s function; this bypasses the considerable complexities of computing the periodic Green’s function [30, 13]. We split the representation of the scattered field in the grating structure into near field and far field components. The near field is represented by standard freespace Helmholtz single and doublelayer potentials on the material interfaces, while the far field is taken care by a local expansion (Fourier–Bessel or expansion) whose coefficients are fixed by enforcing the periodic boundary condition explicitly in the linear system. This builds upon recent ideas of the last author and coworkers [9, 10, 18].
Solving for discretized layer densities on each of the inclusion boundaries would introduce an unnecessarily large number of unknowns. Hence, following [22, 32], we precompute the inclusion scattering matrices, then treat the set of outgoing scattering coefficients as a reduced set of unknowns. When particles are subwavelength, and not extremely close to each other, this is highly accurate with only 20 or so unknowns per particle [32]. The full rectangular linear system then couples these to the grating interface densities and periodizing expansion coefficients. By eliminating the last two (via a Schur complement and pseudoinverse) we are left with a square linear system for the particle scattering coefficients, which we precondition with a blockdiagonal matrix and then solve via GMRES with FMM acceleration, with effort scaling linearly in . The result is a robust, efficient, highorder accurate solver that we expect to be useful for design and optimization problems for periodic photonic devices.
The outline of the paper is as follows. Section 2 gives the mathematical formulation of the periodic problem. Section 3 proposes the integral approach for the scattering from a periodic structure without particle inclusions, based on the free space Green’s function. Section 4 reviews classical multiparticle scattering and discusses the evaluation of the scattering matrix. The quasiperiodizing scheme combining all the above techniques is given in Section 5, and numerical experiments are shown in Section 6. We draw conclusions in Section 7.
2 Problem formulation
Consider the planewave incident time harmonic scattering (with time dependence ) from a 2D periodic (or grating) structure with period . As shown in Figure 1, the unit cell consists of three layers, denoted by , and . Let and denote the two smooth interfaces separating the layers. The left and right boundaries of are denoted by and , . Assume the permittivity is given as , and in the three layers respectively. A large number of particles, collectively denoted by , with the same permittivity , are located inside . The permeability is assumed to be constant everywhere.
For TM polarization, in which case the total electric field is , the full time harmonic Maxwell equations
are reduced to the Helmholtz equation:
(1) 
where , and where the wavenumber takes one of four values,
(2) 
In the usual setting of scattering theory, the full wave is , where is the incident wave and is the resulting wave scattered from the periodic structure. The incident wave is a plane wave in , and elsewhere, with the incident angle. Since the wave is propagating in different layers, the continuity condition along various interfaces in TM polarization with constant permeability is:
(3) 
where denotes the jump of a function across the interface, is the normal derivative, and is the total field in each layer [20].
We use the term quasiperiodic if a function (such as ) satisfies
(4) 
where is the incident horizontal wavevector. The factor
is the Bloch phase associated with translation by one period. Since is quasiperiodic, we seek a scattered wave with this same symmetry, hence (4) also holds for the full wave [41]. Restricting to a single unit cell, with left and right walls and respectively, we have matching conditions
(5)  
(6) 
Finally, must satisfy a radiation condition. Let be sufficiently large such that lies between the lines and (see Figure 1). Define , , and , where and the sign of the squareroot is taken as positive real or positive imaginary. Then the plane wave with wavevector is quasiperiodic for each integer , and satisfies the Helmholtz equation at frequency . The radiation condition on is expressed by uniformly convergent outgoing or decaying Rayleigh–Bloch expansions:
(7)  
(8) 
The complex coefficients , , for such that (propagating waves), are the Bragg diffraction amplitudes at the grating orders. For all other these give evanescent components which do not contribute to the far field. A Wood’s anomaly occurs if, for some , , thus (upper halfspace), or if thus (lower halfspace). The radiation conditions ensure that is outgoing except at a Wood’s anomaly, when the th Rayleigh–Bloch mode in (7) or (8) is constant in (a horizontally traveling plane wave). It is also possible to have a double Wood’s anomaly, when for a pair of integers , .
Theorem 2.1.
Henceforth we will assume parameter values for which the solution is unique.
To summarize, we are interested in the solution of (1) together with the continuity condition (3), the unitcell quasiperiodicity (5)–(6), and the radiation condition (7)–(8) satisfied by . In the next section we discuss the solution for the layered periodic structure without inclusions, via the integral equation approach.
Remark 1.
From now on we slightly abuse the notation , and introduced, by replacing them by their truncated versions. In other words, has boundaries , , , and , while has boundaries , , , and . The unit cell will refer to . The artificial upper and lower boundaries and are also called transparent boundaries in some literature [7].
3 Robust solution for empty periodic layered structure
The standard approach to convert a nonperiodic integral equation formulation into the correct periodic one is to use the quasiperiodic Green’s function [34], defined at wavenumber as
(9) 
where is the target and the source point, and is the free space Green’s function at wavenumber , i.e. , where is the outgoing Hankel function of order zero, and is the Euclidean norm. (9) is welldefined away from Wood’s anomalies, and has been successfully applied in many grating problems [2, 14, 23, 41]. However, it has two major practical drawbacks: i) it is expensive to evaluate (requiring either series acceleration, or lattice sums [34]), and ii) it blows up (with an inversesquareroot singularity) at Wood’s anomalies, causing a purely numerical breakdown in the solution of a what remains a wellposed problem. The key idea is that (9) can be rewritten as
(10) 
where is a positive integer, is a fixed origin, is the angle of a vector , is the Bessel function of order , and the lattice sum coefficients can be found by Graf’s addition theorem [40, 10.23(ii)]. The second term accounts for the smooth field due to the infinite set of far sources and from (9). The direct terms account for the near field.
Remark 2.
The sum in (10) converges (and exponentially fast) if and only if the target is closer to the origin than the nearest “far” source, i.e. . This geometric condition is satisfied with for all if the region is not much taller than , and is placed near the center of . Thus, we will use in our numerical experiments. However, if is much taller than (high aspect ratio), then needs to be increased to guarantee uniform convergence.
Using this idea in the scattering setting we represent the scattered field in the first layer as,
(11) 
where is the th periodic translation of , i.e. , and are unknown periodic density functions defined on the interface , and is choice of origin. The coefficients are now unknown and need to be determined by the boundary conditions. and are the usual single and doublelayer potentials [19], which we may define living on a general interface at wavenumber by,
(12)  
(13) 
These representations satisfy the Helmholtz equation at wavenumber in ; thus the representation (11) satisfies the relevant Helmholtz equation in . When restricted to target points on these give the boundary integral operators (which is weakly singular), and (which is continous for smooth, and is to be interpreted in the principal value sense). Here and in what follows, the notation means the operator from source curve to target curve . We will also need the operators corresponding to the target normal derivatives on ,
(14) 
The operator is hypersingular and defined in the Hadamard finite part sense. The books [19, 20] give further details.
We also need the jump relations that relate the limiting values of (12)–(13) to the actions of the above boundary integral operators. Let and be the limiting values and normal derivatives approaching from the positive () or negative () side. Then for all continuous densities and ,
(15)  
(16)  
(17)  
(18) 
Thus the singlelayer potential is continuous for all , whereas the doublelayer is generally discontinuous across .
Turning to the third layer , we similarly represent the scattered field using layer potentials on ,
(19) 
where . The scattered field in the second layer has contribution from both and , thus
(20) 
where .
To determine the unknown densities , , , and the coefficients , and , we enforce the following boundary conditions according to (3) and (5)–(8).
Substituting the representations (11), (19) and (3) into conditions (21)–(23), one reaches a system of coupled integral and functional equations that can only be solved numerically, which means the interfaces must be discretized, and the infinite series truncated. The expansions are truncated up to order , i.e. they retain terms. The centers of the expansion , and should be located roughly in the centers of the domains , and . We discretize the four interfaces , , and using equally spaced points ( and are discretized through equallyspaced nodes in their parametrizations). The left and right boundary and are discretized by Gauss–Legendre nodes. The singular integrals involved in the layer potentials are discretized via the Nyström method [27], with 16thorder Alpert quadrature corrections [1]. A phase correction is applied to account for the Bloch phase factors when the parameter “wraps” around the end of the open curves or . The standard application of quadrature rules, and the use of the Alpert scheme, is described in [18, Sec. 2.5] (see the smooth case only).
In the end, stacking as block rows the equations (21), (22) and (23), the linear system for the discretized layered periodic structure without particle inclusions (the “empty” structure) takes the form,
(24) 
Here the righthand side vector has the form
(25) 
The unknown coefficient vector, which we will call , stacks the discretizations of the paired densities for interfaces , the coefficient vectors for the expansions in layers , and the Rayleigh–Bloch coefficient vectors and from (7)–(8). We may then summarize (24) by
We now describe the matrix blocks in (for readability we do not give formulae for every single block, trusting that their construction is unambiguous from the above; for more detail in a related scheme see [18]). Each block maps unknowns to values and normal derivatives at target nodes, or their phased differences between left and right walls. Matrix entries involve either freespace Green’s functions between source and target nodes, or or Rayleigh–Bloch expansions at target nodes.

and : Nyström selfinteraction matrices for and respectively, including the phased summation over source near neighbors. For instance,
(26) (where here and below, the Nyström discretizations of the various operators are implied), maps at the quadrature nodes on to the field jumps at these same nodes. indicates the by identity matrix. All terms (after subtractions) are compact apart from the identities, meaning that this subsystem is of Fredholm second kind. Recall that the identities originate in the jump relations (15)–(18). This BIE scheme for dielectric interfaces is due in the electromagnetic case to Müller [39], and the acoustic case to Kress–Roach [28] and Rokhlin [42].

, : interaction of ( summed) source densities on with value and normal derivatives on , and vice versa, at wavenumber .

, : values and normal derivatives of the middlelayer expansion on , .

, , and : phased differences of values and normal derivatives between the left and right walls and in (22) (case ), due to the layer potentials on , , and the middlelayer expansion respectively. , : phased differences for upper walls (case ) due to layer potentials on and the upper expansion. , : phased differences for lower walls (case ) due to layer potentials on and the lower expansion.

, (, ): values and normal derivatives of the layer potentials on (), and the upper (lower) expansions, evaluated on the upper boundary (lower boundary ).

, : values and normal derivatives of the Rayleigh–Bloch expansions on the upper and lower boundaries , respectively.
Remark 3.
The matrix is generally rectangular, depending on the specific numbers of discretization nodes. Although the upperleft block of is a square system coming from a Fredholm second kind system of BIEs, the other blocks involve expansions evaluated on interfaces and walls, which make as a whole exponentially illconditioned. Figure 2(a) shows the singular values of : there are many singular values clustered around , although the situation is alleviated a little by rescaling the columns that correspond to the expansions (also shown). However, such illconditioning is not an obstacle as long as the system is consistent: since is not too large (of typical size for corresponding to up to several wavelengths across one period ), we may use direct dense linear algebra for a smallnorm leastsquares solution. We use the mldivide command in MATLAB. Figure 2(b) shows the resulting real part of the scattered field with a total flux error (whose definition will be given in section 6) of .
4 Multiparticle scattering
Multiparticle scattering in free space has been discussed extensively in [22, 32, 31], with applications including climatology, remote sensing, and design of composite materials. Recently, we have developed a fast solver for finding the field scattered from a large number of particles located in a layered medium, by combining the Sommerfeld integral and multiple scattering theory. Here we briefly review the method introduced in [32], and then combine it with the periodic grating. See [32, 22] for more details. Note that our approach can be seen as a simple version of a reduced basis method [16].
4.1 Scattering matrix of a single particle
Consider for now a dielectric inclusion with wavenumber surrounded by uniform dielectric with . When the inclusion is a disk of radius centered at the origin, it is well known that the solution can be represented using separation of variables, with
(27) 
in the interior, and scattered wave
(28) 
in the exterior. Here, are the polar coordinates of a point in the plane, and is the Hankel function of the first kind of order .
Definition 4.1.
The mapping between the incoming coefficients and outgoing coefficients is referred as the scattering matrix. It will be denoted by , with matrix elements .
For a single disk, the scattering matrix is diagonal and is easily found analytically [19]; this is not true for an arbitrary inclusion shape. We instead seek a solution via BIE, using the Müller–Kress–Roach–Rokhlin scheme from the previous section. Suppose that the inclusion has boundary and is enclosed by a disk centered at the origin. Given the incident wave and the boundary conditions (3), the exterior scattered field and the field within have the following representations [19]:
(29)  
(30) 
Here and are unknown single and doublelayer densities on , and, in this section without ambiguity we drop the subscripts . Enforcing the interface conditions (3) and taking appropriate limits using jump relations (15)–(18) yields a system of Fredholm integral equations of the second kind:
(31)  
(32) 
4.2 Multiple inclusions
Suppose now that we have inclusions that are identical up to translation and rotation, and are well separated in the sense that each inclusion lies within a disk of radius such that the disks are not overlapping (see Fig. 3). The incident wave for the th inclusion may now be expanded as
(35) 
where are defined to be the polar coordinates of relative to the center of . We will denote by the set of incoming coefficients, and by the outgoing coefficients, for the th particle. Thus
(36) 
where denotes the truncated scattering matrix acting on the truncated expansion about the center of the th particle. If the particles were only translated, we would have for all , where is the truncated scattering matrix from the previous section with elements . We allow general rotations of particles; rotation of the th particle by angle introduces phase factors, so .
Multiparticle scattering has a key difference from single particle scattering, namely that the incoming field experienced by each particle consists of two parts: the (applied) incident field , and the contribution to the scattered field from all of the other particles. We denote by the translation operator (or M2L in FMM language [43]) that maps the outgoing coefficient vector from particle to their contribution to the local expansion coefficients centered at particle . With this operator in place, the incoming coefficients for the th particle are
(37) 
where is the (truncated) local expansion (35) of the incident wave relative to particle .
Combining (36) and (37), one can easily eliminate the incoming coefficients to obtain the following linear system that only involves the outgoing coefficients:
(38) 
where
The system (38) can be solved iteratively, using GMRES. Since each translation operator is dense, a naive matrixvector product requires operations, where is the order of the truncated expansion. The cost can be reduced to work by FMM acceleration, for which we refer the reader to [43, 17]. Furthermore, there exists an effective preconditioner for the system (38). Leftmultiplying by the block diagonal matrix results in the preconditioned system matrix . This significantly reduces the number of iterations.
The advantage of using the oneparticle scattering matrices over boundary integral equations is clear: the number of degrees of freedom per inclusion is only rather than the number of nodes needed to discretize the domain boundaries . For complicated inclusions, this permits a vast reduction in the number of degrees of freedom required, forming the basis for the socalled FMPS method [22]. Moreover, the blockdiagonal preconditioned multiple scattering equations are much better conditioned than the BIE (31)–(32), while FMM acceleration is particularly fast in this setting.
Remark 4.
It is straightforward to extend the method to more than one type of particles as long as the assumption that the enclosed circles are well separated still holds. The additional cost is simply the bookkeeping for the different scattering matrices.
5 Multiparticle scattering in the periodic layered medium
We now combine the schemes of the previous two sections. The field in the middle layer is the periodic layered contribution from (3), plus the scattered field from the inclusions and, cruicially, their neighboring nearfield phased copies. We need the notation for polar coordinates relative to the origin of the th particle translated by . Then,
(39) 
It only remains to set up the interactions between the layered periodic structure and the inclusion structure. We denote the translation matrix mapping layer densities and the middlelayer expansion to the incoming coefficients of all particles by , and the translation matrix mapping (phased summed) outgoing particle multipole coefficients to data on the layers and walls by . Adhering to the ordering of unknowns and conditions in (24), one can show that they have the forms (since the particles interact only with the middle layer),
(40)  
(41) 
where here denotes the blockwise transpose. , and map layer densities on , and the expansion to the incoming coefficients of the particles. , , and map the outgoing coefficients from particles to the values and normal derivatives on interfaces , and the discrepancy from to .
It is easy to construct the elements in , , and by direct evaluation of multipole expansions. Note that involves cancellations which mean that only particlewall interactions over distances greater than survive, as discussed in Remark 6.
To obtain the elements in , and , one again uses Graf’s addition theorem. In particular, the translation submatrix that maps the coefficients from one expansion to another expansion (the localtolocal or L2L operator in the FMM) is constructed through the following lemma.
Lemma 5.1 ([43]).
Let disk be centered at and let disk be centered at . Then the local expansion
(42) 
induces a field on disk of the form
(43) 
where
Let us denote by the multipole coefficients for all particles in . Combining the matrices (24), (38), (40) and (41), we obtain the final system:
(44) 
Since (44) is a rectangular illconditioned matrix, we cannot easily solve this whole system iteratively. However, we now present a Schur complement scheme to generate a smaller, wellconditioned square linear system for which an iterative solution is efficient. Since has size of order in both dimensions, it is much smaller than , and we can eliminate the unknowns via , where is the pseudoinverse of . We precompute