Efficient and accurate numerical-projection of electromagnetic multipoles for scattering objects

Abstract

In this paper, we develop an efficient and accurate procedure of electromagnetic multipole decomposition by using the Lebedev and Gaussian quadrature methods to perform the numerical integration. Firstly, we briefly review the principles of multipole decomposition, highlighting two numerical projection methods including surface and volume integration. Secondly, we discuss the Lebedev and Gaussian quadrature methods, provide a detailed recipe to select the quadrature points and the corresponding weighting factor, and illustrate the integration accuracy and numerical efficiency (that is, with very few sampling points) using a unit sphere surface and regular tetrahedron. In the demonstrations of an isotropic dielectric nanosphere, a symmetric scatterer, and an anisotropic nanosphere, we perform multipole decomposition and validate our numerical projection procedure. The obtained results from our procedure are all consistent with those from Mie theory, symmetry constraints, and finite element simulations.

Graphical abstract

Keywords

Multipole decomposition / Numerical quadrature / Light scattering

1 Introduction

Electromagnetic multipoles play an indispensable role across different sub-branches of optics and photonics [1]. In classical electrodynamics, the field excited by localized spatial charge and current distribution can be decomposed into electric and magnetic multipole fields of all orders [2]. A visualization of field decomposition is presented in Fig. 1. The total scattered field can be decomposed into a series of multipole modes of all orders, including dipole, quadrupole, octupole, and high-order multipoles. In optics, the multipole decomposition method has been widely used in many nanophotonics scenarios for the analysis of light scattering of both a single scatterer and of periodic arrays of nanostructures [3–10]. The electric and magnetic multipoles excited in the scatterer contain valuable information about its optical response, including resonance, far-field patterns, and scattering cross-section [11–15]. Moreover, high-order multipoles such as electric and magnetic octupoles can be structurally engineered to tune optical resonances, absorption and scattering [16–18].

Apart from symmetric scatterers such as spheres or cylinders, normally there are no analytical solutions to the electromagnetic multipoles of irregular scatterers. Hence, numerical projection of electromagnetic multipoles is desirable for two purposes: (1) for understanding the underlying principles of many exotic phenomena such as ultra-directional scattering, perfect reflection and transmission, anapole effects [19–24]; (2) for designing photonic devices such as metasurfaces and plasmonic arrays [25–31]. Several routes are feasible for realizing numerical projection of electromagnetic multipoles, as is well documented in recent works [32–34]. For instance, Alaee et al. proposed a multipole expansion method based on exact multipole moments beyond long-wavelength approximation [32]. Evlyukhin et al. studied discrete dipole approximation to calculate the multipole components of non-spherical scatterers [33]. Evlyukhin et al. also combined analytical and numerical methods to calculate the multipoles of scatterers with anisotropic optical properties [34]. The approach for multipole decomposition considered in this paper is based on multipole modes in spherical coordinates, rather than multipole moments in Cartesian coordinates [16, 17, 23, 24, 35]. Evidently, numerical integration is one of the key ingredients to perform multipole decomposition, and can be carried out using surface or volume integral techniques, provided that the scattered field is known from full-wave simulation. Despite the relevance of numerical integration in multipole decomposition, previous studies provide no systematic treatment of how the numerical integration is carried out, including the selection of sampling points, numerical accuracy and efficiency.

In this paper, we propose an efficient and accurate numerical projection procedure for multipole decomposition, which is complementary to the existing works in multipole projection. To the best of our knowledge, the Lebedev and Gaussian quadrature methods are introduced for the first time to process the numerical projection of multipoles, wherein the integration sampling points can be selected efficiently and accurately based on the solid framework given by the quadrature rule, i.e., in accordance with the interpolation functions used in full-wave simulations such as finite element method. Moreover, high-order electric and magnetic multipoles up to the 8th order are accurately coded in our program. The codes for MATLAB implementation of our new numerical projection procedure are open-source and available on GitHub [36].

The paper is organized as follows. In Sect. 2, we review the general approach to multipole decomposition for scattering problems. Then the Lebedev and Gaussian quadrature methods for surface and volume integration are introduced and tested in detail. In Sect. 3, we validate the numerical projection procedure using benchmarks of an isotropic dielectric nanosphere, a symmetric scatterer, and an anisotropic nanosphere. Finally, Sect. 4 summarizes the paper.

2 Theory

In this section, we briefly review the principles of the multipole decomposition method. Two types of numerical projection methods are highlighted, including surface integration based on the scattered field and volume integration based on the induced current density. To perform surface and volume integration in the numerical projection procedure, the Lebedev and Gaussian quadrature methods are adopted and discussed in detail.

The general working flow of the algorithm, following the common procedure of multipole decomposition based on numerical projection, is depicted in Fig. 2. Firstly, for a given structure, the scattered and total field in the whole space is computed using full-wave simulation such as finite element method (FEM) or finite-difference time-domain (FDTD) method. Based on the field data, the multipole coefficients of the scattered field are obtained by numerical projection. There are two approaches for numerical projection, each suitable for different situations. The first approach is surface integration based on the scattered field outside the scatterer, which is mainly used for the multipole analysis of a single structure. The second approach is volume integration based on the induced current density inside the scatterer, which is applicable to the building blocks in periodic structures such as metasurfaces and photonic crystals. Examples show that the surface integration method has faster convergence than the volume integration method with less computational resources in practice.

