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

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


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][4][5][6][7][8][9][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][12][13][14][15].Moreover, high-order multipoles such as electric and magnetic octupoles can be structurally engineered to tune optical resonances, absorption and scattering [16][17][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 phenomenon such as ultradirectional scattering, perfect reflection and transmission, anapole effects [19][20][21][22][23][24]; (2) for designing photonic devices such as metasurfaces and plasmonic arrays [25][26][27][28][29][30][31].Several routes are feasible for realizing numerical projection of electromagnetic multipoles, as is well documented in recent works [32][33][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 Section 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 Section 3, we validate the numerical projection procedure using benchmarks of an isotropic dielectric nanosphere, a symmetric scatterer, and an anisotropic nanosphere.
Finally, Section 4 summarizes the paper.

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.

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   =  0  (⋅−) , where  is the angular frequency,  0 is the complex amplitude of electric field, and  is the wave vector in the background material.The total electric/magnetic field is denoted by /, which is the sum of the incident field   /  and scattered field   /  , i.e.,  =   +   and  =   +   .In the source-free region, the wave equations for total electric and magnetic field are written as follows: (1) The wave vector  = √ ⋅  =  √  , 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  and  of order (, ) are expressed as follows [37]: In physics,   and   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   can be expanded in terms of the electric and magnetic multipole modes, where   and   are electric/magnetic multipole coefficients, which represent the contributions of electric/magnetic multipoles to the scattered field.The constant coefficient   is given by: . Considering the orthogonality of vector spherical harmonics, the multipole coefficients   and   are computed as follows [37]: 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   is computed by full-wave simulation, and the integration surface should not go beyond the simulation region.

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: where   is the permittivity distribution in the whole space with presence of the scatterer,   is the permittivity of background material, and () is the total electric field.Based on Maxwell equations, the relations between the scattered field   /  and the induced current density   are expressed as follows: The wave equations for the scattered field   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 (, ) [38], where Π  () =   () 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, where is the vacuum impedance.The two integrals in Eq. ( 9) are evaluated over the volume  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.

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 fullwave simulation.

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  with a unit volume, the numerical quadrature of a function  over the tetrahedron  is in the following form [39], where    ( = 1,2,3,4) are the natural coordinates of the  th quadrature point in the tetrahedron, and   is the corresponding weight of the th quadrature point.According to the definition of natural coordinates [40], the constraint for the quadruplet ( ) 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  of the quadrature method.For the integration over a tetrahedron, Gaussian quadrature formulas of degree  = 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  =  2 +  2 +  2 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  = 4 for the computation, resulting in an integral value of

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  2 which is described by  2 +  2 +  2 = 1, the integration of a function  on  2 is approximated using Lebedev quadrature in the following form [43], ()d ≈  1 ∑   [43].To determine the corresponding weights of quadrature points in Lebedev quadrature of order , the quadrature method must exactly integrate over all spherical harmonics   of order  ≤ .
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   ,...,  and quadrature points for order 9 ≤  ≤ 17 can be found in Ref. [42].
Without loss of generality, we demonstrate an example of accurate integration for a function  =  2 +  3 −  4 , 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

Results and discussion
In this section, we study three scattering objects, including an isotropic dielectric nanosphere, a symmetric scatterer with  3ℎ 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].

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. 5(a).
The center of the nanosphere is located at the origin of the coordinate system.The nanosphere has a refraction index of  = 3.5 and a radius of  = 210 nm.The incident light is a plane wave propagating along the  axis with linear polarization along the  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   and   , the total scattering cross-section  sca is determined by the following expression, where  sca is normalized based on the cross-sectional area π 2 of the sphere.As shown in Fig. 5(a), 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. 5(b), 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.

Scattering by a symmetric scatterer with 𝑫 𝟑𝒉 symmetry
The second benchmark evaluates the multipole coefficients of eigenmodes of a symmetric scatterer with  3ℎ 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  3ℎ group, and each eigenmode belongs to one of these representations.
Specifically, we focus on two eigenmodes belonging to different representations  1 ″ and  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  1 ″ and  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  3ℎ group can be found in Ref. [47].The multipole coefficients for the two eigenmodes in Fig. 7  As shown in Fig. 8, we evaluate the expression  points.The surface integration-based numerical projection is adopted.

Scattering by an anisotropic nanosphere
Lastly, we use an anisotropic nanosphere to verify the accuracy of our numerical projection

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.

Figure 1 .
Figure 1.Schematic of multipole decomposition.The leftmost field pattern represents the total scattered field, which can be irregular and complex.The field patterns on the right show the multipole modes of order 1,2,3 and 4.

Figure 2 .
Figure 2. Working flow of the algorithm for multipole decomposition.

0. 485352892045444
up to 16 digits of accuracy, which precisely matches the analytical solution of

Figure 3 .
Figure 3. Geometry of the integration region and selected Gaussian quadrature points.

Figure 4 .
Figure 4. Geometry of the integration region and selected Lebedev quadrature points.

Figure 5 .
Figure 5.The multipole scattering cross-section of the dielectric nanosphere as a function of the incident wavelength.The analytical curves are derived from Mie theory, while the

Figure 6 .
Figure 6.Convergence versus the number of quadrature points.Subplot (a) shows the convergence of magnetic octupole at 818 nm using surface integration.Subplot (b) shows the convergence of magnetic quadrupole at 1000 nm using volume integration.Black dashed line shows the analytical value from Mie theory.
Fig. 7(c), the mode profile is unchanged under the operation 3  , indicating that this eigenmode satisfies the requirement of the  2 ″ representation.The corresponding constraints for the  2 ″ representation require  3,3 =  3,−3 and  4,3 = − 4,−3 , which also agree with our numerical results as shown in Fig. 7(d).

Figure 7 .
Figure 7. Two eigenmodes of the  3ℎ scatterer and the normalized (with respect to itself) multipole coefficients of each eigenmode.The eigenmode shown in subplot (a)/(c) belongs to  1 ″ / 2 ″ representation, whose multipole coefficients are presented in subplot (b)/(d).

Figure 8 . 3 𝑏 3 ,− 3 of the 𝐴 1 ″
Figure 8. Convergence of  3,−3 + 3,3  3,−3 of the  1 ″ eigenmode versus the number of quadrature 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  axis.The anisotropic dielectric tensor of the nanosphere is ‾  = (   i 0 −i   0 0 0  ), where   = 5.2,  = 3, and i = √−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.

Figure 9 .
Figure 9. Scattering cross-sections of an anisotropic nanosphere.The incident light of subplot (a) and (b) is LCP and RCP light respectively.The triangles represent the results calculated by COMSOL, and the circles represent the results from the numerical projection algorithm.
Because of the symmetry of tetrahedron, all possible permutations of the

Table 1 .
The natural/Cartesian coordinates and weights of Gaussian quadrature points of degree  = 4 for the computation of integrals over a regular tetrahedron.
,   ,   and   are the corresponding weights.To be more specific, for each group of quadrature points with the same , they are invariant under the octahedral rotation group with inversion  8 * and share the same weight

Table 2 .
The Cartesian coordinates and weights of Lebedev quadrature points of order  = 7 on a unit sphere.

Table 3 .
The character table of  3ℎ group representation  1 ″ and  2 ″ , and constraints on multipoles coefficients of the eigenmodes of two representations. is an integer.