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

TL;DR

This paper introduces a finite element based Green's function integral equation method for light scattering, improving geometric approximation and numerical accuracy over traditional methods, demonstrated through complex scatter examples.

Contribution

The paper integrates finite element techniques into Green's function integral equations, enhancing geometric modeling and numerical integration accuracy for light scattering simulations.

Findings

Improved geometric approximation of scatterers.

Enhanced numerical integration accuracy.

Effective modeling of complex, large-scale scatterers.

Abstract

We revisit the volume Green's function integral equation for modelling light scattering with discretization strategies as well as numerical integration recipes borrowed from finite element method. The merits of introducing finite element techniques into Green's function integral equation are apparent. Firstly, the finite element discretization provides a much better geometric approximation of the scatters, compared with that of the conventional discretization method using staircase approximation. Secondly, the accuracy of numerical integral inside one element associated with all different types of Green's function integral equations can be greatly improved by introducing quadrature rules. Within the standard framework of Green's function integral equation, we seamlessly introduce finite element techniques into the Green's function integral equation by introducing the auxiliary variables that confines the singular integrand inside each element, leading to a better and more flexible approximation of the geometry of the scatters and a more accurate numerical integral. We then illustrate the advantages of our finite element based Green's function integral equation method via a few concrete examples in modelling light scattering by optically large and complex scatters.

Tables

Table 1. Table 1: Assembling transformation

Triangulation matrix	Assembling transformation matrix

Equations

$$\mathbf{E}\left(\mathbf{r}\right)=-i\omega {\mu }_0{\int }_{{\Omega }_s}\bar{\mathbf{g}}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right){\mathbf{J}}_{\mathbf{s}}\left({\mathbf{r}}^{\prime }\right)dV,$$

$$[\nabla \times \nabla \times -k_0^2{\bar{\mathbf{\epsilon }}}_{bg}]{\mathbf{E}}_0(\mathbf{r})=0,$$

$$[\nabla \times \nabla \times -k_0^2{\bar{\mathbf{\epsilon }}}_r(\mathbf{r})]{\mathbf{E}}_{tot}(\mathbf{r})=0.$$

$${\mathbf{E}}_{tot}\left(\mathbf{r}\right)={\mathbf{E}}_0\left(\mathbf{r}\right)+{\int }_{{\mathbf{r}}^{\prime }\in \Omega }\bar{\mathbf{g}}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)k_0^2\left(\epsilon \left({\mathbf{r}}^{\prime }\right)-{\epsilon }_{ref}\right)\mathbf{E}\left({\mathbf{r}}^{\prime }\right)dv.$$

$$E_z\left(\mathbf{r}\right)=E_0^z\left(\mathbf{r}\right)+{\int }_{{\mathbf{r}}^{\prime }\in \Omega }g\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)k_0^2\left(\epsilon \left({\mathbf{r}}^{\prime }\right)-{\epsilon }_{ref}\right)E_z\left({\mathbf{r}}^{\prime }\right)dv,$$

$$E_i=E_{0,i}+\sum _j{g}_{ij}{k}_0^2\left({\epsilon }_j-{\epsilon }_{ref}\right)E_j{A}_j$$

$$\left({\begin{array}[]{*{20}{c}}{\alpha_{(i,\bm{r}_{1})}^{l}u_{l}^{i}}\\ {\alpha_{(i,\bm{r}_{2})}^{l}u_{l}^{i}}\\ {\alpha_{(i,\bm{r}_{3})}^{l}u_{l}^{i}}\end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{E_{0}^{(i,\bm{r}_{1})}}\\ {E_{0}^{(i,\bm{r}_{2})}}\\ {E_{0}^{(i,\bm{r}_{3})}}\end{array}}\right)+k_{0}^{2}\Delta\varepsilon\frac{{{A_{i}}}}{3}\left({\begin{array}[]{*{20}{c}}{g_{(i,\bm{r}_{1})}^{(i,\bm{r}_{p})}\alpha_{(i,\bm{r}_{p})}^{l}u_{l}^{i}}\\ {g_{(i,\bm{r}_{2})}^{(i,\bm{r}_{p})}\alpha_{(i,\bm{r}_{p})}^{l}u_{l}^{i}}\\ {g_{(i,\bm{r}_{2})}^{(i,\bm{r}_{p})}\alpha_{(i,\bm{r}_{p})}^{l}u_{l}^{i}}\end{array}}\right)+\sum\limits_{j\neq i}k_{0}^{2}\Delta\varepsilon\frac{{{A_{j}}}}{3}\left({\begin{array}[]{*{20}{c}}{g_{(i,\bm{r}_{1})}^{(j,\bm{r}_{p})}\alpha_{(j,\bm{r}_{p})}^{l}u_{l}^{j}}\\ {g_{(i,\bm{r}_{2})}^{(j,\bm{r}_{p})}\alpha_{(j,\bm{r}_{p})}^{l}u_{l}^{j}}\\ {g_{(i,\bm{r}_{3})}^{(j,\bm{r}_{p})}\alpha_{(j,\bm{r}_{p})}^{l}u_{l}^{j}}\end{array}}\right),$$