2.1 Surface integration based on the scattered field

The first approach to numerical projection of multipoles is based on the scattered field outside the scatterer. Consider a monochromatic incident plane wave E_i = E₀e^{i(k·x−ωt)}, where ω is the angular frequency, E₀ is the complex amplitude of electric field, and k is the wave vector in the background material. The total electric/magnetic field is denoted by E/H, which is the sum of the incident field E_i/H_i and scattered field E_s/H_s, i.e., E = E_i + E_s and H = H_i + H_s. In the source-free region, the wave equations for total electric and magnetic field are written as follows:

\begin{array}{}\text{(1)}& \begin{array}{l}{\nabla }^{2}E+k^{2}E=\mathrm{0,}\\ {\nabla }^{2}H+k^{2}H=0.\end{array}\end{array}

The wave vector k=\sqrt{k\cdot k}=\omega \sqrt{\epsilon \mu }, where ε and μ are the permittivity and permeability of isotropic background. Adapting the spherical coordinate system, the eigen solutions of wave equations Eq. (1) are vector spherical harmonics. The concrete forms of vector spherical harmonics M and N of order (l, m) are expressed as follows [37]:

\begin{array}{}\text{(2)}& \begin{array}{ll}{M}_{lm}=\hfill & \left[{\text{iπ}}_{lm}(\theta ){\widehat{e}}_{\theta }-{\tau }_{lm}(\theta ){\widehat{e}}_{\varphi }\right]z_l(kr){\text{e}}^{\text{i}m\varphi }\mathrm{,}\hfill \\ N_{lm}=\hfill & l(l+1)P_l^m(\text{cos}\theta )\frac{z_l(kr)}{kr}{\text{e}}^{\text{i}m\varphi }{\widehat{e}}_r\hfill \\ \hfill & +\left[{\tau }_{lm}(\theta ){\widehat{e}}_{\theta }+{\text{iπ}}_{lm}(\theta ){\widehat{e}}_{\varphi }\right]\frac{{\left[kr{z}_l(kr)\right]}^{\prime }}{kr}{\text{e}}^{\text{i}m\varphi }\mathrm{,}\hfill \end{array}\end{array}

where P_l^m (cosθ) is the associated Legendre polynomial, z_l(kr) is the spherical Bessel function, and [krz_l(kr)]′ denotes the differentiation with respect to the argument kr. The auxiliary functions τ_lm(θ) and π_lm(θ) are defined as:

\begin{array}{}\text{(3)}& \begin{array}{c}{\tau }_{lm}(\theta )=\frac{\text{d}}{\text{d}\theta }{P}_l^m(\text{cos}\theta )\mathrm{,}\\ {\text{π}}_{lm}(\theta )=\frac{m}{\sin \theta }{P}_l^m(\text{cos}\theta )\mathrm{.}\end{array}\end{array}

In physics, N_lm and M_lm represent the field excited by electric and magnetic multipoles, respectively, located at the origin in the spherical coordinate system. Because of the completeness of vector spherical harmonics, the scattered field E_s can be expanded in terms of the electric and magnetic multipole modes,

\begin{array}{}\text{(4)}& E_s={\displaystyle {\sum }_{l=1}^{\infty }}{\displaystyle {\sum }_{m=-l}^l}{E}_{lm}\left(a_{lm}{N}_{lm}+b_{lm}{M}_{lm}\right)\mathrm{,}\end{array}

where a_lm and b_lm are electric/magnetic multipole coefficients, which represent the contributions of electric/magnetic multipoles to the scattered field. The constant coefficient E_lm is given by: E_{lm}=\frac{i^{l+1}(2l+1)}{2}\sqrt{\frac{(l-m)!}{l(l+1)(l+m)!}}. Considering the orthogonality of vector spherical harmonics, the multipole coefficients a_lm and b_lm are computed as follows [37]:

\begin{array}{}\text{(5)}& \begin{array}{l}{a}_{lm}=\frac{{\displaystyle \iint }{N}_{lm}^{*}\cdot E_s\text{d}s}{E_{lm}{\displaystyle \iint }{N}_{lm}^{*}\cdot N_{lm}\text{d}s}\mathrm{,}\\ b_{lm}=\frac{{\displaystyle \iint }{M}_{lm}^{*}\cdot E_s\text{d}s}{E_{lm}{\displaystyle \iint }{M}_{lm}^{*}\cdot M_{lm}\text{d}s}\mathrm{,}\end{array}\end{array}

where the asterisk * denotes complex conjugation. The surface integration in Eq. (5) is taken over a spherical surface which encloses the scattering object. In practice, the scattered field E_s is computed by full-wave simulation, and the integration surface should not go beyond the simulation region.

2.2 Volume integration based on the induced current density

The second approach to numerical projection of multipoles is based on the induced current density inside the scatterer, which is defined as follows:

\begin{array}{}\text{(6)}& J_s(r)=-\text{i}\omega {\epsilon }_0({\epsilon }_r-{\epsilon }_b)E(r)\mathrm{,}\end{array}

