Decomposition of scattered electromagnetic fields into vector spherical harmonics on surfaces with general shapes

Decomposition of scattered electromagnetic fields into vector spherical harmonics on surfaces with general shapes


Decomposing the field scattered by an object into vector spherical harmonics (VSH) is the prime task when discussing its optical properties on more analytical grounds. Thus far, it was frequently required in the decomposition that the scattered field is available on a spherical surface enclosing the scatterer; being with that adapted to the spatial dependency of the VSHs but being rather incompatible with many numerical solvers. To mitigate this problem, we propose an orthogonal expression for the decomposition that holds for any surface that encloses the scatterer, independently of its shape. We also show that the orthogonal relations remain unchanged when the radiative VSH used for the expansion of the scattered field are substituted by the VSH used for the expansion of the illumination as test functions. This is a key factor for the numerical stability of our decomposition. As example, we use a finite-element based solver to compute the multipole response of a nanorod illuminated by a plane wave and study its convergence properties.


I Introduction

The vector spherical harmonic (VSH) decomposition, or multipole expansion, is a useful tool to study electromagnetic scattering phenomena. Vector spherical harmonics are well known solutions to the time harmonic Maxwell’s equations in homogeneous media Jackson (1999); Mishchenko et al. (2002), forming a complete basis for the electromagnetic field. The field scattered by an object immersed in a homogeneous medium upon interacting with an incident field can always be decomposed into the radiative VSHs and as


where the elements of the basis correspond to the field created by a point multipole with definite properties: refers to the total angular momentum and to the angular momentum along a chosen axis. and are multipolar fields with different parity symmetry and their exact definition can be found in literature Mishchenko et al. (2002). They correspond to the electric field of electric multipoles and the electric field of magnetic multipoles, respectively. The complex amplitudes and in the expansion express then the contribution of the respective multipolar field to the total scattered field. It is the purpose of the multipole expansion to identify these amplitudes.

These amplitudes contain valuable information about the interaction of light with the scatterer. In consequence, the decomposition is used in many streams of research. Prime examples would be the study of optical nanoantennas Coenen et al. (2014); Rusak et al. (2014); Stout et al. (2011); Filonov et al. (2012); Pors et al. (2010), the study of meta-atoms and metamaterials Liu and Kivshar (2017); Grahn et al. (2012); Powell (2017); Petschulat et al. (2010) or the analysis of the interaction of scatterers with isolated molecules Vercruysse et al. (2014); Shen et al. (2016); Heaps and Schatz (2017). Using the multipole expansion we can also construct the T-matrix of a scatterer that entirely expresses how an arbitrary incident field is scattered by the pertinent object Waterman (1965); Demésy et al. (). Once the T-matrix of different individual scatterers is known, the interaction of light with a larger cluster of heterogeneous particles can be solved using a multiple-scattering algorithm Doicu et al. (2006); Gimbutas and Greengard (2013).

For spherical particles, the multipole expansions can be calculated analytically using Mie theory Wriedt (2012). However, for more complex structures, the use of numerical solvers is needed in order to obtain the scattered field first and to decompose it afterwards Fruhnert et al. (2017). Among other methods, FEM is specially suitable for the first task, as it can accurately deal with structures possessing complex shapes. When the scattered field () produced by a particle under a certain illumination has been calculated, the multipole coefficients can be obtained thanks to the orthogonality relations of the VSHs


where the integrals are on the surface of a sphere of radius R centered at the origin of coordinates.

The above expressions hold for the radiative VSHs (J=3), which fulfill the radiation condition and can be used for decomposing a scattered field, and for the regular VSHs (J=1), used for example for expressing an illumination field in terms of multipoles.

Therefore, coefficients and in Eq. 1 can be obtained by computing the integrals


However, the fact that the domain of integration has to be a perfect sphere has some drawbacks. First, the numerical solution must be interpolated across the sphere in order to perform the decompositions  5-6. This interpolation produces accuracy losses and it is also computationally expensive, considerably increasing the total computation time.