$$\left({\begin{array}[]{*{20}{c}}\bar{\bm{\alpha}}_{11}&\bm{0}&\bm{0}&\bm{0}\\ \bm{0}&\bar{\bm{\alpha}}_{22}&\bm{0}&\bm{0}\\ \bm{0}&\bm{0}&\bar{\bm{\alpha}}_{33}&\bm{0}\\ \bm{0}&\bm{0}&\bm{0}&\bar{\bm{\alpha}}_{44}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}\bm{U}_{1}\\ \bm{U}_{2}\\ \bm{U}_{3}\\ \bm{U}_{4}\end{array}}\right){\rm{=}}\left({\begin{array}[]{*{20}{c}}\bm{E}_{01}\\ \bm{E}_{02}\\ \bm{E}_{03}\\ \bm{E}_{04}\end{array}}\right){\rm{+}}\left({\begin{array}[]{*{20}{c}}\bm{\bar{g}}_{11}&\bm{\bar{g}}_{12}&\bm{\bar{g}}_{13}&\bm{\bar{g}}_{14}\\ \bm{\bar{g}}_{21}&\bm{\bar{g}}_{22}&\bm{\bar{g}}_{23}&\bm{\bar{g}}_{24}\\ \bm{\bar{g}}_{31}&\bm{\bar{g}}_{32}&\bm{\bar{g}}_{33}&\bm{\bar{g}}_{34}\\ \bm{\bar{g}}_{41}&\bm{\bar{g}}_{42}&\bm{\bar{g}}_{43}&\bm{\bar{g}}_{44}\end{array}}\right)\left({\begin{array}[]{*{20}{c}}\bm{U}_{1}\\ \bm{U}_{2}\\ \bm{U}_{3}\\ \bm{U}_{4}\end{array}}\right),$$

$$\left({\bar{\mathbf{T}}}_{asm}\left[{\bar{\mathbf{\alpha }}}_{ii}\right]{\bar{\mathbf{T}}}_{asm}^T\right){\mathbf{U}}^g={\bar{\mathbf{T}}}_{asm}{\mathbf{E}}^{inc}+\left({\bar{\mathbf{T}}}_{asm}\left[{\bar{\mathbf{g}}}_{ij}\right]{\bar{\mathbf{T}}}_{asm}^T\right){\mathbf{U}}^g,$$

$$g_{ii}=\frac{1}{A_i}{\int }_{A_i}g\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)d{A}^{\prime }$$

$$g_{ii}\approx \frac{1}{2i{\left(k_0{n}_{ref}a\right)}^2}\left(k_0{n}_{ref}a{H}_1^{\left(2\right)}\left(k_0{n}_{ref}a\right)-\frac{2i}{\pi }\right)$$

$$g^{ff}\left(\mathbf{r},{\mathbf{r}}^{\prime }\right)=\frac{1}{4}\sqrt{\frac{2}{\pi kr}}{e}^{\frac{-i\pi }{4}}{e}^{-ikr}{e}^{ik\frac{\mathbf{r}}{r}\cdot {\mathbf{r}}^{\prime }}$$

$$E_{SC}^{ff}\left(\mathbf{r}\right)=\frac{1}{4}\sqrt{\frac{2}{\pi kr}}{e}^{\frac{-i\pi }{4}}{e}^{-ikr}\int k_0^2\left(\epsilon \left(\mathbf{r}\right)-{\epsilon }_{ref}\right)E\left({\mathbf{r}}^{\prime }\right)e^{ik\frac{\mathbf{r}}{r}\cdot {\mathbf{r}}^{\prime }}d{A}^{\prime }$$

