Finite element based Green&#x02019;s function integral equation for modelling light scattering

Wen Li; Dong Tan; Jing Xu; Shubo Wang; Yuntian Chen

doi:10.1364/OE.27.016047

1. Introduction

Green’s function integral equation (GIE) is often used to calculate scattering of light [1–9]. The boundary conditions of GIE are satisfied directly and the unknowns are only distributed in the region of the scattering object, which means that only the scatterer is needed to be discretized and field values at all the other locations can be obtained by overlapping the integration between the Green’s function and the field inside the scatterer [10–23]. This is an appealing advantage compared with other methods, e.g., finite element method (FEM) and finite difference time-domain (FDTD) method, which requires to calculate the region inside the scatterers as well as outside the scatterer. The introduction of FEM techniques, e.g., advanced discretization strategies [24–26] and quadrature-rule [4, 27], not only yields an improved geometric approximation but also improves the numerical accuracy of the integral. However, even with the introduction of FEM techniques, singular behavior of the integrand might be problematic, which limits improvement of the integration accuracy [27].

Our paper is motivated by many efforts devoted to introducing FEM techniques into GIE [28–30]. In particular, we revisit finite element based Green’s function integral equations (FEGIE) by borrowing the mature techniques of FEM (triangle discretization and quadrature rule) into GIE. In contrast to aforementioned work [15, 24, 31–36], FEGIE can be discretized using either dependent variables or auxiliary variables. We find that auxiliary variable based FEGIE performs well in terms of convergence stability and computational load. Auxiliary variables are field values of the coordinates of the quadrature points inside each element. We can easily interpolate auxiliary variables by using dependent variables. The purpose of introducing the auxiliary variables is to confine the singular behavior of integrand inside its own triangle element such that the singular term can be further approximated analytically. Notably, the introduction of auxiliary variables can improve singular contribution from its own element as well as other elements. Lastly, our FEGIE scheme with auxiliary discretizing GIE is more efficient and shows a stable and fast convergence with less dependent variables. We also note that the singular behavior of Green’s function in 2D scattering is rather weak, thus extending our approach to the 3D scattering problems with strong singular behaviors, though beyond the scope of this paper, is certainly interesting to investigate in the future.

The paper is organized as follows. In Section 2, we give a brief yet complete description of Green’s function integral equation, accompanied with two different implementation strategies, i.e., type-1 FEGIE and type-2 FEGIE [37]. In Section 3, we use FEGIE to calculate the far-field pattern of single circular rod, 3-layered rod, unidirectional 5-layered rod, and a complex scatterer with corrugated surface. We find that the two schemes in FEGIE yield considerably accurate far-field pattern compared with that obtained from Mie theory [38–41] or commercial software package: Multiphysics COMSOL [42], and yield more accurate far-field pattern compared with the standard discretization scheme, i.e., staircase approximated with square grids in the case of weak singularity. Finally, Section 4 concludes the paper.

2. Principles of finite element based Green’s function integral equation

Without losing the generality, we illustrate the seamless integration of finite element techniques into the volume GIE with a better geometric approximation as well as smaller number of degrees of freedom (DOFs). Notably, the extension to other types of Green’s function integral equation, i.e., surface GIE equation, is straightforward and can be implemented in a similar fashion.

2.1. Introduction of Green’s function integral equation in the modelling of light scattering

We give a brief introduction to the volume GIE for completeness. Provided a current source with its current density $J_{s} (r)$ distributed in the domain Ω_s, the electric field $E (r)$ in the whole computational domain Ω can be given by

E (r) = - i ω μ_{0} \int_{Ω_{s}} \bar{g} (r, r^{'}) J_{s} (r^{'}) d V,

where