One possible solution to overcome these drawbacks is to perform the expansion based on volume integrals of the currents induced in the scatterer structures Grahn et al. (2012); Fernandez-Corbaton et al. (2015); Alaee et al. (2018). In this work, we propose a different approach, using the orthogonality property for VSH extended to any closed surface. In this way, any surface which encloses the scatterer can be used to perform the decomposition, including the boundary of the computational domain.

Implementing fast and highly-accurate surface integration methods is normally easier over the natural boundaries of the problem, such as the surfaces of the particles or the surface of the computational domain. Furthermore, most of the available FEM and FDTD software have already efficient built-in implementations for these integrals that we can profit from.

Ii Decomposition of the scattered field over a general surface

We derive in the following an expression to decompose the scattered field into vector spherical harmonics via an integration over a surface with a non-constrained shape. We assume that the scatterer is localized in space and surrounded by a homogeneous, isotropic, and lossless medium.

Let us denote by the field solution of our time-harmonic scattering problem and by another field solution of Maxwell’s equations. Under the above conditions, Maxwell’s equations can be written as


with being the wave number , and are the magnetic permeability and the electric permittivity of the homogeneous medium, respectively, and is the angular frequency of the field.

We consider that both fields fulfil an outwards Silver-Müller radiation condition Kirsch and Hettlich (2016)


We want to obtain a scalar product, which shows an orthogonality relation similar to expressions 2 to 4, but featuring an integral over a surface with a shape not restricted to a sphere. In order to obtain this expression, we start by multiplying Eq. 7 with the complex conjugate of the field


In this derivation we omit for brevity the explicit position dependency of the fields , , and but it is implicitly assumed all the time.

Figure 1: Illustration of the problem considered in this contribution. An object of arbitrary shape generates some scattered field everywhere in the outer region. We consider the integration of Eqn. 10 in a volume surrounding the scattered object. The volume is delimited by two surface boundaries and . Eventually, only the integral across is used to express the multipole coefficients. The shape of this surface can be arbitrary.

Integrating Eq. 10 over a volumetric region that surrounds but does not contain any part of the scatterer (see Fig. 1) we get


Applying two partial integrations over the integration volume


and exploiting the fact that fulfils Maxwell’s equations (Eqn. 8)


we get


where is the boundary of the volume V.

Now we split the surface integral into the contribution of an inner surface and an outer surface (see Fig. 1),


Integrals in the previous equation are defined to have elements pointing inwards and outwards the volume in and , respectively, as seen in Fig. 1.

As we have not imposed any condition on or , the integral values must hold independently of the shape of the surfaces. Since we want an expression that equals relations 2-4, we will consider in the following the exterior surface to be a sphere of radius R ().

Now we will apply the circular shift invariance of the scalar triple product to the terms of . The reason for this change will be clear a few lines below. It reads as


Since , must be independent of the radius of the sphere. We can assume it to be big enough for the far-field approximation to hold. Using spherical coordinates, the Silver-Müller radiation boundary condition given in Eqn. 9 reads


By expressing the differential surface element of a sphere in spherical coordinates,


we can formulate the Silver-Müller radiation condition as


Both fields and fulfil the outward radiation boundary condition Eqn. 9. fulfils the inward radiation boundary condition. We use Eqn. 19 in Eqn. II to get


Note that the last expression holds only in the far-field.

Finally, substituting by or into Eqn. 20, we get from Eqns. II and 20


In comparison with Eqns. 5 and 6, we see that the last integral provides the coefficient of the decomposition of the scattered field into the basis of VSHs times the norm of the VSH in the far-field, i.e.


where we have used that the integral of the norm of every VSH has a value of in the far-field (cf. Ref. Mishchenko et al., 2002, Eqn. C.152). We, therefore, conclude from Eqns. 21 and 22 that we can decompose the scattered field into the vector spherical harmonics as the integral over the boundary


where we have used the well known properties Mishchenko et al. (2002) and to replace by .