Full text

††journal: oe

1 Introduction

Green’s function integral equation (GIE) is an elegant and efficient approach of calculating the scattering of waves [1, 2], and has been widely applied in modelling light scattering in nano-structures [3, 9, 4, 5, 6, 7, 8] and acoustic scattering [10]. An appealing advantage of GIE is that only the scatter is needed to be discretized, in contrast to other approaches, e.g., finite element method (FEM) and finite difference time-domain (FDTD) method, which require both the scatter and the background to be discretized. This leads to the reduction of the computational overhead, which can be further reduced by transforming the volume GIE [11, 12] into a surface GIE [13, 14, 15, 20, 16, 19, 17, 18, 21, 22, 23, 24], where only the boundary of the scatter needs to be discretized. However, the challenge associated with volume GIE and surface GIE is the singular behavior of integrand, which has attracted considerable attention in the last decade. Several schemes have been developed to address the challenge, include: singularity subtraction [25], improved numerical integration [26, 27, 28, 29], and analytical approximation of the singular terms [30, 31, 32].

Interestingly, considerable efforts have been dedicated to introducing FEM techniques [33, 34, 26], e.g., advanced discretization strategies [32, 35] as well as efficient algorithms of carrying out numerical integrals [5, 36], into GIE. Indeed, the finite element discretization could not only yield an improved geometric approximation, but also improve the numerical accuracy of the integral by introducing the quadrature-rule based integration techniques that are widely used in FEM. For instance, triangle discretization has been introduced to approximate the two-dimensional (2D) scatters in volume GIE [32], or surfaces of 3D scatters in surface GIE [35, 37, 5]. Apart from the various advantages, the accuracy of GIE based on finite element discretization and quadrature integration rule is still compromised by the singular integrand, as studied in Ref [36], where the link between the numerical quadrature order and the integral accuracy is investigated in details. In particular, the authors in Ref [36] find that the singular behavior of Green’s tensor at the source causes a large variation in field amplitude over triangle elements close to it. Such large variation nearby the sources is due to the use of the divergence-conforming elements equipped with Rao-Wilton-Glison basis function [35, 36, 5], which are defined on a pair of triangles sharing the same edge.

In this paper, we elaborate the combination of FEM techniques, i.e., triangle discretization and quadrature rule, into GIE by introducing auxiliary variables. In our proposed finite element based Green’s function integral equation (FEGIE), we follow the same course of implementing GIE that is extensively studied by Søndergaard [39, 32, 40, 41, 42, 16, 43, 38]. In contrast to Søndergaard’s work, we use the auxiliary variables instead of the dependent variables to discretize GIE. The auxiliary variables are the field values with carefully selected coordinates, which coincide with the coordinates of quadrature points inside each element. Given the dependent variables of the triangles and the associated basis functions, the auxiliary variables can be easily interpolated. The purpose of introducing the auxiliary variables is to confine the singular behavior of integrand inside its own triangle, such that the singular term can to a good approximation be integrated analytically one by one. The second benefit of introducing auxiliary variables is the improved integral contributions from other triangles, apart from the singular contribution from its own triangle. Lastly, the auxiliary variables can be seamlessly and conveniently integrated into the current framework of GIE with an efficient assembly strategy. As an illustration of this idea, we use triangle elements to discretize the scatter, and provide two types of implementation strategy to realize FEGIE [44].

The paper is organized as follows. In Section 2, we give a brief yet complete description of volume Green’s function integral equation, accompanied with two different implementation strategies, i.e., type-1 FEGIE and type-2 FEGIE. 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 scatter with corrugated surface. We find that the two schemes in FEGIE yield considerably accurate far-field pattern compared with that obtained from Mie theory [45, 46, 47, 48, 49] or commercial software package: Multiphysics COMSOL [50], and yield more accurate far-field pattern compared with the standard discretization scheme, i.e., staircase approximated with square grids. 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 volume 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 $\bm{J}_{s}(\bm{r})$ distributed in the domain $\Omega_{s}$, the electric field $\bm{E}(\bm{r})$ in the whole computational domain $\Omega$ can be given by

[TABLE]

