A fast and robust solver for the scattering from a layered periodic structure containing multi-particle inclusions

A fast and robust solver for the scattering from a layered periodic structure containing multi-particle inclusions


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 high-order 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 quasi-periodic Green’s function—fails at Wood’s anomalies. We instead use the free-space 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 single-inclusion scattering matrix, this is solved iteratively in time using a fast matrix-vector product. Numerical experiments show that a diffraction grating containing inclusions per period can be solved to 9-digit accuracy in under 5 minutes on a laptop.

Keywords: fast solver, periodic scattering, multi-particle 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, meta-materials, plasmonics, and other micro-scale structures, are becoming key to efficient devices, including lasers, sensors, anti-reflective surfaces and absorbers [21], and solar cells [3]. For instance, in thin-film 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 multi-particle wave scattering problems arise in acoustics and elastodynamics, and in general whenever a super-cell 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 frequency-domain 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 two-dimensional (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 high-order accuracy with appropriate surface quadratures. However, the resulting linear system is often dense, making a naive matrix-vector 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 free-space Helmholtz single- and double-layer 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 co-workers [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 sub-wavelength, 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 block-diagonal matrix and then solve via GMRES with FMM acceleration, with effort scaling linearly in . The result is a robust, efficient, high-order 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 multi-particle scattering and discusses the evaluation of the scattering matrix. The quasi-periodizing 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 plane-wave 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:


where , and where the wavenumber takes one of four values,

Figure 1: A 2D grating scattering geometry: a plane wave incident on a three-layered periodic structure with periodicity . A large number of identical dielectric obstacles are embedded in the middle layer. We use to denote the set of all these particles and to denote the th particle. The vertical dotted lines indicate the unit cell walls , , while the top and bottom dotted lines indicate the fictitious boundaries at . The three layers are denoted by , and .

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:


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 quasi-periodic if a function (such as ) satisfies


where is the incident horizontal wavevector. The factor

is the Bloch phase associated with translation by one period. Since is quasi-periodic, 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


which are in fact equivalent to (4) [11, Sec. 3].

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 square-root is taken as positive real or positive imaginary. Then the plane wave with wavevector is quasi-periodic 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:


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 half-space), or if thus (lower half-space). 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 , .

The following theorem [4, 6, 11] describes the well-posedness of the boundary value problem.

Theorem 2.1.

Fixing an incident angle , there exists a unique solution to (1)–(3) and (5)–(8) for all but a discrete set of .

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 unit-cell quasi-periodicity (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 non-periodic integral equation formulation into the correct periodic one is to use the quasi-periodic Green’s function [34], defined at wavenumber as


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 well-defined 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 inverse-square-root singularity) at Wood’s anomalies, causing a purely numerical breakdown in the solution of a what remains a well-posed problem. The key idea is that (9) can be rewritten as


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,


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 double-layer potentials [19], which we may define living on a general interface at wavenumber by,


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 ,


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 ,


Thus the single-layer potential is continuous for all , whereas the double-layer is generally discontinuous across .

Turning to the third layer , we similarly represent the scattered field using layer potentials on ,


where . The scattered field in the second layer has contribution from both and , thus


where .

To determine the unknown densities , , , and the coefficients , and , we enforce the following boundary conditions according to (3) and (5)–(8).

  • On and , the continuity condition (3) is imposed, giving,


    Note that since the representation of in the three layers involves layer potentials, the limits must be taken from the appropriate side of and using jump relations (15)–(18).

  • On and , where , the quasi-periodicity condition (4) is imposed:


    The left hand sides (phased differences) is sometimes known as the discrepancy [10].

  • On the parts of the artificial boundaries and lying in the unit cell (denote this part of the boundary by and ), the radiation conditions (7) and (8) are imposed for values and normal derivatives:

Figure 2: Scattering from a periodic structure without inclusions, with , , and . The interfaces are given by the graph , and the graph . (a) The singular value spectrum of matrix , vs index : original matrix (blue dots), vs after rescaling the columns of the blocks in corresponding to expansions (red dots). (b) Real part of the scattered field computed to 13-digit accuracy.

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 equally-spaced 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 16th-order 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,


Here the right-hand side vector has the form


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 free-space Green’s functions between source and target nodes, or or Rayleigh–Bloch expansions at target nodes.

  • and : Nyström self-interaction matrices for and respectively, including the phased summation over source near neighbors. For instance,


    (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 middle-layer -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 middle-layer -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 lower-right block of in (24) involves only the effect of the auxiliary periodizing degrees of freedom ( and Rayleigh–Bloch expansions) on the auxiliary matching conditions (discrepancies and radiation conditions). In related work this block is given the symbol [10, 18].

The matrix is generally rectangular, depending on the specific numbers of discretization nodes. Although the upper-left 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 ill-conditioned. 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 ill-conditioning 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 small-norm least-squares 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 Multi-particle scattering

Multi-particle 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


in the interior, and scattered wave


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


Here and are unknown single- and double-layer 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:


Let and denote the solution to (31)–(32) for . We may then precompute the scattering matrix elements as the multipole expansion coefficients (truncated up to terms)


from the densities and via Graf’s addition theorem, giving the standard formula [43, 17]

Figure 3: Two inclusions and their enclosing disks. The scattering matrix for each inclusion with wavenumber is defined as the map from an incoming field on to the corresponding outgoing field.

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


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


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 .

Multi-particle 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


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:



The system (38) can be solved iteratively, using GMRES. Since each translation operator is dense, a naive matrix-vector 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). Left-multiplying by the block diagonal matrix results in the preconditioned system matrix . This significantly reduces the number of iterations.

The advantage of using the one-particle 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 so-called FMPS method [22]. Moreover, the block-diagonal 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 Multi-particle 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 near-field phased copies. We need the notation for polar coordinates relative to the origin of the th particle translated by . Then,


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 middle-layer 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),


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 particle-wall 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 local-to-local 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


induces a field on disk of the form



Let us denote by the multipole coefficients for all particles in . Combining the matrices (24), (38), (40) and (41), we obtain the final system:


where , as in (38), and is the right-hand side vector (25).

Since (44) is a rectangular ill-conditioned matrix, we cannot easily solve this whole system iteratively. However, we now present a Schur complement scheme to generate a smaller, well-conditioned 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