A 3D spatial spectral integral equation method for electromagnetic scattering from finite objects in a layered medium

The generalization of a two-dimensional spatial spectral volume integral equation to a three-dimensional spatial spectral integral equation formulation for electromagnetic scattering from dielectric objects in a stratified dielectric medium is explained. In the spectral domain, the Green function, contrast current density, and scattered electric field are represented on a complex integration manifold that evades the poles and branch cuts that are present in the Green function. In the spatial domain, the field-material interactions are reformulated by a normal-vector field approach, which obeys the Li factorization rules. Numerical evidence is shown that the computation time of this method scales as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(N \log N)$$\end{document}O(NlogN) on the number of unknowns. The accuracy of the method for three numerical examples is compared to a finite element method reference.


3
206 Page 2 of 22 structures, where the number of unknowns is often very large, there is a demand for fast solvers, for which the computational complexity scales well for large numbers of unknowns. A good strategy to find a potentially efficient algorithm is to exploit symmetries. For stratified media such a symmetry is the translation symmetry in the layered background medium. This symmetry can be exploited through the use of a Fourier representation in a volume integral formulation.
In a stratified dielectric medium, an analytic expression exists for the Green function in the electromagnetic case, as a function of one spatial coordinate in the direction of stratification and two spectral coordinates in the two directions perpendicular to the stratification (Kong 1975;Felsen and Marcuvitz 1973;Chew 1995;Michalski and Mosig 1997). It is advantageous to use the stratified-medium Green function, since it incorporates the response of the multilayer medium analytically. Therefore, little computation time or memory is used for computing the scattered electromagnetic field throughout the entire layered stack, since the electric field on a domain slightly larger than the scattering object suffices. It is possible, using Sommerfeld (1909) or Fourier integrals, to transform the Green function completely to the spatial domain and then use it in an integral equation method (Chew 1995, Chapter 8;Felsen and Marcuvitz 1973, Chapter 5;Kong 1975, Chapter 4;Wait 1970, Chapter 2). However, these Sommerfeld integrals are often tedious to compute, because of poles and branch cuts present in the Green function that can be located on or close to the integration path. Since the Green function has to be re-calculated for every modification in the multilayer medium, caching the Green function in a library is only advantageous when exactly the same multilayered medium is used many times.
It is also possible to use the Green function directly in the spectral domain, where it is known analytically. For a periodically repeating object, the Green function decomposes into a discrete set of modes as derived in for example (Beurden 2011(Beurden , 2012. Problems with poles and oscillations along branch cuts in the Green function (Chew 1995;Felsen and Marcuvitz 1973) can be avoided on such a discrete set of modes since the modes and locations of the poles will most likely not coincide. However, for a finite scatterer the spectral domain is continuous and now the poles and oscillations along the branch cuts are hard to discretize Beurden 2016, 2017). Deformations of the Sommerfeld integration path to a complex-plane path (Ruiter 1981;Newman and Forrai 1987;Hochman and Leviatan 2010;Michalski and Mosig 2016) can help to evade these poles and branch cuts. In  an algorithm for two-dimensional electromagnetic scattering with TE polarization in a multilayered medium is presented, where both contrast-current density and scattered field are represented on a path in the complex plane of the spectral domain. It is this path that allows for the use of Gohberg and Koltracht (1985) fast, flexible and recursive Green-function convolution in the stratification direction.
The first challenge in three dimensions is that, instead of one, now two directions perpendicular to the stratification direction need to be handled. The complex integration path is turned into a complex integration manifold and since the transformation from the spatial domain to the complex integration manifold is part of the core of the algorithm, transformations back and forth need to be computationally efficient. We show an integration plane consisting of nine regions of three distinct types and show transformations to and from the spatial domain that can be computed in O(N log N) time, where N is the number of spectral unknowns.
The second challenge is that the discontinuity of both the permittivity and the electric field at material interfaces leads to poor convergence in spectral formulations (Li and Haggans 1993). This effect was also observed for a Gabor-frame based solver for TM-polarized scattering . For periodic scattering problems with a discrete spectral expansion a Page 3 of 22 206 reformulation of the field-material interactions corrects this poor rate of convergence (Granet and Guizal 1996;Lalanne and Morris 1996), which is explained in more detail in Li (1996) introducing the so-called Li-rules. In Dilz et al. (2017) it is shown that the same mechanism can also be used for a continuous spectral expansion and the algorithm of  is extended to efficiently deal with the discontinuous field-material interaction in a way that does abide these Li-rules. Here, we propose a generalization of this method to three dimensions. Inspired by Beurden and Setija (2017), we show that a normal-vector field formulation (Popov and Nevire 2001) can be used for three-dimensional scattering to replace the field-material interaction.
We start by a short formulation of the volume integral equation. Subsequently, we give a more detailed explanation of the discretization, with emphasis on the complex-plane spectral domain representation, followed by a short summary of the normal-vector field framework. The applicability of the present algorithm is highlighted by three numerical examples, with numerical evidence that the computation time scales as O(N log N) with the number of unknowns and comparison against a finite-element reference calculation.

The volume integral equation
Consider a stratified dielectric medium where layers with different relative permittivities are stacked in the z-direction. Layer n is located between z n and z n+1 and has relative permittivity rb,n . Index n = 0 coincides with the top half-space, z < 0 , and index n = N L with the half-space z > z N L +1 below all layers, an example of which is also illustrated in Fig. 1 object is characterised by a relative permittivity function r ( ) , with = (x, y, z) , or more conveniently by the contrast function which is nonzero only in the object. An incident electromagnetic field originates from the upper half-space at arbitrary angles and polarization. The electric field in presence of the multilayered background medium rb,n , but in absence of the dielectric object can be readily calculated (Chew 1995;van Kraaij 2011) and is denoted as i ( ) . The dielectric object generates a scattered field s ( ) that together with the incident field i ( ) adds up to the total electric field ( ) , i.e.
The scattered field s ( ) can be calculated via the multilayer Green tensor through where the contrast current density ( ) is given by the field-material interaction which is again nonzero only in the scattering object. Combining Eqs.
(2)-(4) yields the integral equation that we propose to solve However, for an efficient numerical scheme several refinements have to be made.

The Green function
The three-dimensional integral in Eq. (3) yields, when implemented naively, an O(N 2 ) matrix-vector product, with N the total number of unknowns and by employing an iterative solver. Analoguous to Beurden (2016, 2017), , we represent the Green function, the contrast current density , and scattered field s in the spectral domain in the transverse xy-plane. We denote coordinates in the transverse plane as T = (x, y) , and in the spectral domain as T = (k x , k y ) . We use a Fourier transformation defined as where we distinguish functions in the spectral domain by arguments containing k x , k y and T and in the spatial domain by the arguments x and y and T . In the spectral domain, a spatial convolution can be executed with O(N x N y ) complexity, with N the number of unknowns used in direction . The transverse convolution in Eq. (5) can be carried out efficiently. The remaining integration in the z-direction can be calculated in O(N z ) time via the recursive algorithm proposed by Gohberg and Koltracht (1985). The multilayer Green tensor in Eq. (3), can be separated in a homogeneous-medium part yielding s,h and reflected waves moving up, u , and down, d . The homogeneousmedium part of the scattered field is given by where the homogeneous-medium Green tensor is given in Cartesian components (x, y, z), respectively, as Here, k 0 is the wave number k 0 = √ 0 0 and is defined as = Note that the factor exp(− |z − z � |) propagates the electric field over a distance |z − z � | , and will therefore be referred to as the propagation function. Now the scattered field s can found by adding reflected waves u∕d from the layer interfaces to the homogeneous scattered field s,h , where u and d refer to waves moving up or down respectively. Consequently, we have with ℛ , ( T , z) the three-dimensional effective reflection coefficient that contains both h and e polarization, which can be calculated from the effective reflection coefficients for h polarization r , h ( T ) , and for e polarization r , e ( T )  as This matrix projects the e and h polarized parts of the electric field to effective transmission coefficients r , e and r , h , respectively. The definition of these effective reflection coefficients is given in , , which is based on the expositions about multilayer media in Kong (1975, Chapter 4), (Wait 1970, Chapter 2), van Kraaij (2011. Since the field-material interaction in Eq. (4) is calculated in the spatial domain and the Green-function operation in Eq. (9) in the spectral domain, we need a fast and efficient means of transforming the current density ( T , z) to the spectral domain and the scattered field s ( T , z) back to the spatial domain. We propose to use a two-dimensional Gabor-frame in the transverse plane, since a Gabor frame is efficient to represent the operation of Fourier transformation. It can be represented analytically by a mere transposition of the coefficient matrix in O(N) operations (Dilz and Beurden 2016).

The Gabor frame
We use a Gabor frame with Gaussian window function with width X in the x-direction and Y in the y-direction. This defines the oversampled twodimensional Gabor frame as with two-dimensional indices = (m x , m y ) and = (n x , n y ) . Here, the spectral spacing is K x = 2 ∕X and K y = 2 ∕Y and oversampling x x < 1 and y y < 1 . The number of coefficients in and is allowed for to be different in the two directions. Gabor coefficients can be calculated as with dual frame There is freedom of choice for the dual window function (x) , but we choose the dual frame function calculated via the Moore Penrose inverse (Feichtinger and Strohmer 1998;Bastiaans 1995), since it is widely used and exhibits a convenient exponential decay in both the spatial and spectral domain. We use the Fourier transformation of Eq. (12) to discretize functions in the spectral domain. This has the advantage that the operation of Fourier transformation reduces to merely a tranposition of coefficients. Details on operations such as Fourier transformation and multiplication of Gabor-represented functions can be found in Dilz and Beurden (2016) for one dimension and the generalization to two dimensions is straightforward.

A complex-plane deformation of the integration manifold
In the z-direction, the integration with the Green tensor in Eq. (7) is discretized completely in the spatial domain. Since it was shown that a piecewise-linear approximation in the z-direction is effective Beurden 2016, 2017;, we propose to use it here again. In the z-direction, the basis functions are then defined as with z the discretization step in the z-direction.
For the discretization in the xy plane, a method similar to the two-dimensional cases in ,  is proposed. The Green function contains poles due to the effective reflection coefficients and many oscillations along the branch cuts may occur. Both these poles and oscillations cannot be represented efficiently in a Gabor frame representation. In the two-dimensional case, these problems can be circumvented by representing the Green-function in the transverse direction in Eq. (9) on a path in the spectral complex plane. For three-dimensional problems, this path can be generalized to a two-dimensional integration manifold in the transverse T coordinates on which the transformation back to the spatial domain takes place. In the k x -direction, the complex spectral path is defined by the function x (k x ) , with k x ∈ ℝ and x ∈ ℂ as and similarly for y (k y ) with A y defining the path along the k y -direction. Here, A x and A y are constants that can be chosen individually. Numerical experiments show that a choice such that A x W x and A y W y are in the range 2 … 5 , yields optimal accuracy, with W x and W y as in Fig. 1. With the coordinate change from T to T , Eqs. (7) and (9) contain smooth functions and these can be used in combination with the Gabor-frame discretization. This complex spectral manifold divides the complex T domain into nine regions as depicted in Fig. 2. All functions in the spectral domain will be represented on this T manifold. Using Jordan's lemma, the Fourier transformation to the spatial domain can be carried out over the T manifold. Closing the contour at T → ∞ is not needed, since the representation using Gabor frames converges to zero rapidly.
It should be noted that the complex integration path of (16) is not the only possible choice. For example, in Dilz and Beurden (2018) a different continuous path in one dimension is chosen. The current choice for the path in (16) was made because the large horizontal stretches allow for fast transformations to and from the spectral domain. Different choices might not allow for such computational efficiency.

Discretization in regions of type 1
Most information is contained in regions of Type 1, since A x and A y are relatively small compared to the complete spectral range to be discretized. The contrast current density is transformed to the complex spectral integration manifold via for the northeast (NE) quadrant, i.e. k x ≥ A x ∧ k y ≥ A y , and similarly for the other regions of Type 1. Analogously, the transformation of the scattered field back to the spatial domain is obtained as with the cut-off function c NE ( T ) equalling 1 on the NE region and zero elsewhere. The Fourier transformation can be performed in O(N) operations and the operation of multipication in O(N log N) operations, for functions represented by N Gabor coefficients. Therefore, the total of these operations allows for an O(N log N) computational complexity. All this means that the scattered electric field s is represented by a five-dimensional array of coefficients s,NE , ,l , with m x , n x and m y , n y corresponding to the Gabor frame on the coordinates, k x + jA x and k y + jA y respectively. The index corresponds to a piecewise-linear (PWL) representation in the z-direction. The scattered electric field in region NE is then approximated as The Green function consists of several parts, some of which are depending on the complex propagation constant ( T ) = √ − rb,i k 2 0 + 2 T . On the real k x k y -plane ( T ) touches, but does not cross, two branch cuts at t = (0, 0) in the case of lossless media. For lossy media the branchcuts are located at some distance from the origin. For both cases, the path passes just in between these two branchcuts. However, when a Type-1 region such as the NE-region is continued to the complete (k x + jA x , k y + jA y ) plane, a branch cut is crossed just outside the NE region, as illustrated in Fig. 3. The branch cut is located on a straight line through T = (0 + jA x , 0 + jA y ) and direction of the line depends on the choice of A x and A y . The continuous nature of a Gabor-frame representation does not allow for an abrupt stop of the discretization domain at the borders of a Type-1 region. Therefore, such a Gabor-frame representation of the Green function exhibits significant Gibbs ringing from the branch cut that spreads into the Type-1 regions. For a two-dimensional case, this is described in , where a linear continuation of the Green function is proposed that suppresses strong Gibbs ringing.
In three dimensions, this issue can also be resolved by making a first-order continuation of the functions to eliminate the branch cut. Since the branch cut can be located close to the k x = 0 or k y = 0 axes, and the function values are needed at k x > A x and k y > A y , we start the continuation of the functions in the middle at k c x = A x ∕2 and k c y = A y ∕2 . Then the Gibbs phenomenon from the discontinuous second derivative will be at a short distance from k x = A x and k y = A y , where Region NE begins. For the continuation of a function f (k x , k y ) along the k y -axis we choose for k x < A x ∕2 and k y > A y ∕2 and similarly f c,y (k x , k y ) can be constructed for k x > A x ∕2 and k y < A y ∕2 , which is illustrated in Fig. 3. The Gaussian factor is added to make the continuation decay to zero slowly. The third part, Note that this expression equals the expression obtained from continuing f c,x onto this domain. The derivative of the function f is calculated using a forward finite-difference method, with a difference of 10 −4 A x or 10 −4 A y for the x and y direction, respectively. For most functions, = min(X 2 , Y 2 ) is a good choice. However, for one part of the Green function, notably the propagation function e − z , care has to be taken that its absolute value does not exceed one in the continuation. By increasing the value of , this condition can always be satisfied. More details can be found in .
A general remark about the importance of this continuation is in place. In principle, the Gibbs phenomenon in a Gabor frame representation is not of much significance, unless two functions with discontinuities at the same position are multiplied. The Li-rules (1996) state that when two functions with spatial discontinuities at the same position are multiplied to form a convolution in the spectral domain, the convergence of this convolution is poor. The Li-rules also apply to Gabor frames ) and since the spatial and spectral domain are both represented by a Gabor frame, a spatial version of the Li-rules is also applicable to the Gabor frame. These spatial Li-rules state that when two spectral functions with discontinuities are multiplied in a Gabor-frame representation, a poor convergence is observed. Now when the NE region of the electric field with its branch cut is multiplied by the cut-off function c NE ( T ) in Eq. (18), which is discontinuous at k x = A x and at k y = A y , functions are multiplied for which the locations of discontinuities almost touch each other. This leads to near-violation of the spatial Li-rules. Since the discontinuities are not exactly at the same location, a high sampling would in principle solve this issue. However, this would require an excessive sample density that is avoided by the continuation of the Green function parts proposed in Eqs. (20) and (21).

Discretization in regions of type 2
First we will approximate the contrast current density in k x around k x = 0 with a Taylor expansion that is found through a Vandermonde matrix. This Taylor expansion is then applied to find corresponding values of the contrast current density on the line Im( x ) = Re( x ) , on which a PWL basis is used as a discretization. This PWL basis consists of 2N s + 1 sampling points, p(1 + j)A∕N s , with p ∈ {−N s , … , N s } . Afterwards, we give a means to directly Fourier transform from the discretized N region to spatial-domain Gabor coefficients. We will only consider the northern (N) region of the complex spectral integration manifold since the E, W, and S region follow by analogy.
For the calculation of the current density in the N region, function values of N are available at the lines x = k x ± jA x , which were calculated via the Gabor representation in the NE and NW region. The analyticity of allows to produce a Taylor expansion of around k x = 0 from values at the lines Im( x ) = ±A x . Afterwards, this Taylor expansion is used to calculate values of the contrast current density at the line Re( x ) = Im( x ) , where they are needed for discretization in the N region, as is shown in Fig. 4. Close to k x = 0 , N can be approximated as where n (k y , z) = n k x N (k x , k y , z) , and 4N v + 2 is the total number of terms in this Taylor expansion. Values for N (k x ± jA x , k y ) can be obtained from the results for the NE and NW region, by using a fast Gabor transform Dilz and Beurden (2016), Bastiaans (1995) that yields values at k x = n k x ± jA x for n ∈ ℤ , with k x the spectral sample spacing corresponding to the Gabor frame. Values for n can be found by solving a small Vandermonde system (Press et al. 2007, Chapter 2.8). By constructing the vector k of k x values as Vandermonde system can be written as a matrix equation K ⋅ = (k y , z) = N (k, k y , z) . The element K mn of the m'th row and n'th column of matrix K is given by K mn = (k m ) n , the n-th power of element m in k . We solve this system using the inverse of K , i.e. Now that it is possible to express the Taylor coefficients n in terms of the 2N + 1 samples on the NW-region, i.e. (k x − jA x , k y + jA y , z) , and the 2N + 1 samples in the NE-region, i.e. (k x + jA x , k y + jA y , z) , they can be used to evaluate the Taylor expansion in Eq. (22) on the N-region, where Im( x ) = Re( x ) . We will write this as a matrix-vector product using the matrix T . The matrix T transforms from a Taylor series to an equidistant sampling on We use the array of numbers N p,m y ,n y , to represent the current in the N region of the complex integration domain. Index p ∈ {−N s , N s } points to the set of piecewise-linear basis functions that are used in the k x -direction on the line x ((1 + j)pA∕N s ) . In the k y -direction we use a Gabor frame, denoted here by indices m y and n y , therefore the y dependence in Eq. (24) is replaced by this set of Gabor indices. Again, a set of N z PWL functions is used in the z-direction denoted by the index .
Having dealt with the transformation to the N region, we will now deal with the transformation from the N region back to its spatial-domain counterpart. After multiplication of the contrast current density N p,m y ,n y , with the Green function (see Sect. 3.1), the contribution of the North part of the scattered electric field yields s,N p,m y ,n y , . From this array we can make an approximation on the N region of the scattered electric field Here s x ,p are piecewise-linear (PWL) basis functions To transform the N region of the scattered electric field s,N p,m y ,n y , back to the spatial domain where it is discretized in the Gabor frame with coefficients s,N , ,l , we use the Fourier transforms of the PWL functions in Eqs. (26) and (25), i.e.
Since the x direction is discretized using Gabor coefficients in the spatial domain, the x-dependence of this function I N p (x) must be Gabor-transformed into Gabor coefficients I N p,m x ,n x . These Gabor coefficients are calculated during initialization of the algorithm via e.g. Eq. (13) or a fast Gabor transform. Now the contribution of the N region to the scattered electric field in the spatial domain is given by where the dots indicate the contributions from the other eight regions to the scattered field.
Similar to regions of Type 1, some parts of the Green function are discretized using a continuation such as in Eq. (20), to avoid a branch cut. For example, for the N region the continuation is only needed in the y-direction, since a Gabor frame is employed in this direction only and a PWL discretization does not suffer from Gibbs ringing. The construction for a one-dimensional continuation is described in more detail in Dilz and Beurden (2017).

Discretization in the region of type 3
For the middle (M) region, a two-dimensional version of the construction for the N region is used. Since the generalization is fairly straightforward, we will not write it down explicitly. The only difference here is that we use a total number of 2N m + 1 PWL functions per direction. We use a different number of PWL functions in this region since, depending on the simulation parameters, the accuracy can depend significantly on the choice for N m . Since the middle part contains information of waves with small T , it contains information about waves traveling almost parallel to the z-direction. Especially for scatterers that are larger in the z-direction, a larger N m is required.
An important remark on the use of Vandermonde matrices is that they are generally illconditioned when a uniform sampling is used, such as is the case in the NE and NW regions. In principle, this could lead to a poor conditioning of the K matrix and therefore to an unstable inverse when the matrix is increased in size. However, the amount of information on the interval T ∈ [−A x , A x ] × [−A y , A y ] is so small that large matrices are not needed.
There are two reasons that a relatively large number of PWL basis functions (typically N m > 10 and N s > 10 ) is needed in regions of Type 2 and 3. The first is that a PWL basis is relatively inefficient compared to a Gabor frame. For the second reason we have to look at both the spatial and the spectral domain. Since the contrast current density is confined to a finite region only, its Fourier transform is fairly smooth. However, the scattered electric field s is not confined to the simulation domain, and therefore its Fourier transform is much less smooth. On the Type-1 and Type-2 regions this lack of smoothness is compensated by a representation on complex spectral coordinates , where the Green function is much smoother. However, the Type-3 region is not shifted as far into the complex plane as the Type 1 and Type 2 regions, and therefore the Green function is less smooth in this region. Since the Green function is implemented recursively, for intermediate results, i.e. the scattered field in between z min and z max , this lack of smoothness should be represented accurately. Afterwards, when the transformation to the spatial domain is performed, this roughness on the M region corresponds to contributions outside the simulation domain, but ignoring the roughness is not an option since it leads to accumulating errors in the recursive handling of the Green function. This is especially important when z max − z min is large compared to the wavelength.

Correspondence between simulation parameters and accuracy
Since there are many simulation parameters, it is not trivial to find values for these parameters that produce both a good accuracy and short computation time. This list is intended to clarify which simulation parameters influence which part of the algorithm. This list is intended as a general guideline for optimal results.
1. Start with a Gabor frame with X = Y = , the wavelength of the light-source, and = = √ 2∕3. 2. C h o o s e m x > 3 + W x ∕ X a n d s i m i l a r l y m y > 3 + W y ∕ Y . C h o o s e n x = n y > 5 max x∈ (1 + (x)) , which guarantees at least 11 unknowns per wavelength per direction. Test whether a function (e.g. a Gaussian with width X), can be represented with the required accuracy over the entire simulation domain . When the accuracy is too low everywhere, increase , when the accuracy is too low at the boundary of only, increase . 3. Start at A x = 3∕W x and A y = 3∕W y , N v = 1 , N s = 10 and N m = 10 . Test whether a set of spatially and spectrally localized functions (e.g. modulated Gaussians that are shifted along the entire simulation domain) can be transformed to the complex spectral integration manifold and back again with the required accuracy. Note that the exponential function in Eq. (17) reduces the accuracy of the Gabor frame. Therefore, a simultaneous decrease of A x and increase of improves the accuracy in the transformation between the spatial and spectral domain. Especially when a high accuracy is needed, N might need to be increased when the error in the N, E, S, W and M regions is too large. Also an increase in N s and N m can be considered when the error in the PWL interpolation is found to be too large. 4. The Green function (Eq. 8) contains a factor −1 , that has a strongly peaked behavior around � T � = √ rb,i k 0 . Test whether the function −1 is represented accurately enough by the Gabor frame with the current parameters. Otherwise A x can be increased (which decreases the accuracy in the previous step) or can be increased (which increases the computation time, but leaves the accuracy in the previous step invariant). 5. Especially for large z max − z min , a lot of information is stored in the M region, which contains information about waves traveling in a narrow cone around the z axis. The main culprit here is the function e − (z max −z min ) in Eq. (8). Choose N m such that this function can be well approximated in the M region.
6. Another simulation can be run with lower values. When both results agree well, convergence in the spectral range has been reached, otherwise should be increased for higher accuracy.

Efficient field-material interaction
Formulating the field-material interaction as proposed in Eq. (4) yields poor convergence since it violates the Li-rules (1996). We propose to use a normal-vector field approach (Popov and Nevire 2001;Beurden and Setija 2017). In  it is shown that when the Li-rules are satisfied, good convergence is reached in a continuous spectral discretization in a formulation similar to the RCWA formulation by Granet and Guizal (1996). We follow the same approach as Popov and Nevire (2001), Beurden (2011, Beurden and Setija (2017), van Beurden (2013) in constructing normal-vector fields and the following is intended as a short summary of that method. When the permittivity is discontinuous at a material interface, it is observed that the electric field normal to the surface is discontinuous, but the electric flux density normal to such a surface is continuous. Therefore, in the field-material interaction in Eq. (4), both and the normal component of are discontinuous and multiplication of two discontinuous functions represented by Gabor shows poor convergence , since it violates the Li-rules (1996). An auxiliary field is introduced that is composed of in the direction normal to every surface of discontinuous and parallel to each of those surfaces. Since this fixes the choice of only at the boundaries of dielectric objects, there is much freedom in choosing it away from the interfaces. Normal-vector fields (Popov and Nevire 2001;Rafler et al. 2008;Götz et al. 2008;Beurden and Setija 2017) can be a good tool to systematically construct an auxiliary field .
Since we use Gabor coefficients only in the transverse plane, we apply the normal-vector field formulation only in the transverse plane. For objects with interfaces that are not aligned with the z or transverse plane, a staircasing approximation is needed. When T ( ) is a vector field of unit amplitude that is directed normal to the transverse part of all discontinuous surfaces in and when ( ) is a scalar function that equals one at these discontinuities, these functions can be used to construct the desired auxiliary field as The field-material interaction in Eq. (4) can be re-written as and the electric field can be recovered from where the Cartesian component i of the electric field due to the Cartesian component j of auxiliary field F is calculated by employing operator C defined as with ij denoting the Kronecker delta and similarly Using this normal-vector field approach, there is a great freedom for the shape of the scatterer. Depending on the scatterer shape, a suitable T and must be chosen. It is possible to choose T and for basic geometric elements, such as rectangular blocks, triangular prisms, or circular cylinders. More intricate objects can then be constructed from

Numerical results
We have chosen three testcases to validate the present algorithm. As a reference to validate our results we use the commercial FEM code JCMWave (Burger et al. 2013). Our goal is to achieve an accuracy of 10 −3 , which is sufficient in our applications, e.g. due to measurement noise or fabrication tolerances. To achieve high accuracy in the validation, we use a relatively small, low-contrast scatterer in the first testcase. A small dielectric cube is embedded in a dielectric medium as shown in Fig. 5a, together with the remaining details of the setup. The incident wave is characterized by a Cartesian wave vector with components = (−k 0 sin(70 • ), 0, k 0 cos(70 • )) , with the electric field polarized in the xz-plane and with unit amplitude. We choose a Gabor frame with X = Y = 80 nm, = = √ 2∕3 . For the highest accuracy we use a Gabor frame, Eq. (12), restricted to m x , m y ∈ {−7, … , 7} and n x , n y ∈ {−10, … , 10} , which equals one basis function per 3.1 nm. In the z-direction we use a step size of 2.5 nm. For the sampling of the regions of Type 2 and 3 we use N = 2 , N s = N m = 15.
With these simulation parameters, the simulation domain in the xy-plane extends over a larger region than the scatterer itself, as is visible in Fig. 6a. In this figure, the norm of is shown on the plane z = 100 nm. In Fig. 6b the absolute difference between results from the present algorithm and the JCMWave reference are shown. Over large regions of the simulation domain the absolute difference is smaller than 10 −5 . From y = −50 nm to y = 50 nm and close to the edges of the cube agreement between both simulations is not as good as on the rest of the simulation domain. This is caused by Gibbs-ringing at the discontinuities in the electric field, especially in the x-component. The field-material interaction operators Fig. 7 The far field for the case in Fig. 5a as a function of the transverse wavenumber T ∕k 0 , scattered back into the half-space z < 0 . In a the modulus | s | of the scattered electric field is shown. In b the difference between a JCMWave validation run and the present algorithm is shown. An average relative error of 4 × 10 −5 was observed. Since an interpolation of the reference data is used that is not accurate at the edge of the radiation circle, the far field data is truncated for large T C and C of Sect. 5 exhibit this Gibbs ringing because they are truncated in the spectral domain. This Gibbs ringing is therefore also present in the scattered electric field.
Since this Gibbs ringing has a very small spatial period, it does not radiate into the far field. Because the far field is the most interesting for the application, we use the far field as a reference for the accuracy of the method. As can be observed in Fig. 7, the error in the far field is much lower than in the near field, because the absence of the Gibbs ringing. The average relative difference with an ℒ 2 -norm in the far field data equals 4 × 10 −5 . Clearly, the far-field results agree much better than the near-field results. The small size of the scatterer and its low contrast results in a far field pattern that does not vary much with the angle. This example is therefore somewhat uninteresting, however, it has the advantage that the FEM reference could achieve a high accuracy in a multilayered scattering problem.
In Fig. 8, we show how both the accuracy and the computation time scale with the number of unknowns used in the calculation. The horizontal axis in Fig. 8a contains the sample density, which was lowered by decreasing the range of the -index in Eq. (12), where n x , n y ∈ {−r, … , r} from r = 10 down to r = 1 . This corresponds to a sample density 1∕ x = 1∕ y = (2r + 1)∕ √ ∕ X . The other simulation parameters were kept (a) (b) Fig. 8 In Figure a both the computation time and the relative error in the far field, computed as the average of Fig. 7b for a range of n x and n y for the Gabor frame. In b the same is shown, but now for different sampling z in the z-direction . Since FFTs are used, we expect an O(N x N y log(N x N y )) behaviour, but the logarithms are apparently negligible compared to other parts of the algorithm at this simulation size. Fig. 8b shows a clear O(N z ) behaviour, which is expected from Gohberg and Koltracht's recursion (1985). The second example for which we provide computational data consists of a dielectric cylinder embedded in a multilayered medium as is described in Fig. 5b. In Fig. 9, the electric field is shown at z = 10 nm for X = Y = 100 nm, = = √ 2∕3 and m x , m y ∈ {−4, … , 4} and n x , n y ∈ {−7, … , 7} , which equals one basis function per 5.1 nm. In the z-direction, 41 basis functions are used with stepsize z = 2.5 nm. For the sampling of the regions of Type 2 and 3 we use N = 2 , N s = N m = 15 . From Fig. 9b it is clear that the difference in the result from the present algorithm compared to the JCMWave reference is somewhat larger than for the previous case. However, this is due to a less accurate reference that was calculated with a lower order p-refinement. A higher order p-refinement was not feasible with the available computational resources.
In Fig. 10, the absolute value of the far field reflected into the upper half-space is plotted. The relative error for the simulation in the far field is 2.8 × 10 −3 , which is significantly larger than for the cubic scatterer, because of the reference, with lower accuracy.
We have calculated the far field with a finer discretization, to show the convergence of the algorithm. We emphasize that the present algorithm does not excel at small computational domains. However, for a small computational domain a more accurate validation result was feasible than for a very large computational domain. The reason that the present algorithm is relatively slow for small simulation domains is that a lower limit on the number of unknowns in the x and y-direction exist, since the Gabor frame (as we choose it) is inaccurate over a range of at least three window widths X (see Eq. (12)) distance from the point at which the Gabor frame is truncated to zero. Therefore, at least seven windows are needed for an accuracy of 10 −3 in the middle of the computational domain, both for the spatial index m x as well as for the spectral index n x in Eq. (12). Consequently, at least seven values for index m x and seven values for index n x are needed, which amounts to a total of 49 coefficients per direction at minimum. Since Fig. 10 The far field for the case in Fig. 5b as a function of the transverse wavenumber T ∕k 0 , scattered back into the half-space z < 0 . In a the modulus | s | of the scattered electric field is shown. In b the difference between a JCMWave validation run and the present algorithm is shown. An average relative error of 2.8 × 10 −3 was observed  Fig. 12 The far field for the case in Fig. 5c as a function of the transverse wavenumber T ∕k 0 , scattered back into the half-space z < 0 . In a the modulus | s | of the scattered electric field is shown. In b the difference between a JCMWave validation run and the present algorithm is shown. An average relative error of 2.5 × 10 −3 was observed we use a Gabor frame in two dimensions, a minimum of 2401 unknowns is needed for a minimum simulation domain. Note that we assumed the use of a Gabor frame with = = √ 2∕3 and the Moore-Penrose inverse to calculate the dual window, where the accuracy increases exponentially with the distance to the truncated coefficients for this choice (Böcskei and Janssen 2003) for sufficiently smooth functions. This effect only exists for small simulation domains. For large simulation domains, the number of unknowns at the edge of the simulation domain is negligible.
The third and final example for which we computed the scattered electric field consists of six dielectric blocks of 350 × 500 × 100 nm deposited on a slightly lossy dielectric substrate as is shown in Fig. 5c. In Fig. 11, the electric field is shown at z = 10 nm for X = Y = 500 nm, = = √ 2∕3 and m x ∈ {−7, … , 7} , m y ∈ {−4, … , 4} and n x , n y ∈ {−7, … , 7} , which equals one basis function per 27 nm. In the z-direction, 21 basis functions are used with stepsize z = 5 nm. For the sampling of the regions of Type 2 and 3 we use N = 2 , N s = N m = 20 . Since the scatterer is larger than in the previous examples, it was efficient to choose larger window widths X and Y, which results a coarser sampling. From Fig. 11b it is clear that this coarser sampling generates a somewhat more pronounced Gibbs-ringing from the edges. However, in the far field, which is shown in Fig. 12, the average relative difference with the JCMWave reference calculation is similar to that for the cylinder case, i.e. 2.5 × 10 −3 , where the estimated relative accuracy of the reference calculation was on the order 2 × 10 −3 . Even though the scatterer extends much wider in the xy plane, the number of unknonws in the xy-plane was increased by only a factor 5 / 3, while the accuracy in the far field remained similar. This clearly shows that the present method performs better for scatterers larger than a wavelength in size.

Conclusion
A volume integral equation for 3D scattering from finite dielectric objects embedded in a dielectric layered medium was presented in the mixed spatial and spectral domain and an algorithm based on Gabor frames was presented for the discretization. The algorithm employs a mixed spatial spectral formulation and Gabor frames for the discretization. A representation of the Green function, contrast current density, and scattered electric field on a complex integration manifold is employed in the spectral domain. A normal-vector field formulation in the transverse spatial domain is employed to improve the convergence in the field-material interaction.
The accuracy of the present algorithm was compared to a FEM algorithm. The results of both algorithms in the far field agree with each other up to the 4 × 10 −5 in one small numerical example. In the other two examples an agreement up to 2.8 × 10 −3 and 2.5 × 10 −3 were observed, because the FEM algorithm did not fully converge with the computational resources at hand. Numerical evidence was presented that the computational complexity of the present algorithm scales as O(N log N) with the number of unknowns.