where $\bar{\bm{g}}\left({{\bm{r}},{\bm{r^{\prime}}}}\right)$ is the dyadic Green’s function determined by $[\nabla\times\nabla\times-k_{0}^{2}\bar{\bm{\varepsilon}}_{r}(\bm{r})]\bar{\bm{g}}(\bm{r},\bm{r}^{\prime})=\bar{\bm{I}}_{0}\delta(\bm{r}-\bm{r}^{\prime})$.

Next, we consider a plane wave ($\bm{E}_{0}(\bm{r})$) incident upon a scatter (dielectric function given by $\bar{\bm{\epsilon}}_{s}(\bm{r})$) embedded in the hosting medium ($\bar{\bm{\epsilon}}_{bg}$), and we intend to compute the scattering field $\bm{E}_{s}(\bm{r})$ in $\Omega$. Importantly, both the incident field $\bm{E}_{0}$ and the total field $\bm{E}_{tot}(\bm{r})=\bm{E}_{0}(\bm{r})+\bm{E}_{s}(\bm{r})$ satisfy the following homogeneous wave equation,

[TABLE]

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

[TABLE]

Equation (3) is the volume GIE in the full vector form, describing the scattering of three-dimensional (3D) scatters. 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 Eq. (3) can be reduced into a scalar form,

[TABLE]

where ${g\left({{\bm{r}},{\bm{r^{\prime}}}}\right)}$ 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 scatter as well as the numerical discretization of GIE. The standard procedure of implementing FEM calculations basically contains 4 steps [33]: (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. (4). The singular terms are generated at $\bm{r}$ = $\bm{r}^{\prime}$ 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 not necessarily the same as the dependent variables. This subtle difference between auxiliary variables and dependent variables that are overlooked in literature [32, 35, 36, 5] 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.

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. (4) into discretized GIEgiven by

[TABLE]

where $A_{j}$ is the area of the $j$th triangle. As $i=j$, the integrand in Eq. (4) 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}}{\pi}}$. This treatment leads to a relatively accurate approximation of $g_{ii}$ by analytically integrating the singular integrand of Eq. (4) within the red circle. Once the descretized GIE, i.e., Eq. (5), 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 $g_{ii}$ and the far-field pattern calculation in appendix A).

2.2.2 Type-2 FEGIE

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 Fig. 1 (a) and Fig. 1 (b) is identical for the same circular disk, the different treatments of the dependent variables and auxiliary variables lead to completely different implementations of FEGIE, which are the main concerns of this paper. 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 scatter (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 and first order interpolating function.

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 $i$th triangle is given by u^{i}(x,y)=\left({\begin{array}[]{*{18}{c}}1&x&y\end{array}}\right){\left({\begin{array}[]{*{17}{c}}1&{{x_{1}^{i}}}&{{y_{1}^{i}}}\\ 1&{{x_{2}^{i}}}&{{y_{2}^{i}}}\\ 1&{{x_{3}^{i}}}&{{y_{3}^{i}}}\end{array}}\right)^{-1}}\left({\begin{array}[]{*{17}{c}}{{u_{1}^{i}}}\\ {{u_{2}^{i}}}\\ {{u_{3}^{i}}}\end{array}}\right), where ${u^{i}_{l}}$ is the dependent variable at the $l$th vertex in the $i$th 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}(\bm{r}_{p})=\alpha_{(i,{\bf r}_{p})}^{l}u_{l}^{i}$, where $\bm{r}_{p}$ represents the coordinate of the $p$th quadrature point, and $\alpha_{(i,{\bf r}_{p})}^{l}$ denotes the $l$th interpolating function corresponding to $u^{i}_{l}$, and $l\in(1,2,3)$. We follow Einstein summation convention for the indices of $l$, which are repeated both in the superscript of $\alpha$ and subscript of $u$.

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

[TABLE]

where Einstein summation convention is applied for the indices of $l$ and $\bm{r}_{p}$, and $(i,\bm{r}_{p})$ denotes the quadrature point ($\bm{r}_{p}\in\left(\bm{r}_{1},\bm{r}_{2},\bm{r}_{3}\right)$) in the $i$th triangle, the label of interpolating function $l\in\left({1,2,3}\right)$, $i/j$ indexes the triangle and does not follow Einstein summation rule. In Eq. (6), the 3 auxiliary variables in triangle $i$ on the left hand side equal to the sum of three terms from the right hand side: the first term is the incident field, the second term is the self-term contribution from the same triangle $i$, the last term is the sum of the contribution from all the other triangles. $A_{j}$ is the area of the triangle $j$. Applying Eq. (6) to the simple scatter 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,