where ε_r is the permittivity distribution in the whole space with presence of the scatterer, {\epsilon }_b is the permittivity of background material, and E(r) is the total electric field. Based on Maxwell equations, the relations between the scattered field E_s/H_s and the induced current density J_s are expressed as follows:

\begin{array}{}\text{(7)}& \begin{array}{l}\nabla \cdot E_s=-\frac{\text{i}\eta }{k}\nabla \cdot J_s\mathrm{,}\\ \nabla \cdot H_s=\mathrm{0,}\\ \nabla \times E_s=\text{i}k\eta H_s\mathrm{,}\\ \nabla \times H_s=-\frac{\text{i}k}{\eta }{E}_s+J_s\mathrm{.}\end{array}\end{array}

The wave equations for the scattered field E_s can be derived based on Eq. (7). Then one can further determine the exact forms of the induced current density of electric and magnetic multipoles of order (l, m) [38],

\begin{array}{}\text{(8)}& \begin{array}{ll}{S}_{lm}=\hfill & \left[{\text{Π}}_l(kr)+{\text{Π}}_l^{″}(kr)\right]P_l^m(\text{cos}\theta ){\text{e}}^{-\text{i}m\varphi }{\widehat{e}}_r\hfill \\ \hfill & +\left[{\tau }_{lm}(\theta ){\widehat{e}}_{\theta }-{\text{iπ}}_{lm}(\theta ){\widehat{e}}_{\varphi }\right]\frac{{\text{Π}}_l^{\prime }(kr)}{kr}{\text{e}}^{-\text{i}m\varphi }\mathrm{,}\hfill \\ T_{lm}=\hfill & \left[{\text{iπ}}_{lm}(\theta ){\widehat{e}}_{\theta }+{\tau }_{lm}(\theta ){\widehat{e}}_{\varphi }\right]{\text{e}}^{-\text{i}m\varphi }{J}_l(kr)\mathrm{,}\hfill \end{array}\end{array}

where Π_l(kr) = krJ_l(kr) is the Riccati-Bessel function. As in the case of the multipole decomposition based on the scattered field, one can also calculate the multipole coefficients based on the induced current density in the following form,

\begin{array}{}\text{(9)}& \begin{array}{l}{a}_{lm}=\frac{{(-\text{i})}^{l-1}{k}^2\eta Q_{lm}}{\sqrt{\text{π}(2l+1)}}{\displaystyle {\iiint }_V}{J}_s(r)\cdot S_{lm}\text{d}\text{Ω}\mathrm{,}\\ b_{lm}=\frac{{(-\text{i})}^{l+1}{k}^2\eta Q_{lm}}{\sqrt{\text{π}(2l+1)}}{\displaystyle {\iiint }_V}{J}_s(r)\cdot T_{lm}\text{d}\text{Ω}\mathrm{,}\end{array}\end{array}

where Q_{lm}=\frac{1}{\sqrt{l(l+1)}}\sqrt{\frac{2l+1}{4\text{π}}\frac{(l-m)!}{(l+m)!}} is a constant coefficient and \eta =\sqrt{\frac{{\mu }_0}{{\epsilon }_0}} is the vacuum impedance. The two integrals in Eq. (9) are evaluated over the volume V of the scattering object. The volume integration method applies to either a single scatterer or the building blocks in a periodic array. Since the unit building block expands in the entire full-wave simulation region in the direction with periodic boundary conditions, only the volume integration method, which extracts the induced current density inside the unit building block, can be used to calculate the multipole components.

It is worth mentioning that the multipole coefficients obtained using either the scattered field-based or induced current density-based numerical projections should be identical for a given scatterer. The expressions for the two methods can be transformed into each other using the relationship between the source current and far-field [38].

As a side remark, if the center of mass of the scatterer is not the same as the coordinate origin, then the multipole decomposition is complicated. Notably, the numerical projection procedures, i.e., surface or volume integration-based numerical projection as used in this paper, still apply and the resulting multipole components are not the same. However, the observable physical quantity obtained by the combination of different kinds of multipoles should be the same. If multipole decomposition is performed with respect to the symmetric center of the scatterer, there are fewer multipole components. For example, for a spherical scatterer, if the sphere center and the coordinate origin are identical, the projected multipole components have the fewest terms. In contrast, if the sphere center has an off-set from the coordinate origin, one can still do multipole decomposition. In such a case, there are more multipole components than is the case in the overlapped scenario, but the observable physical quantities are the same. For calculation of multipoles of asymmetrical or compound scatterers, it is beneficial to choose the center of mass as the reference point to get minimum number of multipole components.

2.3 Numerical implementations of integration

To compute the multipole coefficients accurately, certain numerical quadrature rules are needed to implement the surface or volume integration described in the above sections. In order to achieve higher computational accuracy with fewer function evaluations, we first adopt the Gaussian and Lebedev quadrature methods to compute the volume and surface integration in the numerical projection, respectively.

The Gaussian quadrature method is widely used in the numerical integration in the practical implementation of electromagnetic finite element computation. As such, the integration sampling points as well as the interpolation functions used in finite element full-wave simulation are also carefully selected according to the quadrature rule and can be directly used in our proposed method for multipole projection, wherein the integration sampling points can be cut down to the minimum, as constrained by finite element full-wave simulation.

2.3.1 Gaussian quadrature method for volume integration