We note that it is also possible to reach Eqn. 23 by starting from Lemma 6.38 of Ref. Colton and Kress, 2012 and applying it to and the VSH basis functions.

Expression 23 lets us implement the multipole expansion in an easier way, without the need to perform any expensive interpolation to a spherical domain. It also allows us to use numerical methods in order to improve the numerical stability of the integration as we will explain in next section.

After implementing expression 23 into an FEM solver, we observed large errors in higher order multipoles when decomposing the fields caused by small scatterers. We found that the reason for this numerical error are the singularities of and at the origin of the coordinates. The numerical error produced by the singular fields increases the closer the integration surface is to the origin of the coordinates and the higher the multipole order is.

Furthermore, these singularities make the decomposition very sensitive to numerical noise, amplifying the noise frequently present in the solutions of the scattered fields computed using numerical solvers.

To analyze and to illustrate this problem, we artificially generated a scattered field composed of the multipoles and and computed their analytical expansion. We chose the amplitude of the higher order multipole such that there is a strong near field with an amplitude of around 20 V/m at a distance of away from the origin of coordinates. The actual field considered was V/m. This field was then decomposed using a cubical surface surrounding the origin of coordinates. The implementation of expression 23 was done using a trapezoidal integration with a regular grid of 40.000 points for each of the faces of the cube. We then tried to decompose the scattered field into the multipole contributions , , , . The error of the decomposition is shown in the top graph of Fig. 2 as a function of the length of the cube edge. It can be seen that for the amplitudes absent in the actual source field, the error exponentially increases the closer the integration surface is to the origin of coordinates. Since our final interest is to have an efficient tool for the decomposition of the field obtained from numerical solvers, it is important that the integration surface can be close to the scatterer in order to reduce the computational domain.

Figure 2: Top: Decomposition of the scattered field into the multipoles , , and using expression 23. The integration surface consist of a cube of length centered at the origin of coordinates. Bottom: Same decomposition using expressions 24 and 25.

To mitigate the numerical stability problems, we need to replace the basis functions and by some regular fields, which fulfil an orthogonal relation equivalent to expression 21. The most natural solution are the regularMishchenko et al. (2002) vector spherical harmonics and . These fields also fulfil Maxwell’s equations in homogeneous media but not the Silver-Müller radiation condition, i.e. they are frequently used to expand the incident field. In an appendix we demonstrate the following relations that allow to calculate the amplitudes of the multipole moments as


Note that and are the coefficients of the expansion of into the radiative fields and , even though and are used in the expressions. These equations constitute the final result of this work.

Computing again the same decomposition with the new expression, we got the results shown in the bottom graph of Fig. 2: the problem with the exponentially increasing error is solved with these new expressions. The remaining existing errors are due to the discretization in the trapezoidal integrations.

Iii Finite Element Method implementation

In order to obtain the multipole decomposition of scatterers with complex shapes, we have implemented the derived expressions 24 and 25 into the commercial finite element solver JCMsuite JCMwave GmbH (). The surface integral is computed on the boundary of the computational domain.

In order to test the implementation, we calculate the decomposition of the scattered field generated by a sphere when illuminated with a plane wave. The scattered field produced by a homogeneous sphere is a very well known solution and can be computed using the analytical expressions from Mie theory Wriedt (2012). In the top part of Fig. 3 we show the spectral multipole decomposition of the scattered field produced by the sphere. The results were obtained with the FEM implementation (lines) and with Mie theory (dots). In the bottom part of the figure a convergence test of the FEM solution is shown for an illumination wavelength of 450 nm. The errors were calculated with respect to Mie theory. The sphere was placed into the center of a cubic computational domain. The edges of the cube had a length of and the decomposition was performed on its boundary.