[TABLE]

where $\bm{U}_{i}=\left[u^{i}_{1},u^{i}_{2},u^{i}_{3}\right]^{T}$ is a $3\times 1$ column vector, i.e., the 3 dependent variables for triangle $i$, and $\bar{\bm{\alpha}}_{ii}$ and $\bm{\bar{g}}_{ij}$ are $3\times 3$ matrices as given by Eq. (6), explicitly \bar{\bm{\alpha}}_{ii}=\left({\begin{array}[]{*{16}{c}}a_{(i,{\bm{r}_{1}})}^{1}&a_{(i,{\bm{r}_{1}})}^{2}&{a_{(i,{\bm{r}_{1}})}^{3}}\\ a_{(i,{\bm{r}_{2}})}^{1}&a_{(i,{\bm{r}_{2}})}^{2}&{a_{(i,{\bm{r}_{2}})}^{3}}\\ a_{(i,{\bm{r}_{3}})}^{1}&a_{(i,{\bm{r}_{3}})}^{2}&{a_{(i,{\bm{r}_{3}})}^{3}}\end{array}}\right), and \bm{\bar{g}}_{ij}=k_{0}^{2}\Delta\varepsilon\frac{{{A_{j}}}}{3}\sum\limits_{\bm{r}_{l}\in({\bm{r}_{1},\bm{r}_{2},\bm{r}_{3}})}\left({\begin{array}[]{*{16}{c}}g_{(i,{\bm{r}_{1}})}^{(j,{\bm{r}_{l}})}&g_{(i,{\bm{r}_{1}})}^{(j,{\bm{r}_{l}})}&g_{(i,{\bm{r}_{1}})}^{(j,{\bm{r}_{l}})}\\ g_{(i,{\bm{r}_{2}})}^{(j,{\bm{r}_{l}})}&g_{(i,{\bm{r}_{2}})}^{(j,{\bm{r}_{l}})}&{g_{(i,{\bm{r}_{2}})}^{(j,{\bm{r}_{l}})}}\\ g_{(i,{\bm{r}_{3}})}^{(j,{\bm{r}_{l}})}&g_{(i,{\bm{r}_{3}})}^{(j,{\bm{r}_{l}})}&{g_{(i,{\bm{r}_{3}})}^{(j,{\bm{r}_{l}})}}\end{array}}\right)\left({\begin{array}[]{*{16}{c}}a_{(j,{\bm{r}_{l}})}^{1}&a_{(j,{r_{l}})}^{2}&{a_{(j,{\bm{r}_{l}})}^{3}}\\ a_{(j,{\bm{r}_{l}})}^{1}&a_{(j,{r_{l}})}^{2}&{a_{(j,{\bm{r}_{l}})}^{3}}\\ a_{(j,{\bm{r}_{l}})}^{1}&a_{(j,{r_{l}})}^{2}&{a_{(j,{\bm{r}_{l}})}^{3}}\end{array}}\right).

In the last, the matrix form of GIE in Eq. (7) is obtained simply by stacking the triangles one by one, thus the shared vertices among the neighbour triangles are not taken into account. Therefore, $\bm{U}=[\bm{U}_{1},\bm{U}_{2},\bm{U}_{3},\bm{U}_{4}]^{T}$ is a $12\times 1$ column vector, in contrast to 6 DOFs in total shown in Fig. 1 (c). By introducing a global indexing of dependent variables ${\bm{U}}^{g}$ as well as the assembling transformation matrix $\bm{\bar{T}}_{asm}$, i.e, $\bm{U}^{g}=\bm{\bar{T}}_{asm}\bm{U}$, the Eq. (7) denoted as $\left[\bm{\bar{\alpha}}_{ii}\right]\bm{U}=\bm{E}^{inc}+\left[\bm{\bar{g}}_{ij}\right]\bm{U}$, can be reformulated into the final matrix form that can be readily solved using the standard linear algebra,

[TABLE]