Since the boundary of the integration region can be irregular, finite elements are often used to mesh the integration region. The most commonly used shape for meshing is the tetrahedron in three-dimensional space and the triangle for two-dimensional surfaces. To compute the integral over a tetrahedron, the Gaussian quadrature method is adopted due to its high accuracy. For a single tetrahedron V with a unit volume, the numerical quadrature of a function f over the tetrahedron V is in the following form [39],

\begin{array}{}\text{(10)}& {\displaystyle {\iiint }_V}f\left({\alpha }_1\mathrm{,}{\alpha }_2\mathrm{,}{\alpha }_3\mathrm{,}{\alpha }_4\right)\text{d}\text{Ω}\approx {\displaystyle {\sum }_i}{\omega }_{i}f\left({\alpha }_1^i\mathrm{,}{\alpha }_2^i\mathrm{,}{\alpha }_3^i\mathrm{,}{\alpha }_4^i\right)\mathrm{,}\end{array}

where {\alpha }_j^i (j = 1, 2, 3, 4) are the natural coordinates of the i th quadrature point in the tetrahedron, and ω_i is the corresponding weight of the i th quadrature point. According to the definition of natural coordinates [40], the constraint for the quadruplet \left({\alpha }_1^i\mathrm{,}{\alpha }_2^i\mathrm{,}{\alpha }_3^i\mathrm{,}{\alpha }_4^i\right) is {\displaystyle {\sum }_{j=1}^4}{\alpha }_j^i=1. Because of the symmetry of tetrahedron, all possible permutations of the quadruplet \left({\alpha }_1^i\mathrm{,}{\alpha }_2^i\mathrm{,}{\alpha }_3^i\mathrm{,}{\alpha }_4^i\right) are equivalent quadrature points sharing the same weight. The coordinates of the quadrature points and weights for the quadrature expression Eq. (10) are determined by integrating polynomials of particular forms exactly, leading to a set of nonlinear equations which are relatively difficult to solve [39]. The choice of polynomial order and form for exact integration depends on the desired degree N of the quadrature method. For the integration over a tetrahedron, Gaussian quadrature formulas of degree N = 4 to 8 are given in Ref. [41].

To illustrate the Gaussian quadrature intuitively, we present an example of computing the integration of the function f = x² + y² + z² over a regular tetrahedron with a unit volume with only 11 sampling points. The geometry of the tetrahedron and selected quadrature points are shown in Fig. 3. The exact coordinates of these quadrature points and the corresponding weights are listed in Table 1. We utilize the Gaussian quadrature method with degree N = 4 for the computation, resulting in an integral value of 0.485352892045444 up to 16 digits of accuracy, which precisely matches the analytical solution of \frac{7}{10\times \sqrt[3]{3}}=0.485352892045444.

2.3.2 Lebedev quadrature method for surface integration

As mentioned above, Gaussian quadrature is compatible with finite elements such as tetrahedrons or triangles. It is straightforward to adapt Gaussian quadrature to compute the surface integration over a sphere, as long as the quadrature points and weights are computed based on the symmetry of triangles. However, note that a sphere has a higher degree of symmetry than a triangle, resulting in more quadrature points that are equivalent to each other. For instance, quadrature points which are equivalent under the octahedron rotation and inversion on the sphere share the same weight [42]. Therefore, a smaller number of independent quadrature points can achieve the same accuracy as Gaussian quadrature. This efficient quadrature method for a sphere is called Lebedev quadrature, which is adopted in this paper.

Given a unit sphere S² which is described by x² + y² + z² = 1, the integration of a function f on S² is approximated using Lebedev quadrature in the following form [43],

\begin{array}{}\text{(11)}& \begin{array}{ll}I=\hfill & {\displaystyle {\iint }_{S^2}}f(r)\text{d}\sigma \approx A_1{\displaystyle \sum _{i=1}^6}f(a_i^1)\hfill \\ \hfill & +A_2{\displaystyle \sum _{i=1}^{12}}f(a_i^2)+A_3{\displaystyle \sum _{i=1}^8}f(a_i^3)\hfill \\ \hfill & +{\displaystyle \sum _{k=1}^{N_1}}{B}_k{\displaystyle \sum _{i=1}^{24}}f(b_i^k)+{\displaystyle \sum _{k=1}^{N_2}}{C}_k{\displaystyle \sum _{i=1}^{24}}f(c_i^k)\hfill \\ \hfill & +{\displaystyle \sum _{k=1}^{N_3}}{D}_k{\displaystyle \sum _{i=1}^{48}}f(d_i^k)\mathrm{,}\hfill \end{array}\end{array}

where the quadrature points a_i^k, b_i^k, c_i^k and d_i^k are determined according to the octahedronal symmetry, and A_k, B_k, C_k and D_k are the corresponding weights. To be more specific, for each group of quadrature points with the same k, they are invariant under the octahedral rotation group with inversion G_8^{*} and share the same weight [43]. To determine the corresponding weights of quadrature points in Lebedev quadrature of order n, the quadrature method must exactly integrate over all spherical harmonics Y_lm of order l ≤ n. In addition, it is sufficient to determine the quadrature weights by requiring the exact integration of the polynomials that are invariant under the octahedral rotation with inversion. As in the case of the determination of Gaussian quadrature, the exact integration of polynomials or spherical harmonics also leads to a set of nonlinear equations to be solved. The computed value of weights A_k, …, D_k and quadrature points for order 9 ≤ n ≤ 17 can be found in Ref. [42].