Figure 3: Top. Spectrally resolved multipole expansion of the field scattered by a silver sphere illuminated with a linearly polarized plane wave. The sphere has a radius of 112.5 nm. The plane wave is propagating in direction and the electric field is polarized in direction. The dots correspond to the results from Mie theory. Bottom. Convergence study of the multipolar decomposition with respect to the mesh size of the FEM discretization. The maximum side length of the tetrahedra which compose the mesh is being the refractive index of the different media. The study has been done at a wavelength of 450 nm. The error of the scattering cross section for each multipole order is relative to the value of its cross section.

Finally, we compute the spectral multipole contribution of the scattered field of a silver nanorod illuminated with a linearly polarized plane wave. The nanorod is modelled as a body of revolution with a major and minor semi-axis of 225 nm and 37.5 nm, respectively. The schematic of the simulated structure is shown in the inset of Fig. 4.

For this specific scatterer, the use of expressions 24-25 leads to a reduction in the computational costs of around 70% with regards to Eqns. 5-6. This is because the nanorod is four times longer in the axis than in the and axis. Therefore, the volume of the computational domain needed to surround the nanorod can be approximately four times smaller than if it had to contain a sphere surrounding the scatterer. To be precise, the computational domain consisted of a cuboid with side lengths of 562.5 nm, 300 nm, and 300 nm respectively. The incident plane wave is propagating in the direction, with an amplitude of 1 V/m and is polarized in direction. The results are shown in Fig. 4.

Figure 4: Multipole expansion of the scattered field produced by a silver nanorod. The nanorod is illuminated with a plane wave propagating in -z direction and a wavelength of 450 nm. The incident electric field is linearly polarized in x. The nanorod is modeled as a body of revolution with an elliptical shape where the major axis is parallel to the axis. The major and minor semi-axis have a length of 225 and 37.5 nm, respectively.

We can see how the resonance shifted to longer wavelengths compared with the sphere and how the scattered field is mainly an electric dipolar field at these frequencies. The total scattering cross section at the resonance wavelength is 0.09 . For the sphere, the resonance peak of the total scattering cross section has a value of 0.27 at a wavelength of 364 nm and it has 4 mainly multipole contributions. The geometrical cross sections are 0.027 and 0.039 for the nanorod and the sphere respectively.

Iv Conclusions

We have proposed an orthonormal relation in order to implement the multipole expansion of scattered fields by integrals on generic surfaces. This is particularly suited for use in combination with numerical solvers relying on methods such as FDTD or FEM. The main advantages of this implementation are the convenience and a better tradeoff between accuracy and computational costs, as no interpolations of the scattered field needs to be computed. We have also shown how to use the regular VSH in order to decompose the scattered field into the radiative VSH. The use of regular fields significantly improves the numerical stability of the implementation. Finally, we have tested the implementation using Mie theory and shown one example for a structure which can not be decomposed with analytical means.


Appendix A Replacing radiative VSH by regular VSH

In this appendix, we show that substituting the radiative VSH by the regular VSH as test functions in expression II ends in a relation proportional to Eqn. 21. This new relation lets us express the decompositions 5 and 6 as an integral over a general surface involving the scattered field and the regular VSH. In order to simplify the expressions, in this section we will just use now the VSH . The same procedure can be applied to the fields giving equivalent results.

Let’s start from Eqn. II by substituting with ,


The regular VSH can be decomposed into a combination of radiative and vector spherical harmonics, fulfilling an outwards and inwards radiation condition, respectively,


By plugging Eqn. 27 into Eqn. 26 we get,


By applying now Silver-Müller radiation conditions 28-29, we get


Finally, using Eqn. II


and using the fact that equals , we obtain


An analogous procedure gives an equivalent relation for fields ,