\bar{g} (r, r^{'})

is the dyadic Green’s function determined by

[\nabla \times \nabla \times - k_{0}^{2} {\bar{ε}}_{r} (r)] \bar{g} (r, r^{'}) = {\bar{I}}_{0} δ (r - r^{'})

.

Next, we consider a plane wave ( $E_{0} (r)$ ) incident upon a scatterer (dielectric function given by ${\bar{ϵ}}_{s} (r)$ ) embedded in the hosting medium ( ${\bar{ϵ}}_{b g}$ ), and we intend to compute the scattering field $E_{s} (r)$ in Ω. Importantly, both the incident field $E_{0}$ and the total field $E (r) = E_{0} (r) + E_{s} (r)$ satisfy the following homogeneous wave equation,

[\nabla \times \nabla \times - k_{0}^{2} {\bar{ε}}_{b g}] E_{0} (r) = 0,

[\nabla \times \nabla \times - k_{0}^{2} {\bar{ε}}_{r} (r)] E (r) = 0,

According to the linearity of vector wave equation, one can derive the relation between total and incident field via the background Green’s function,

E (r) = E_{0} (r) + \int_{r^{'} \in Ω} \bar{g} (r, r^{'}) k_{0}^{2} (ε (r^{'}) - ε_{r e f}) E (r^{'}) d v .

Equation (4) is the volume GIE in the full vector form, describing the scattering of three-dimensional (3D) scatterers. In order to illustrate the basic concepts of how the FEM techniques can be useful for GIE, we consider a 2D scattering object under the illumination of TE polarization, in which case Eq. (4) can be reduced into a scalar form,

E_{z} (r) = E_{0}^{z} (r) + \int_{r^{'} \in Ω} g (r, r^{'}) k_{0}^{2} (ε (r^{'}) - ε_{r e f}) E_{z} (r^{'}) d v,

where

g (r, r^{'})

is the 2D Green’s function, whose explicit expression and its usage to calculate the far-field pattern given in Appendix A.

2.2. Finite element based Green’s function integral equations

In the following, we proceed to discuss the geometric discretization of the scatterer as well as the numerical discretization of GIE. The standard procedure of implementing FEM calculations basically contains 4 steps [29]: (1) using the grids or elements to approximate complex structures; (2) selecting proper polynomial functions with order n to interpolate the field within an element; (3) using Galerkin method to formulate FEM problems defined by certain partial differential equations; (4) Assembling the matrix and solving the set of linear equations. Our proposed FEGIE method has the same 4 steps, except that the FEM formulation in step (3) is replaced by the formulation of the discretized GIE. In FEM, the dependent variables are discretized values of the field, the coordinates of which are carefully chosen according to the grids and the interpolated functions. With the assistance of quadrature rules, the numerical integration in FEM can be highly accurate by interpolating the field values at quadrature points. In contrast, the convergence of numerical integration in the discretized GIE is more difficult to achieve than that in FEM, due to the existence of singular terms in the integrand of Eq. (5). The singular terms are generated at $r$ = $r^{'}$ in the Green’s function. This is in conflict with the standard FEM implementation, since the dependent variables are defined at the vertices shared by a number of triangles, which causes the sharing of singular behavior of the integrand among neighbour triangles.

To overcome this difficulty, the auxiliary variables, i.e., the field values used in the discretization of GIE, is introduced to remove the singular behavior of integrand from the nearby triangles. Once the singular behavior is limited inside its own triangle, analytical approximation or singularity subtraction techniques can be applied to obtain a relatively accurate integration. The auxiliary variables can be interpolated from the dependent variable, but are not necessarily the same as the dependent variables. This subtle difference between auxiliary variables and dependent variables that has not been reported in literature [4, 24, 25, 27] leads to two different strategies of implementing FEGIE, i.e., named as type-1 FEGIE and type-2 FEGIE, which will be discussed in this paper. In type-1 FEGIE, auxiliary variables and dependent variables are indeed identical, which is easier to implement but contains redundant DOFs. In type-2 FEGIE, auxiliary variables and dependent variable are different, which leads to the significant reduction of DOFs and conceptual difficulty that will be addressed shortly.

Fig. 1 A circular disk is discretized into 15 triangles which are further used in (a) type-1 FEGIE and (b) type-2 FEGIE. (c) A small portion, indicated by light gray, of the circular disk of (b) with only 4 triangles and 6 vertices. In type-2 FEGIE, the dependent variables are the vertices indicated by the solid blue dots, while the auxiliary variables associated with each triangle are the quadrature points indicated by red diamonds inside the corresponding triangles.

Download Full Size | PDF

2.2.1. Type-1 FEGIE

We first consider a simple scenario, i.e., type-1 FEGIE. As an example shown in Fig. 1(a), the locations of dependent variables and auxiliary variables in type-1 FEGIE coincide at the origins, i.e., indicated by the blue dots, of the in-circles of 15 triangles. With discretized triangles labelled by i, one can reformulate Eq. (5) into discretized GIE given by

E_{i} = E_{0, i} + \sum_{j} g_{i j} k_{0}^{2} (ε_{j} - ε_{r e f}) E_{j} A_{j}

where A_j is the area of the jth triangle. As

i = j

, the integrand in calculating g_ii becomes singular and can be calculated by approximating the triangle by a circle with the same area. Taking the gray triangle in Fig. 1(a) for example, the equivalent circle is shown as the red circle, which has the radius given by

r = \sqrt{\frac{A_{10}}{π}}

. This treatment leads to a relatively accurate approximation of g_ii by analytically integrating the singular integrand of Eq. (5) within the red circle, if the three sides of the triangle have approximately the same length. Once the descretized GIE, i.e., Eq. (6), is known, it is straight forward to calculate the dependent variable E_i, and thus the field values at all the other coordinates (see explicit expression for g_ii and for the far-field pattern calculation in appendix A).

2.2.2. Type-2 FEGIE and its assembling strategy

In contrast to type-1 FEGIE, we consider the second scenario where the dependent variables and the auxiliary variables are different, as illustrated in Fig. 1(b). As shown in Fig. 1(b), the dependent variable are field values at the vertices, while auxiliary variables that are used to discretize GIE are the field values at the red diamonds. Though the triangulation in Figs. 1(a) and 1(b) is identical for the same circular disk, the different treatments of the dependent variables and auxiliary variables lead to different implementations of FEGIE, especially assembling the matrices of individual triangles into a single matrix with global index, which is one of the main concerns of this subsection. To further illustrate our strategy to implement the discretized GIE, i.e., the aforementioned step (3) similar to FEM modelling, we consider the simplest discretization of a scatterer (containing only 6 vertices and 4 triangles), as shown in Fig. 1(c), and show how to construct type-2 FEGIE through this simple example with first order quadrature rule, accompanied with an efficient assembling strategy.

Firstly, we use six dependent variables, i.e., the field values at the 6 vertices, to interpolate field values at all the other locations, including the auxiliary variables with the locations indicated by the red diamonds. For instance, as shown in Fig. 1(c), the field value $u^{i} (x, y)$ at position $(x, y)$ inside the ith triangle is given by $u^{i} (x, y) = (\begin{matrix} 1 & x & y \end{matrix}) {(\begin{matrix} 1 & x_{1}^{i} & y_{1}^{i} \\ 1 & x_{2}^{i} & y_{2}^{i} \\ 1 & x_{3}^{i} & y_{3}^{i} \end{matrix})}^{- 1} (\begin{matrix} u_{1}^{i} \\ u_{2}^{i} \\ u_{3}^{i} \end{matrix})$ , where $u_{l}^{i}$ is the dependent variable at the lth vertex in the ith triangle (l can also be considered as the local node number), $l \in (1, 2, 3)$ and $i \in (1, 2, 3, 4)$ . Accordingly, the 3 auxiliary variables inside one triangle can be given by $E_{i} (r_{p}) = α_{(i, r_{p})}^{l} u_{l}^{i}$ , where $r_{p}$ represents the coordinate of the pth quadrature point, and $α_{(i, r_{p})}^{l}$ denotes the lth interpolating function corresponding to $u_{l}^{i}$ , and $l \in (1, 2, 3)$ . We follow Einstein summation convention for the indices of l, which are repeated both in the superscript of α and subscript of u.

Secondly, we construct the discretized equation of Eq. (6) for the auxiliary variables, i.e., field value at three quadrature points inside a single triangle as follows,

(\begin{matrix} α_{(i, r_{1})}^{l} u_{l}^{i} \\ α_{(i, r_{2})}^{l} u_{l}^{i} \\ α_{(i, r_{3})}^{l} u_{l}^{i} \end{matrix}) = (\begin{matrix} E_{0}^{(i, r_{1})} \\ E_{0}^{(i, r_{2})} \\ E_{0}^{(i, r_{3})} \end{matrix}) + \sum_{j} k_{0}^{2} Δ ε \frac{A_{j}}{3} (\begin{matrix} g_{(i, r_{1})}^{(j, r_{p})} α_{(j, r_{p})}^{l} u_{l}^{j} \\ g_{(i, r_{2})}^{(j, r_{p})} α_{(j, r_{p})}^{l} u_{l}^{j} \\ g_{(i, r_{3})}^{(j, r_{p})} α_{(j, r_{p})}^{l} u_{l}^{j} \end{matrix}),

where Einstein summation convention is applied for the indices of l and

r_{p}

, and

(i, r_{p})

denotes the quadrature point (

r_{p} \in (r_{1}, r_{2}, r_{3})

) in the ith triangle, the label of interpolating function

l \in (1, 2, 3)

,

i / j

indexes the triangle and does not follow Einstein summation rule. In Eq. (7), the 3 auxiliary variables in triangle i on the left hand side equal to the sum of two terms from the right hand side: the first term is the incident field, the second term is the sum of the contribution from all triangles including the same triangle while

j = i

and the other triangles while

j \neq i

. A_j is the area of the triangle j.

Lastly, our final goal is to assemble the matrices of individual triangles into a single matrix with global index, which further involves two steps. Step 1: Applying Eq. (7) to the simple scatterer shown in Fig. 1(c) for the 4 triangles one by one, one shall obtain the following discretized GIE in the matrix form as follows,

(\begin{matrix} {\bar{α}}_{11} & 0 & 0 & 0 \\ 0 & {\bar{α}}_{22} & 0 & 0 \\ 0 & 0 & {\bar{α}}_{33} & 0 \\ 0 & 0 & 0 & {\bar{α}}_{44} \end{matrix}) (\begin{matrix} U_{1} \\ U_{2} \\ U_{3} \\ U_{4} \end{matrix}) = (\begin{matrix} E_{01} \\ E_{02} \\ E_{03} \\ E_{04} \end{matrix}) + (\begin{matrix} {\bar{g}}_{11} & {\bar{g}}_{12} & {\bar{g}}_{13} & {\bar{g}}_{14} \\ {\bar{g}}_{21} & {\bar{g}}_{22} & {\bar{g}}_{23} & {\bar{g}}_{24} \\ {\bar{g}}_{31} & {\bar{g}}_{32} & {\bar{g}}_{33} & {\bar{g}}_{34} \\ {\bar{g}}_{41} & {\bar{g}}_{42} & {\bar{g}}_{43} & {\bar{g}}_{44} \end{matrix}) (\begin{matrix} U_{1} \\ U_{2} \\ U_{3} \\ U_{4} \end{matrix}),

where

U_{i} = {[u_{1}^{i}, u_{2}^{i}, u_{3}^{i}]}^{T}

is a 3 × 1 column vector, i.e., the 3 dependent variables for triangle i, and

{\bar{α}}_{i i}

and

{\bar{g}}_{i j}

are 3 × 3 matrices as given by Eq. (7), explicitly

{\bar{α}}_{i i} = (\begin{matrix} α_{(i, r_{1})}^{1} & α_{(i, r_{1})}^{2} & α_{(i, r_{1})}^{3} \\ α_{(i, r_{2})}^{1} & α_{(i, r_{2})}^{2} & α_{(i, r_{2})}^{3} \\ α_{(i, r_{3})}^{1} & α_{(i, r_{3})}^{2} & α_{(i, r_{3})}^{3} \end{matrix})

, and

{\bar{g}}_{i j} = k_{0}^{2} Δ ε \frac{A_{j}}{3} \sum_{r_{l} \in (r_{1}, r_{2}, r_{3})} (\begin{matrix} g_{(i, r_{1})}^{(j, r_{l})} & g_{(i, r_{1})}^{(j, r_{l})} & g_{(i, r_{1})}^{(j, r_{l})} \\ g_{(i, r_{2})}^{(j, r_{l})} & g_{(i, r_{2})}^{(j, r_{l})} & g_{(i, r_{2})}^{(j, r_{l})} \\ g_{(i, r_{3})}^{(j, r_{l})} & g_{(i, r_{3})}^{(j, r_{l})} & g_{(i, r_{3})}^{(j, r_{l})} \end{matrix}) (\begin{matrix} α_{(j, r_{l})}^{1} & α_{(j, r_{l})}^{2} & α_{(j, r_{l})}^{3} \\ α_{(j, r_{l})}^{1} & α_{(j, r_{l})}^{2} & α_{(j, r_{l})}^{3} \\ α_{(j, r_{l})}^{1} & α_{(j, r_{l})}^{2} & α_{(j, r_{l})}^{3} \end{matrix})

.

Table 1. Assembling Transformation

View Table

In step 1, the matrix form of GIE in Eq. (8) is obtained simply by stacking the triangles one by one, thus the shared vertices among the neighbour triangles are not taken into account. Therefore, $U = {[U_{1}, U_{2}, U_{3}, U_{4}]}^{T}$ is a 12 × 1 column vector, in contrast to 6 DOFs in total shown in Fig. 1(c). In step 2, wedirectly rearrange the simply stacked local matrices, i.e., Eq. (8) into a global indexed matrix via the assembling transformation. To that end, we introduce a global indexing of dependent variables $U^{g}$ as well as the assembling transformation matrix ${\bar{T}}_{a s m}$ , i.e, $U^{g} = {\bar{T}}_{a s m} U$ . Thus, Eq. (8) denoted as $[{\bar{α}}_{i i}] U = E^{i n c} + [{\bar{g}}_{i j}] U$ , can be reformulated into the final matrix form that can be readily solved using the standard linear algebra,

({\bar{T}}_{a s m} [{\bar{α}}_{i i}] {\bar{T}}_{a s m}^{T}) U^{g} = {\bar{T}}_{a s m} E^{i n c} + ({\bar{T}}_{a s m} [{\bar{g}}_{i j}] {\bar{T}}_{a s m}^{T}) U^{g},

where

[{\bar{α}}_{i i}]

(

[{\bar{g}}_{i j}]

) is the full matrix of the block matrix

{\bar{α}}_{i i}

(

{\bar{g}}_{i j}

). The assembling transformation matrix

{\bar{T}}_{a s m}

can be obtained from the triangulation matrix, as tabulated in Table 1. The triangulation matrix indicates how the vertices (the global node number) are connected to the local node number in each triangle. The column number of

{\bar{T}}_{a s m}

equals to the total number of elements of the triangulation matrix

E_{t r i}

, while the row number equals to the total number of vertices. By flattening

E_{t r i}

row by row, one can construct

{\bar{T}}_{a s m}

by mapping the row to the element of the flattened

E_{t r i}

, with the location of nonzero element in the same row corresponding to the value of that element in the flattened

E_{t r i}

.

Fig. 2 (a) Triangle discretization of the dielectric rod. (b) Far-field pattern calculated using type-1 FEGIE, type-2 FEGIE, staircase GIE, and Mie theory. The radius of the scatterer is 0.6 μm, and the vacuum wavelength is 0.633 μm. The numbersof DOFs in type-1 FEGIE, type-2 FEGIE, and staircase are 1080, 595, and 1085, respectively.

Download Full Size | PDF

3. Results and discussions

We examine the far-field pattern of 2D light scattering using our proposed method, i.e., type-1 FEGIE and type-2 FEGIE, in comparison with staircase GIE [8, 24] and Mie theory [41]. As for staircaseGIE, we assign the average value of the dielectric constant for the squares located on the boundary to improve the accuracy. In the first example, we study light scattering by a dielectric rod (infinitely long along z-axis, ϵ_r=4), with the background being air. The incident light is TE polarized, propagating along x-axis with vacuum wavelength 0.633 μm. Without explicit mention, the same incident conditions apply to the scattering problems throughout this section. The dielectric rod is discretized with triangle elements shown in Fig. 2(a), and the farfield pattern calculated using Mie theory (red dotted line), type-1 FEGIE (blue dashed line), type-2 FEGIE (black dotted line), staircase GIE (magenta dash-dot line). As for the scatterer with simple geometric shape, our type-1 FEGIE and type-2 FEGIE yield reasonably accurate far-field calculations, as benchmarked against the results obtained from Mie theory as well as staircase GIE. Staircase GIE captures the main lobes slightly better thanType-1 FEGIE and TYPE-2 GIE, but a close examination also shows that type-2 FEGIE captures the sharp features of the far-field pattern much better than that of type-1 FEGIE and staircase GIE, which can be attributed to the more integration points used in the numerical integral for each triangles. As a side remark, the number of DOFs in type-1 FEGIE (staircase GIE) is the same as the number of the triangles (squares), i.e. approximately 1085. While the number of DOFs used in type-2 FEGIE is equal to the number of vertices, i.e., 595, which is less than those in type-1 FEGIE and staircase GIE. The smaller number of DOFs can be considered as the second advantage of type-2 FEGIE, apart from the ability of capturing the sharp features in the farfield pattern.

Fig. 3 Far-field pattern of the light scattering by (a) the 3-layered rod and (b) the 5-layered rod calculated using type-1 FEGIE, type-2 FEGIE, and Mie theory. For the 3-layered (5-layered) rod, the numbers of DOFs in type-1 FEGIE and type-2 FEGIE are 3186 (4812) and 1674 (2507).

Download Full Size | PDF

In the second example as shown in Fig. 3, we proceed to discuss light scattering of 3-layered rod and 5-layered rod using type-1 FEGIE (blue dashed line) and type-2 FEGIE (black dotted line), which are again compared with finite element calculations from Mie theory (red dotted line). The staircase GIE calculation of this problem is poorly converged, which is not shown here. As shown in the inset of (a), the radii of the three layers rod are specified to be r₁=0.3 μm, r₂=0.7 μm, r₃=1.0 μm, and the dielectric constants of the three layers are given as $ϵ_{1} = 3$ , $ϵ_{2} = 4$ , $ϵ_{3} = 5$ . The radii of the 5-layered rod [43], as sketched in the inset of (b), are given by r₁=0.0158 μm, r₂=0.0285 μm, r₃=0.0475 μm, r₂=0.0601 μm, r₅=0.633 μm. The corresponding dielectric constants are given as $ε_{1} = 0.9814$ , $ε_{2} = 4.7676$ , $ε_{3} = 1.7075$ , $ϵ_{4} = 1.4711$ , $ε_{5} = 1.2025$ . We examine the differences of the far-field patterns for the 3/5-layered rod calculated from type-1 FEGIE and type-2 FEGIE, as shown in Fig. 3(a/b). In Fig. 3(a), both type-1 FEGIE and type-2 FEGIE can capture the overall features of far-field pattern, including the angular positions and amplitudes of the main lobes and side lobes (SLs), in consistency with calculations from Mie theory. Notably, the far-field pattern from type-2 FEGIE yields excellent agreement with that of FEM calculations, except slight deviations in the main lobe and a few side lobes (SL1, SL2, SL10, SL11). In contrast, the far-field pattern from type-1 FEGIE deviates in most of the side lobes, and close examination also reveals that the far-field pattern (blue dashed line)is not symmetric with respect with the x-axis. As a side remark, the number of DOFs used in type-1 FEGIE and type-2 FEGIE are 3186 and 1674 respectively. The light scattering of a 5-layered rod is examined in Fig. 3(b). The particular choices of the geometric parameters and dielectric functions for the 5 layers rod yields highly directional radiation pattern [43], which can be well captured using our proposed type-1 FEGIE and type-2 FEGIE. Importantly, our proposed approaches yield accurate results and using much less DOFs, which is highly promising for modelling light scattering of optically large structures.

Fig. 4 (a) Far-field patterns of the light scattering by a complex scatterer with corrugated surface calculated using type-1 FEGIE, type-2 FEGIE, staircase GIE, and Mie theory. (b/c) Triangle/Staircase approximation of the surface-corrugated scatterer. The numbers of DOFs in type-1 FEGIE, type-2 FEGIE, staircase, and COMSOL are 2586, 1407, 2590, and 280525 respectively.

Download Full Size | PDF

As the third example, we study the far-field pattern of a geometrically complicated scatterer with corrugated surface as shown in Fig. 4. The far-field pattern is calculated using type-1 FEGIE (blue dashed line), type-2 FEGIE (black dotted line), staircase FIE (magenta dash-dot line)) and COMSOL (red dotted line). Similar to the staircase GIE used in the first example shown in Fig. 2, the averaged dielectric constant for the squares on the boundary is used in the staircase GIE. The scatterer with corrugated surface discretized with triangle and square elements are shown in Figs. 4(b) and 4(c), which have comparable number of elements, i.e., 2586 triangles versus 2590 squares. Figure 4(a) shows a few side lobes of the far-field, as indicated by the dashed rectangle shown in the inset of the full far-field pattern. As can be seen from the inset and main panel of Fig. 4(a), the main lobe and two side lobes (SL3 and SL5) calculated from type-2 FEGIE shows excellent agreements withthat calculated from finite element method using COMSOL. In contrast, the type-1 FEGIE and staircase GIE yield considerable deviation of the main lobe and all the other side lobes. As regarding to other side lobes shown in Fig. 4(a), type-2 FEGIE yet performs better in calculating far-field pattern than type-1 FEGIE and staircase GIE. It is worthy to mention that the DOFs used in type-1 FEGIE, type-2 FEGIE, staircase GIE and COMSOL are 2586, 1407, 2590 and 280525.

Fig. 5 Farfield convergence versus the number of DOFs.

Download Full Size | PDF

Lastly, Fig. 5 shows the convergence of the farfield plotted in Fig. 4(a) at θ=0, for type-1 FEGIE (green-square), type-2 FEGIE (blue-circle), and staircase GIE (violet-triangle). The red dash line indicates the reference value of farfield obtained from COMSOL with 1.8 million DOFs. Evidently, our Type-2 FEGIE has a stable and fast convergence as the number of DOFs increases, while Type-1 FEGIE shows oscillations and a slower convergence than that of Type-2 FEGIE due to the redundant DOFs. Staircase GIE has the slowest convergence due the poor approximation of the geometry as well as the redundant DOFs. Thus, as evident from the far-field calculations and the listed DOFs used in the four approaches, our proposed type-2 FEGIE indeed shows the advantages of higher accuracy and computational efficiency due to the significant reduction of DOFs.

4. Conclusion

In conclusion, we introduce finite element techniques into the Green’s function integral equation, with a better and more flexible approximation of the geometry of the scatterers and a more accurate numerical integral. The ameliorative performance is due to the fact that more quadrature points chosen inside each element yield the numerical integrals more accurately. On the other hand, the introduction of auxiliary variables leads to an easier treatment of singular terms which are confined inside each element and less number of DOFs. We have calculated the farfield patterns of several concrete scatterers with two types of FEGIE. We find that the two discretization schemes in FEGIE yield considerably accurate far-field pattern compared with that obtained using COMSOL, and yield more accurate far-field pattern compared with the standard discretization scheme, i.e., stair-cased approximated with square grids, which is highly promising for modelling light scattering of optically large and complex structures. Besides, the number of DOFs used in type-2 FEGIE is less than those in type-1 FEGIE and staircase GIE.

Appendix A Scalar Green’s function, analytical approximation of self-term integration, and the far-field calculations in 2D Green’s function integral equation

The scalar Green’s function of the TE polarization in 2D light scattering problems is given by $g (r, r^{'}) = \frac{1}{4 i} H_{0}^{(2)} (k_{0} n_{r e f} | r - r^{'} |)$ . As $i = j$ in Eq. (6), g_ii represents the interaction of a discrete element with itself, it can be approximated with

g_{i i} = \frac{1}{A_{i}} \int_{A_{i}} g (r, r^{'}) d A^{'}

If we adopt square discrete elements, a square element can be approximated by circular one with the same area and same center position, Eq. (4) can be evaluated analytically resulting in

g_{i i} \approx \frac{1}{2 i {(k_{0} n_{r e f} a)}^{2}} (k_{0} n_{r e f} a H_{1}^{(2)} (k_{0} n_{r e f} a) - \frac{2 i}{π})

For far field, since $| r^{'} | / | r | ≪ 1$ , the Green’s function can be approximated with

g^{f f} (r, r^{'}) = \frac{1}{4} \sqrt{\frac{2}{π k r}} e^{\frac{- i π}{4}} e^{- i k r} e^{i k \frac{r}{r} \cdot r^{'}}

Then the far field can be obtained as

E_{S C}^{f f} (r) = \frac{1}{4} \sqrt{\frac{2}{π k r}} e^{\frac{- i π}{4}} e^{- i k r} \int k_{0}^{2} (ε (r^{'}) - ε_{r e f}) E (r^{'}) e^{i k \frac{r}{r} \cdot r^{'}} d A^{'}

Funding

National Natural Science Foundation of China (NSFC) (11874026, 61735006, 61775063); Fundamental Research Funds for the Central Universities (HUST: 2017KFYXJJ027); National Key Research and Development Program of China (2017YFA0305200); Research Grants Council of Hong Kong SAR (CityU 21302018).

Acknowledgment

We thank Thomas M. Søndergaard for fruitful discussions and a careful reading of the manuscript.

References

1. E. N. Economou, Green’s Functions in Quantum Physics(Springer, 1983). [CrossRef]

2. A. V. Lavrinenko, J. Lægsgaard, N. Gregersen, F. Schmidt, and T. Søndergaard, Numerical Methods in Photonics (CRC, 2014). [CrossRef]

3. L. Novotny and B. Hecht, Principles of Nano-Optics(Cambridge University, 2012). [CrossRef]

4. W. C. Chew, M. S. Tong, and B. Hu, Integral Equation Methods for Electromagnetic and Elastic Waves (Morgan & Claypool Publishers, 2009).

5. O. J. Martin, A. Dereux, and C. Girard, “Iterative scheme for computing exactly the total field propagating in dielectric structures of arbitrary shape,” J. Opt. Soc. Am. A 11(3), 1073–1080 (1994). [CrossRef]

6. J. Jung, T. Søndergarrd, T. G. Pedersen, K. Pedersen, A. N. Larsen, and B. B. Nielsen, “Dyadic Green’s functions of thin films: Applications within plasmonic solar cells,” Phys. Rev. B 83(8), 085419 (2011). [CrossRef]

7. Y. T. Chen, Y. Zhang, and A. F. Koenderink, “General point dipole theory for periodic metasurfaces: magnetoelectric scattering lattices coupled to planar photonic structures,” Opt. Express 25(18), 21358–21378 (2017). [CrossRef] [PubMed]

8. T. Søndergaard, “Modeling of plasmonic nanostructures: Green’s function integral equation methods,” Phys. Status Solidi (b) 244(10), 3448–3462 (2007). [CrossRef]

9. S. N. Chandler-Wilde and D. C. Hothersall, “Efficient calculation of the Green function for acoustic propagation above a homogeneous impedance plane,” J. Sound Vibr. 180(5), 705–724 (1995). [CrossRef]

10. M. A. Yurkin and M. I. Mishchenko, “Volume integral equation for electromagnetic scattering: Rigorous derivation and analysis for a set of multilayered particles with piecewise-smooth boundaries in a passive host medium,” Phys. Rev. A 97(4), 043824 (2018). [CrossRef]

11. J. P. Kottmann and O. J. F. Martin, “Accurate solution of the volume integral equation for high-permittivity scatterers,” IEEE Trans. Antennas Propag. 48(11), 1719–1726 (2000). [CrossRef]

12. W. C. Gibson, The Method of Moments in Electromagnetics (CRC, 2015).

13. R. F. Harrington, Field Computation by Moment Methods(Macmillan, 1968).

14. N. Morita, N. Kumagai, and J. R. Mautz, Integral Equation Methods for Electromagnetics (Aetech House, 1990).

15. J. Jung and T. Søndergaard, “Green’s function surface integral equation method for theoretical analysis of scatterers close to a metal interface,” Phys. Rev. B 77(24), 245310 (2008). [CrossRef]

16. U. Hohenester and A. Trügler, “MNPBEM-A Matlab toolbox for the simulation of plasmonic nanoparticles,” Comput. Phys. Commun. 183(2), 370–381 (2012). [CrossRef]

17. U. Hohenester, “Simulating electron energy loss spectroscopy with the MNPBEM toolbox,” Comput. Phys. Commun. 185(3), 1177–1187 (2014). [CrossRef]

18. J. Waxenegger, A. Trügler, and U. Hohenester, “Plasmonics simulations with the MNPBEM toolbox: Consideration of substrates and layer structures,” Comput. Phys. Commun. 193, 138–150 (2015). [CrossRef]

19. S. B. Wang, H. H. Zheng, J. J. Xiao, Z. F. Lin, and C. T. Chan, “Fast multipole boundary element method for three dimensional electromagnetic scattering problem,” Int. J. Comput. Mater. Sci. Eng. 1(04), 1250038 (2012).

20. D. W. Prather, M. S. Mirotznik, and J. N. Mait, “Boundary integral methods applied to the analysis of diffractive optical elements,” J. Opt. Soc. Am. A 14(1), 34–43 (1997). [CrossRef]

21. V. Myroshnychenko, E. Carbó-Argibay, I. Pastoriza-Santos, J. Pérez-Juste, L. M. Liz-Marzán, and F. J. García de Abajo“Modeling the optical response of highly faceted metal nanoparticles with a fully 3D boundary element method,” Adv. Mater. 20(22), 4288–4293 (2008). [CrossRef]

22. F. J. García de Abajo and A. Howie“Retarded field calculation of electron energy loss in inhomogeneous dielectrics,” Phys, Rev. B 65(11), 115418 (2002). [CrossRef]

23. R. Kress, “Boundary integral equations in time-harmonic acoustic scattering,” Mathl. Comput. Modelling 15(3–5), 229–243 (1991). [CrossRef]

24. T. Søndergaard, Green’s Function Integral Equation Methods in Nano-Optics (CRC, 2018).

25. A. M. Kern and O. J. F. Martin, “Surface integral formulation for 3D simulations of plasmonic and high permittivity nanostructures,” J. Opt. Soc. Am. A 26(4), 732–740 (2009). [CrossRef]

26. R. Rodríguez-Oliveros and J. A. Sánchez-Gil, “Localized surface-plasmon resonances on single and coupled nanoparticles through surface integral equations for flexible surfaces,” Opt. Express 19(13), 12208–12219 (2011). [CrossRef] [PubMed]

27. T. V. Raziman, W. R. C. Somerville, O. J. F. Martin, and E. C. Le Ru, “Accuracy of surface integral equation matrix elements in plasmonic calculations,” J. Opt. Soc. Am. B 32(3), 485–492 (2015). [CrossRef]

28. D. Chen, W. Cai, B. Zinser, and M. H. Cho, “Accurate and efficient Nyström volume integral equation method for the Maxwell equations for multiple 3-D scatterers,” J. Comput. Phys. 321, 303–320 (2016). [CrossRef]

29. J. Jin, The Finite Element Method in Electromagnetics(Wiley, 2002).

30. Y. T. Chen, T. R. Nielsen, N. Gregersen, P. Lodahl, and J. Mørk, “Finite-element modeling of spontaneous emission of a quantum emitter at nanoscale proximity to plasmonic waveguides,” Phys. Rev. B 81(12), 125431 (2010). [CrossRef]

31. T. Søndergaard and S. Bozhevolnyi, “Strip and gap plasmon polariton optical resonators,” Phys. Status Solidi (b) 245(1), 9–19 (2008). [CrossRef]

32. T. Søndergaard and S. Bozhevolnyi, “Slow-plasmon resonant nanostructures: Scattering and field enhancements,” Phys. Rev. B 75(7), 073402 (2007). [CrossRef]

33. J. Jung, T. Søndergaard, J. Beermann, A. Boltasseva, and S. I. Bozhevolnyi, “Theoretical analysis and experimental demonstration of resonant light scattering from metal nanostrips on quartz,” J. Opt. Soc. Am. B 26(1), 121–124 (2009). [CrossRef]

34. V. Siahpoush and T. Søndegaard, ard J. Jung, “Green’s function approach to investigate the excitation of surface plasmon polaritons in a nanometer-thin metal film,” Phys. Rev. B 85(7), 075305 (2012). [CrossRef]

35. T. Søndergaard, V. Siahpoush, and J. Jung, “Coupling light into and out from the surface plasmon polaritons of a nanometer-thin metal film with a metal nanostrip,” Phys. Rev. B 86(8), 085455 (2012). [CrossRef]

36. T. Søndergaard and S. I. Bozhevolnyi, “Theoretical analysis of plasmonic black gold: periodic arrays of ultra-sharp grooves,” New J. Phys. 15(1), 013034 (2013). [CrossRef]

37. The matlab scripts of the practical impelemntation can be requested via yuntian@hust.edu.cn

38. J. Schäfer, S. C. Lee, and A. Kienle, “Calculation of the near fields for the scattering of electromagnetic waves by multiple infinite cylinders at perpendicular incidence,” J. Quant. Spectrosc. Radiat. Transf. 113(16), 2113–2123 (2012). [CrossRef]

39. S. C. Lee, “Dependent scattering of an obliquely incident plane wave by a collection of parallel cylinders,” J. Appl. Phys. 68(10), 4952–4957 (1990). [CrossRef]

40. E. M. Loebl, Scattering of Light and Other Electromagnetic Radiation (Academic, 1969).

41. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles(Wiley, 2008).

42. Comsol, “Comsol Multiphysics Modeling Software,”https://www.comsol.com/.

43. S. Arslanagić and R. W. Ziolkowski, “Highly subwavelength, superdirective cylindrical nanoantenna,” Phys. Rev. Lett. 120(23), 237401 (2018). [CrossRef]

Finite element based Green’s function integral equation for modelling light scattering

Abstract

1. Introduction

2. Principles of finite element based Green’s function integral equation

2.1. Introduction of Green’s function integral equation in the modelling of light scattering

2.2. Finite element based Green’s function integral equations

2.2.1. Type-1 FEGIE

2.2.2. Type-2 FEGIE and its assembling strategy

3. Results and discussions

4. Conclusion

Appendix A Scalar Green’s function, analytical approximation of self-term integration, and the far-field calculations in 2D Green’s function integral equation

Funding

Acknowledgment

References

Cited By

Figures (5)

Tables (1)

Equations (13)

Optics Express