Without loss of generality, we demonstrate an example of accurate integration for a function f = x² + y³ − z⁴, which is randomly chosen here, over a unit sphere using the Lebedev quadrature method with only 26 sampling points. The sampling points are shown in Fig. 4, and the coordinates and the corresponding weights are also listed in Table 2. Remarkably, the computed integral yields a value of 1.675516081914556, which shows agreement with the analytical solution of \frac{8\text{π}}{15}=1.675516081914556, up to 16 digits of accuracy.

3 Results and discussion

In this section, we study three scattering objects, including an isotropic dielectric nanosphere, a symmetric scatterer with D_3h symmetry, and an anisotropic nanosphere, to validate our numerical projection procedure for multipole decomposition. The scattered field and induced current density are computed by full-wave simulation using COMSOL Multiphysics [44].

3.1 Scattering by an isotropic dielectric nanosphere

We compute the multipole scattering cross-section of an isotropic dielectric nanosphere as discussed in Ref. [45]. The nanosphere and the incident condition are illustrated in Fig. 5a. The center of the nanosphere is located at the origin of the coordinate system. The nanosphere has a refraction index of n = 3.5 and a radius of a = 210 nm. The incident light is a plane wave propagating along the z axis with linear polarization along the x axis. The wavelength of incident light ranges from 800 to 1900 nm. For an isotropic spherical scatterer, Mie theory provides analytical expressions for the calculation of multipole coefficients [1]. Furthermore, we utilize both the surface and volume integration methods of numerical projection to compute the multipole coefficients. Based on the multipole coefficients a_lm and b_lm, the total scattering cross-section σ_sca is determined by the following expression,

\begin{array}{}\text{(12)}& {\sigma }_{\text{sca}}=\frac{2}{k^2{a}^2}{\displaystyle {\sum }_{l=1}^{\infty }}{\displaystyle {\sum }_{m=-l}^l}(2l+1)\left({\left|a_{lm}\right|}^2+{\left|b_{lm}\right|}^2\right)\mathrm{,}\end{array}

where σ_sca is normalized based on the cross-sectional area πa² of the sphere. As shown in Fig. 5a, it is evident that the numerical results from the surface and volume integration methods of our numerical projection procedure both have excellent agreement with the analytical results from Mie theory. The two resonance peaks are characterized accurately in the wavelength range from 900 to 1900 nm. Furthermore, it can be found that within this wavelength range, only the electric/magnetic dipole and quadrupole make significant contributions to the total scattering cross-section. In contrast, as shown in Fig. 5b, the magnetic octupole has a notable impact in the wavelength range from 800 to 900 nm. In particular, the resonance peak at 818 nm is exclusively due to the magnetic octupole.

The convergence analysis of the numerical projection procedure based on both surface and volume integration is shown in Fig. 6. The convergence of magnetic octupole (MO) scattering cross-section at incident wavelength of 818 nm is calculated using surface integration, and the convergence of magnetic quadrupole (MQ) scattering cross-section at 1000 nm is calculated using volume integration. The green region in Fig. 6 shows the convergence range for each case. The convergence thresholds are set to 1% and 0.2% for subplots (a) and (b), respectively. Evidently, Lebedev quadrature has a faster convergence than the trapezoidal rule of surface integration when the number of quadrature points is higher. Similarly, Gaussian quadrature also shows a faster and more stable convergence, while the trapezoidal rule of volume integration shows oscillation near the analytical value from Mie theory. It is also evident that the surface integration method requires fewer quadrature points than are required by the volume integration method. This is because meshing the volume of scatterer needs more finite elements than meshing the scatterer's surface.

As a side remark, the distribution of Lebedev quadrature points in this demonstration is the same as in the benchmark case (a unit sphere) in Sec. 2.3.2, up to a scaling factor. This is because Lebedev quadrature points are only determined by the Lebedev quadrature rule, regardless of the geometrical shape of the scatterer inside the integration sphere. However, to calculate volume integrals, the 3D scatterer needs to be meshed using tetrahedrons, so that the quadrature points are chosen at each tetrahedron. The Gaussian quadrature rule provides the quadrature points of different order in a reference tetrahedron. To apply Gaussian quadrature to real physical tetrahedron meshes, whose geometrical shape is generally irregular, there exists a mapping from the reference tetrahedron to physical tetrahedron. Therefore, the distribution of Gaussian quadrature points in each tetrahedron mesh for scattering particle cases can be different from the benchmark case (a regular tetrahedron) in Sec. 2.3.1, and the quadrature points can be transformed to each other by a mapping.

3.2 Scattering by a symmetric scatterer with D_3h symmetry