These expressions are just the desired relations corresponding to Eqn. 21 but with the regular VSH. They can be used analogously to derive Eqns. 24 and 25 of the main manuscript.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 675745.
This project has received funding from Deutsche Forschungsgemeinschaft (DFG) grant RO 3640/7-1.
We acknowledge the support from the Karlsruhe School of Optics and Photonics (KSOP).
This work is partially funded through the project 17FUN01 ”BeCOMe” within the programme EMPIR. The EMPIR initiative is co-founded by the European Union’s Horizon 2020 research and innovation programme and the EMPIR Participating Countries.
We also thank Philipp Gutsche for his work on the VSH code and for his contribution to test the implementation.


  1. J. D. Jackson, Classical electrodynamics (AAPT, 1999).
  2. M. I. Mishchenko, L. D. Travis,  and A. A. Lacis, Scattering, absorption, and emission of light by small particles (Cambridge university press, 2002).
  3. T. Coenen, F. B. Arango, A. F. Koenderink,  and A. Polman, Nat. Commun. 5, 3250 (2014).
  4. E. Rusak, I. Staude, M. Decker, J. Sautter, A. E. Miroshnichenko, D. A. Powell, D. N. Neshev,  and Y. S. Kivshar, Appl. Phys. Lett. 105, 221109 (2014).
  5. B. Stout, A. Devilez, B. Rolly,  and N. Bonod, J. Opt. Soc. Am. B 28, 1213 (2011).
  6. D. S. Filonov, A. E. Krasnok, A. P. Slobozhanyuk, P. V. Kapitanova, E. A. Nenasheva, Y. S. Kivshar,  and P. A. Belov, Appl. Phys. Lett. 100, 201113 (2012).
  7. A. Pors, M. Willatzen, O. Albrektsen,  and S. I. Bozhevolnyi, J. Opt. Soc. Am. B 27, 1680 (2010).
  8. W. Liu and Y. S. Kivshar, Philos. Trans. Royal Soc. A 375, 20160317 (2017).
  9. P. Grahn, A. Shevchenko,  and M. Kaivola, New J. Phys. 14, 093033 (2012).
  10. D. A. Powell, Phys. Rev. Appl. 7, 034006 (2017).
  11. J. Petschulat, J. Yang, C. Menzel, C. Rockstuhl, A. Chipouline, P. Lalanne, A. Tünnermann, F. Lederer,  and T. Pertsch, Opt. Express 18, 14454 (2010).
  12. D. Vercruysse, X. Zheng, Y. Sonnefraud, N. Verellen, G. Di Martino, L. Lagae, G. A. Vandenbosch, V. V. Moshchalkov, S. A. Maier,  and P. Van Dorpe, ACS Nano 8, 8232 (2014).
  13. H. Shen, R. Y. Chou, Y. Y. Hui, Y. He, Y. Cheng, H. C. L. Tong, Q. Gong,  and G. Lu, Laser Photonics Rev. 10, 647 (2016).
  14. C. W. Heaps and G. C. Schatz, J. Chem. Phys 146, 224201 (2017).
  15. P. C. Waterman, Proc. IEEE 53, 805 (1965).
  16. G. Demésy, B. Stout,  and J.-C. Auger,  arXiv:1802.00596 .
  17. A. Doicu, T. Wriedt,  and Y. A. Eremin, Light scattering by systems of particles: null-field method with discrete sources: theory and programs (Springer, 2006).
  18. Z. Gimbutas and L. Greengard, J. Comput. Phys. 232, 22 (2013).
  19. T. Wriedt, Mie theory: a review (Springer, 2012) pp. 53–71.
  20. M. Fruhnert, I. Fernandez-Corbaton, V. Yannopapas,  and C. Rockstuhl, Beilstein J. Nanotechnol. 8, 614 (2017).
  21. I. Fernandez-Corbaton, S. Nanz, R. Alaee,  and C. Rockstuhl, Opt. Express 23, 33044 (2015).
  22. R. Alaee, C. Rockstuhl,  and I. Fernandez-Corbaton, Opt. Commun. 407, 17 (2018).
  23. A. Kirsch and F. Hettlich, Mathematical Theory of Time-harmonic Maxwell’s Equations (Springer, 2016).
  24. D. Colton and R. Kress, Inverse acoustic and electromagnetic scattering theory (Springer Science & Business Media, 2012).
  25. JCMwave GmbH,
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Comments 0
The feedback must be of minumum 40 characters
Add comment

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question