where $\left[\bm{\bar{\alpha}}_{ii}\right]$ ($\left[\bm{\bar{g}}_{ij}\right]$) is the full matrix of the block matrix $\bm{\bar{\alpha}}_{ii}$ ( $\bm{\bar{g}}_{ij}$). The assembling transformation matrix $\bm{\bar{T}}_{asm}$ 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 $\bm{\bar{T}}_{asm}$ equals to the total number of elements of the triangulation matrix ${\bm{E}}_{tri}$, while the row number equals to the total number of vertices. By flattening ${\bm{E}}_{tri}$ row by row, one can construct $\bm{\bar{T}}_{asm}$ by mapping the row to the element of the flattened ${\bm{E}}_{tri}$, with the location of nonzero element in the same row corresponding to the value of that element in the flattened ${\bm{E}}_{tri}$.

3 Results and discussions

We examine the far-field pattern of 2D light scattering [51] using our proposed method, i.e., type-1 FEGIE and type-2 FEGIE, in comparison with stair-case GIE [32, 9] and Mie theory. In the first example, we study light scattering by a dielectric rod (infinitely long along z-axis, $\epsilon_{r}$=4), with the background being air. The incident light is TE polarized, propagating along x-axis with vacuum wavelength 0.633 $\mu$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 FIE (magenta dash-dot line). As for the scatter 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. Due to improved geometric approximation of the scatter, it is expected our proposed FEGIE is more accurate than that of staircase GIE. Indeed, the farfield pattern calculated from our type-1 FEGIE and type-2 FEGIE shows more accurate results, compared with staircase GIE, where the number of triangles in type-1 FEGIE and type-2 FEGIE are the same, i.e., 1080, and the number of squares in staircase GIE is 1085. A close examination 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, due to the fact that the finite element integration technique in each element, i.e., quadrature rule, is introduced in type-2 FEGIE to improve the numerical integration. As a side remark, the numbers 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.

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_{1}$=0.3 $\mu$m, $r_{2}$=0.7 $\mu$m, $r_{3}$=1.0 $\mu$m, and the dielectric constants of the three layers are given as $\epsilon_{1}=3$, $\epsilon_{2}=4$, $\epsilon_{3}=5$. The radii of the 5-layered rod [52], as sketched in the inset of (b), are given by $r_{1}$=0.0158 $\mu$m, $r_{2}$=0.0285 $\mu$m, $r_{3}$=0.0475 $\mu$m, $r_{2}$=0.0601 $\mu$m, $r_{5}$=0.633 $\mu$m. The corresponding dielectric constants are given as $\varepsilon_{1}=0.9814$, $\varepsilon_{2}=4.7676$, $\varepsilon_{3}=1.7075$, $\epsilon_{4}=1.4711$, $\varepsilon_{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 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 [52], 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.

As the third example, we study the far-field pattern of a geometrically complicated scatter 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). The scatter with corrugated surface discretized with triangle and square elements are shown in Fig. 4(b) and (c), which have comparable number of elements, i.e., 2516 triangles versus 2502 squares. Fig. 4 (a) shows a few side lobes of the far-field, as indicated by the dash 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 with that 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 2502, 1370, 2516 and 103561.

Lastly, Fig. 5 shows the convergence of the farfield plotted in Fig. 4 (a) at $\theta$=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. As a final remark, the accuracy of the staircase GIE used in this paper can be improved by using averaged dielectric constant. The averaged dielectric constant is not implemented in this paper due to the complications of calculating the differences of approximated squares with the true geometry, especially for complex structures. 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 the current paper, is certainly interesting to investigate in the future.

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 scatters and a more accurate numerical integral. We have calculated the farfield patterns of several concrete scatters 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.

Acknowledge

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

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\left({{\bf{r}},{\bf{r^{\prime}}}}\right)=\frac{1}{{4i}}H_{0}^{\left(2\right)}\left({{k_{0}}{n_{ref}}\left|{{\bf{r}}-{\bf{r^{\prime}}}}\right|}\right)$. As $i=j$ in Eq. (5), $g_{ii}$ represents the interaction of a discrete element with itself, it can be approximated with

[TABLE]

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

[TABLE]

For far field, since $\left|{{\bf{r^{\prime}}}}\right|/\left|{\bf{r}}\right|\ll 1$, the Green function can be approximated with

[TABLE]

Then the far field can be obtained as

[TABLE]

Funding

Natural National Science Foundation (NSFC) (11874026, 61735006, 61775063); National Key Research and Development Program of China (2017YFA0305200); Research Grants Council of Hong Kong SAR (CityU 21302018).