The second benchmark evaluates the multipole coefficients of eigenmodes of a symmetric scatterer with D_3h symmetry. In optics, symmetry can be used to classify eigenmodes in structures such as waveguides or photonic crystals [46]. Furthermore, optical symmetries can introduce certain constraints on the optical properties of the structure. For symmetric scatterers, a group theory approach can be applied to analyze the constraints of multipole coefficients imposed by the symmetry of the scatterer [47, 48]. There are six irreducible representations of D_3h group, and each eigenmode belongs to one of these representations. Specifically, we focus on two eigenmodes belonging to different representations A_1^{″} and A_2^{″}. It is demonstrated that there are exact relationships between multipole coefficients of the eigenmodes. The constraints on the multipole coefficients can vary for eigenmodes belonging to different representations of the symmetry group. The character tables for the A_1^{″} and A_2^{″} repersentations, along with the corresponding constraints on the multipole coefficients of each representation are listed in Table 3. A detailed explanation of the symmetry operations in the D_3h group can be found in Ref. [47].

The multipole coefficients for the two eigenmodes in Fig. 7a and c are computed using the surface integration method of our numerical projection procedure. The results are shown in Fig. 7b and d, respectively. The scatterer, represented by the solid black lines in Fig. 7a and c, has a refractive index of n = 8 and is a three-dimensional structure, and only the xy plane is depicted. It can be checked that both eigenmodes are invariant under the operation E and 2C₃, and are reversed under the operation σ_h and 2s₃. For the eigenmode shown in Fig. 7a, the mode profile is reversed under the operation 3σ_v, which represents reflection in the xy plane. Therefore, this eigenmode satisfies the character of the A_1^{″} representation. The constraints listed in Table 3 for the A_1^{″} representation require the multipole coefficients to satisfy b_3,3 = −b_3,−3 and a_4,3 = a_4,−3, which perfectly match the computed multipole coefficients shown in Fig. 7b. Similarly, for the eigenmode in Fig. 7c, the mode profile is unchanged under the operation 3σ_v, indicating that this eigenmode satisfies the requirement of the A_2^{″} representation. The corresponding constraints for the A_2^{″} representation require b_3,3 = b_3,−3 and a_4,3 = −a_4,−3, which also agree with our numerical results as shown in Fig. 7d.

As shown in Fig. 8, we evaluate the expression \frac{b_{\mathrm{3,}-3}+b_{\mathrm{3,3}}}{b_{\mathrm{3,}-3}} of the A_1^{″} eigenmode to analyze the convergence of surface integration-based numerical projection based on this benchmark. According to the symmetry constraints for A_1^{″} eigenmodes, the expression \frac{b_{\mathrm{3,}-3}+b_{\mathrm{3,3}}}{b_{\mathrm{3,}-3}} should be exactly zero. Evidently, Lebedev quadrature also has a faster convergence than the trapezoidal rule of surface integration in this benchmark.

3.3 Scattering by an anisotropic nanosphere

Lastly, we use an anisotropic nanosphere to verify the accuracy of our numerical projection procedure on scatterers with optical anisotropy. The center of the nanosphere is located at the origin of coordinates and has a radius of 200 nm. The incident light is a plane wave with left-handed or right-handed circular polarization (LCP or RCP), and the incident direction is along the positive z axis. The anisotropic dielectric tensor of the nanosphere is \overline{{\epsilon }_r}=\left(\begin{array}{ccc}{\epsilon }_r& \text{i}g& 0\\ -\text{i}g& {\epsilon }_r& 0\\ 0& 0& {\epsilon }_r\end{array}\right), where ε_r = 5.2, g = 3, and \text{i}=\sqrt{-1}. As the incident light is LCP or RCP light, the nanosphere exhibits different scattering characteristics due to its anisotropy [34]. The computation results of the scattering cross-section under LCP and RCP illumination are shown in Fig. 9. The total scattering cross-sections obtained by numerical projection are calculated by the sum of scattering cross-sections of multipoles of each order, which are in excellent agreement with the total scattering cross-sections calculated in COMSOL.

4 Conclusion

In summary, we introduce the Lebedev and Gaussian quadrature methods in the numerical projection procedure of multipole decomposition for various scattering objects, and provide the corresponding open-source program, wherein the electric and magnetic multipoles up to the 8th order are coded. Two approaches for numerical projection, surface and volume integration, are reviewed. To validate our numerical projection procedure, we apply it to analyze the light scattering of an isotropic dielectric nanosphere, a symmetric scatterer, and an anisotropic nanosphere. The results show excellent agreement with those from Mie theory, symmetry constraints, and finite element simulation. Moreover, the Lebedev and Gaussian quadrature methods show faster convergence, than is achieved using the trapezoidal rule, for surface and volume integration in multipole decomposition.

Our proposed numerical projection procedure is highly efficient and accurate in analyzing the scattering properties of various nanostructures, and is a useful complement to existing works on multipole projection. Our numerical method is publicly available and could prove beneficial for those working in the field of nanophotonics, enabling convenient design of nanostructures with special scattering properties such as directional scattering, phase modulation, and high transmittance.

References

[1] Bohren, C.F., Huffman, D.R.: Absorption and scattering of light by small particles. John Wiley & Sons, Hoboken (1998)

[2] Jackson, J.: Classical electrodynamics. John Wiley & Sons, Hoboken (1998)

[3] Evlyukhin, A.B., Fischer, T., Reinhardt, C., Chichkov, B.N.: Optical theorem and multipole scattering of light by arbitrarily shaped nanoparticles. Phys. Rev. B 94(20), 205434 (2016)

[4] Terekhov, P.D., Babicheva, V.E., Baryshnikova, K.V., Shalin, A.S., Karabchevsky, A., Evlyukhin, A.B.: Multipole analysis of dielectric metasurfaces composed of nonspherical nanoparticles and lattice invisibility effect. Phys. Rev. B 99(4), 045424 (2019)

[5] Chen, W., Yang, Q., Chen, Y., Liu, W.: Scattering activities bounded by reciprocity and parity conservation. Phys. Rev. Res. 2(1), 013277 (2020)

[6] Kuznetsov, A.V., Valero, A.C., Shamkhi, H.K., Terekhov, P., Ni, X., Bobrovs, V., Rybin, M.V., Shalin, A.S.: Special scattering regimes for conical all-dielectric nanoparticles. Sci. Rep. 12(1), 21904 (2022)

[7] Chen, W., Yang, Q., Chen, Y., Liu, W.: Extremize optical chiralities through polarization singularities. Phys. Rev. Lett. 126(25), 253901 (2021)

[8] Alaee, R., Rockstuhl, C., Fernandez-Corbaton, I.: Exact multipolar decompositions with applications in nanophotonics. Adv. Opt. Mater. 7(1), 1800783 (2019)

[9] Evlyukhin, A.B., Chichkov, B.N.: Multipole decompositions for directional light scattering. Phys. Rev. B 100(12), 125415 (2019)

[10] Fu, Y.H., Kuznetsov, A.I., Miroshnichenko, A.E., Yu, Y.F., Luk'yanchuk, B.: Directional visible light scattering by silicon nanoparticles. Nat. Commun. 4(1), 1527 (2013)

[11] Chen, W., Chen, Y., Liu, W.: Multipolar conversion induced sub-wavelength high-Q Kerker supermodes with unidirectional radiations. Laser Photonics Rev. 13(9), 1900067 (2019)

[12] van de Haar, M.A., van de Groep, J., Brenny, B., Polman, A.: Controlling magnetic and electric dipole modes in hollow dielectric nanocylinders. Opt. Express 24(3), 2047–2064 (2016)

[13] Geffrin, J.M., García-Cámara, B., Gómez-Medina, R., Albella, P., Froufe-Pérez, L.S., Eyraud, C., Litman, A., Vaillon, R., González, F., Nieto-Vesperinas, M., Sáenz, J.J., Moreno, F.: Magnetic and electric coherence in forward- and back-scattered electromagnetic waves by a single dielectric subwavelength sphere. Nat. Commun. 3(1), 1171 (2012)

[14] Chen, W., Chen, Y., Liu, W.: Singularities and Poincaré Indices of Electromagnetic Multipoles. Phys. Rev. Lett. 122(15), 153907 (2019)

[15] Kuznetsov, A.I., Miroshnichenko, A.E., Brongersma, M.L., Kivshar, Y.S., Luk'yanchuk, B.: Optically resonant dielectric nanostructures. Science 354, 2742 (2016)

[16] Gurvitz, E.A., Ladutenko, K.S., Dergachev, P.A., Evlyukhin, A.B., Miroshnichenko, A.E., Shalin, A.S.: The high-order toroidal moments and anapole states in all-dielectric photonics. Laser Photonics Rev. 13(5), 1800266 (2019)

[17] Terekhov, P.D., Evlyukhin, A.B., Redka, D., Volkov, V.S., Shalin, A.S., Karabchevsky, A.: Magnetic octupole response of dielectric quadrumers. Laser Photonics Rev. 14(4), 1900331 (2020)

[18] Prokhorov, A.V., Terekhov, P.D., Gubin, M.Y., Shesterikov, A.V., Ni, X., Tuz, V.R., Evlyukhin, A.B.: Resonant light trapping via lattice-induced multipole coupling in symmetrical metasurfaces. ACS Photonics 9(12), 3869–3875 (2022)

[19] Liu, W., Zhang, J., Lei, B., Ma, H., Xie, W., Hu, H.: Ultra-directional forward scattering by individual core-shell nanoparticles. Opt. Express 22(13), 16178–16187 (2014)

[20] Liu, W.: Generalized Magnetic Mirrors. Phys. Rev. Lett. 119(12), 123902 (2017)

[21] Yang, Y., Miroshnichenko, A.E., Kostinski, S.V., Odit, M., Kapitanova, P., Qiu, M., Kivshar, Y.S.: Multimode directionality in all-dielectric metasurfaces. Phys. Rev. B 95(16), 165426 (2017)

[22] Liu, W., Miroshnichenko, A.E.: Beam steering with dielectric metalattices. ACS Photonics 5, 1733–1741 (2018)

[23] Miroshnichenko, A.E., Evlyukhin, A.B., Yu, Y.F., Bakker, R.M., Chipouline, A., Kuznetsov, A.I., Luk'yanchuk, B., Chichkov, B.N., Kivshar, Y.S.: Nonradiating anapole modes in dielectric nanoparticles. Nat. Commun. 6, 8069 (2015)

[24] Saadabad, R.M., Huang, L., Evlyukhin, A.B., Miroshnichenko, A.E.: Multifaceted anapole: from physics to applications [Invited]. Opt. Mater. Express 12, 1817–1837 (2022)

[25] Khorasaninejad, M., Chen, W.T., Devlin, R.C., Oh, J., Zhu, A.Y., Capasso, F.: Metalenses at visible wavelengths: diffraction-limited focusing and subwavelength resolution imaging. Science 352(6290), 1190–1194 (2016)

[26] Wang, S., Wu, P.C., Su, V.C., Lai, Y.C., Chen, M.K., Kuo, H.Y., Chen, B.H., Chen, Y.H., Huang, T.T., Wang, J.H., Lin, R.M., Kuan, C.H., Li, T., Wang, Z., Zhu, S., Tsai, D.P.: A broadband achromatic metalens in the visible. Nat. Nanotechnol. 13(3), 227–232 (2018)

[27] Yu, N., Capasso, F.: Flat optics with designer metasurfaces. Nat. Mater. 13(2), 139–150 (2014)

[28] Neshev, D., Aharonovich, I.: Optical metasurfaces: new generation building blocks for multi-functional optics. Light Sci. Appl. 7(1), 58 (2018)

[29] Yang, Y., Wang, W., Moitra, P., Kravchenko, I.I., Briggs, D.P., Valentine, J.: Dielectric meta-reflectarray for broadband linear polarization conversion and optical vortex generation. Nano Lett. 14(3), 1394–1399 (2014)

[30] Decker, M., Staude, I., Falkner, M., Dominguez, J., Neshev, D.N., Brener, I., Pertsch, T., Kivshar, Y.S.: High-efficiency dielectric Huygens' surfaces. Adv. Opt. Mater. 3(6), 813–820 (2015)

[31] Jahani, S., Jacob, Z.: All-dielectric metamaterials. Nat. Nanotechnol. 11(1), 23–36 (2016)

[32] Alaee, R., Rockstuhl, C., Fernandez-Corbaton, I.: An electromagnetic multipole expansion beyond the long-wavelength approximation. Opt. Commun. 407, 17–21 (2018)

[33] Evlyukhin, A.B., Reinhardt, C., Chichkov, B.N.: Multipole light scattering by nonspherical nanoparticles in the discrete dipole approximation. Phys. Rev. B 84(23), 235429 (2011)

[34] Evlyukhin, A.B., Tuz, V.R.: Electromagnetic scattering by arbitrary-shaped magnetic particles and multipole decomposition: Analytical and numerical approaches. Phys. Rev. B 107(15), 155425 (2023)

[35] Nanz, S.: Toroidal multipole moments in classical electrodynamics: an analysis of their emergence and physical significance. Springer, Wiesbaden (2016)

[36] Guo, W., Cai, Z., Xiong, Z., Chen, W., Chen, Y.: HUST-CPO/Multipole-Decomposition-for-Scattering. Availale at the website of github.com/HUST-CPO/Multipole-Decomposition-for-Scattering, GitHub (2023)

[37] Mühlig, S., Menzel, C., Rockstuhl, C., Lederer, F.: Multipole analysis of meta-atoms. Metamaterials 5(2), 64–73 (2011)

[38] Grahn, P., Shevchenko, A., Kaivola, M.: Electromagnetic multipole theory for optical nanomaterials. New J. Phys. 14(9), 093033 (2012)

[39] Yu, J.: Symmetric gaussian quadrature formulae for tetrahedronal regions. Comput. Methods Appl. Mech. Eng. 43(3), 349–353 (1984)

[40] Dunavant, D.A.: High degree efficient symmetrical Gaussian quadrature rules for the triangle. Int. J. Numer. Meth. Eng. 21(6), 1129–1148 (1985)

[41] Keast, P.: Moderate-degree tetrahedral quadrature formulas. Comput. Methods Appl. Mech. Eng. 55(3), 339–348 (1986)

[42] Lebedev, V.I.: Values of the nodes and weights of ninth to seventeenth order Gauss-Markov quadrature formulae invariant under the octahedron group with inversion. USSR Comput. Math. Math. Phys. 15(1), 44–51 (1975)

[43] Lebedev, V.I.: Quadratures on a sphere. USSR Comput. Math. Math. Phys. 16(2), 10–24 (1976)

[44] COMSOL Multiphysics 6.0: a finite element analysis, solver and simulation software. Available at the website of comsol.com (2023)

[45] García-Etxarri, A., Gómez-Medina, R., Froufe-Pérez, L.S., López, C., Chantada, L., Scheffold, F., Aizpurua, J., Nieto-Vesperinas, M., Sáenz, J.J.: Strong magnetic response of submicron dielectric particles in the infrared. Opt. Express 19(6), 4815–4826 (2011)

[46] Sakoda, K.: Optical properties of photonic crystals. Springer, Berlin (2005)

[47] Xiong, Z., Yang, Q., Chen, W., Wang, Z., Xu, J., Liu, W., Chen, Y.: On the constraints of electromagnetic multipoles for symmetric scatterers: eigenmode analysis. Opt. Express 28(3), 3073–3085 (2020)

[48] Poleva, M., Frizyuk, K., Baryshnikova, K., Evlyukhin, A., Petrov, M., Bogdanov, A.: Multipolar theory of bianisotropic response of meta-atoms. Phys. Rev. B 107(4), L041304 (2023)

RIGHTS & PERMISSIONS

The Author(s) 2023