Volume Integral Formulation for the Calculation of Material Independent Modes of Dielectric Scatterers

TL;DR

This paper introduces a volume integral eigenvalue approach to compute material-independent electromagnetic modes of dielectric scatterers, enabling efficient analysis of scattering and resonances across various permittivities.

Contribution

It presents a novel modal basis derived from an auxiliary eigenvalue problem, simplifying the calculation of scattered fields for different permittivities.

Findings

Moderate number of modes suffices for accurate far-field scattering description.

Method effectively analyzes plasmonic and photonic resonances.

Computational efficiency improves over direct volume integral solutions for multiple permittivities.

Abstract

In the frame of volume integral equation methods, we introduce an alternative representation of the electromagnetic field scattered by a homogeneous object of arbitrary shape at a given frequency, in terms of a set of modes independent of its permittivity. This is accomplished by introducing an auxiliary eigenvalue problem, based on a volume integral operator. With this modal basis the expansion coefficients of the scattered field are simple rational functions of the permittivity of the scatterer. We show, by studying the electromagnetic scattering from a sphere and a cylinder of dimensions comparable to the incident wavelength, that only a moderate number of modes is needed to accurately describe the scattered far field. This method can be used to investigate resonant scattering phenomena, including plasmonic and photonic resonances, and to design the permittivity of the object to pursue a prescribed tailoring of the scattered field. Moreover, the presented modal expansion is computationally advantageous compared to direct solution of the volume integral equation when the scattered field has to be computed for many different values of the dielectric permittivity, given the size and shape of the dielectric body.

Tables

Table 1. TABLE I: Values of γ h = σ h − 1 + 1 subscript 𝛾 ℎ superscript subscript 𝜎 ℎ 1 1 \gamma_{h}=\sigma_{h}^{-1}+1 associated to the eigenvalues σ h subscript 𝜎 ℎ \sigma_{h} of the ten modes shown in Fig. 7 .

#	1	2	3	4	5
${\gamma }_h$	$0.45-0.74i$	$0.88-1.2i$	$1.6-1.3i$	$3.6-0.88i$	$3.7-0.49i$
#	6	7	8	9	10
${\gamma }_h$	$4.5-1.4i$	$5.6-0.58i$	$6.0-0.43i$	$7.2-0.24i$	$9.6-1.4i$

Equations

$$\mathbf{J}\left(\mathbf{r}\right)=\{\begin{array}{ll}-i\omega {\epsilon }_0\chi \left(\omega \right)\left({\mathbf{E}}_S^{+}\left(\mathbf{r}\right)+{\mathbf{E}}_i\left(\mathbf{r}\right)\right)& \mathbf{r}\in \mathring{\Omega }\\ 0& \mathbf{r}\in {\mathbb{R}}^3\backslash \Omega \end{array}$$

$$\begin{gathered} {\mathbf{E}}_S^{+}\left(\mathbf{r}\right)=+i\omega {\mu }_0{\iiint }_{\Omega }\mathbf{J}\left({\mathbf{r}}^{\prime }\right)g\left({\mathbf{r}-\mathbf{r}}^{\prime }\right)\text{d}{V}^{\prime } \\ +\frac{1}{i\omega {\epsilon }_0}\nabla {\text{varoiint}}_{\partial \Omega }\mathbf{J}\left({\mathbf{r}}^{\prime }\right)\cdot \hat{\mathbf{n}}\left({\mathbf{r}}^{\prime }\right)g\left({\mathbf{r}-\mathbf{r}}^{\prime }\right)\text{d}{S}^{\prime }\forall \mathbf{r}\in \Omega , \end{gathered}$$

$${\chi }^{-1}\mathbf{J}\left(\mathbf{r}\right)-\mathcal{L}\left\{\mathbf{J}\right\}\left(\mathbf{r}\right)=-i\omega {\epsilon }_0{\mathbf{E}}_i\left(\mathbf{r}\right)\forall \mathbf{r}\in \Omega ,$$

$$\begin{gathered} \mathcal{L}\left\{\mathbf{J}\right\}\left(\mathbf{r}\right)=-\nabla {\text{varoiint}}_{\partial \Omega }\mathbf{J}\cdot \hat{\mathbf{n}}\left({\mathbf{r}}^{\prime }\right)g\left({\mathbf{r}-\mathbf{r}}^{\prime }\right)\text{d}{S}^{\prime } \\ +k_0^2{\iiint }_{\Omega }\mathbf{J}\left({\mathbf{r}}^{\prime }\right)g\left({\mathbf{r}-\mathbf{r}}^{\prime }\right)\text{d}{V}^{\prime }\forall \mathbf{r}\in \Omega . \end{gathered}$$

$$⟨{\mathbf{J}}^{\prime },\mathcal{L}{\mathbf{J}}^{\prime \prime }{⟩}_{\Omega }=⟨\mathcal{L}{\mathbf{J}}^{\prime },{\mathbf{J}}^{\prime \prime }{⟩}_{\Omega }$$

$$⟨\mathbf{A},\mathbf{B}{⟩}_V={\iiint }_V\mathbf{A}\cdot \mathbf{B}\text{mbox}dV.$$

$$\mathcal{L}\left\{\mathbf{u}\right\}\left(\mathbf{r}\right)=\sigma \mathbf{u}\left(\mathbf{r}\right)$$

$${\mathcal{L}}^{†}={\mathcal{L}}^{\ast }$$

$$⟨{\mathbf{u}}_r,{\mathbf{u}}_s{⟩}_{\Omega }=0\text{for}{\sigma }_r\ne {\sigma }_s.$$

$$\text{mbox}Re\left\{{\sigma }_r^{-1}\right\}$$

$$\text{mbox}Im\left\{{\sigma }_r^{-1}\right\}$$

$$\mathbf{J}={\mathbf{J}}_L+{\mathbf{J}}_S,$$

$$\nabla \cdot {\mathbf{J}}_L$$

$$\nabla \cdot {\mathbf{J}}_S$$

$${\mathbf{w}}_k^L\left(\mathbf{r}\right)=\nabla \times {\mathbf{N}}_k\left(\mathbf{r}\right)$$

$$\mathbf{J}=\sum _{k=1}^{N_L}{I}_k^L{\mathbf{w}}_k^L+\sum _{k=1}^{N_S}{I}_k^S{\mathbf{w}}_k^S,$$

$${\chi }^{-1}\mathbf{R}\mathbf{I}-\mathbf{K}\mathbf{I}=\mathbf{V}$$

$$\mathbf{K}=i\omega \mathbf{L}+\frac{1}{i\omega }\mathbf{D},$$

$$\mathbf{R}=\left[\begin{array}[c]{lll}\mathbf{R}_{LL}&\mathbf{R}_{LS}\\ \mathbf{R}_{SL}&\mathbf{R}_{SS}\end{array}\right],\quad\mathbf{L}=\left[\begin{array}[c]{lll}\mathbf{L}_{LL}&\mathbf{L}_{LS}\\ \mathbf{L}_{SL}&\mathbf{L}_{SS}\\ \end{array}\right],\\ \mathbf{D}=\left[\begin{array}[c]{ll}\mathbf{0}&\mathbf{0}\\ \mathbf{0}&\mathbf{D}_{SS}\end{array}\right],\quad\mathbf{I}=\left[\begin{array}[c]{c}\mathbf{I}^{L}\\ \mathbf{I}^{S}\end{array}\right],\quad\mathbf{V}=\left[\begin{array}[c]{c}\mathbf{V}^{L}\\ \mathbf{V}^{S}\end{array}\right].$$

$${\left({\mathbf{R}}_{\alpha \beta }\right)}_{pq}$$

$${\left({\mathbf{L}}_{\alpha \beta }\right)}_{pq}$$

$${\left({\mathbf{D}}_{SS}\right)}_{pq}=\frac{1}{{\epsilon }_0}\times$$

$${\text{varoiint}}_{\partial \Omega }{\text{varoiint}}_{\partial \Omega }\left({\mathbf{w}}_q^S\cdot {\hat{\mathbf{n}}}^{\prime }\right)\left({\mathbf{r}}^{\prime }\right)\left({\mathbf{w}}_p^S\cdot \hat{\mathbf{n}}\right)\left(\mathbf{r}\right)g\left({\mathbf{r}-\mathbf{r}}^{\prime }\right)\text{d}{S}^{\prime }\text{d}S.$$

$${\left({\mathbf{V}}^{\alpha }\right)}_p={\iiint }_{\Omega }{\mathbf{w}}_p^{\alpha }\cdot {\mathbf{E}}_i\text{d}V,\forall \alpha \in \left\{L,S\right\}.$$

$$\mathbf{K}\mathbf{I}=\sigma \mathbf{R}\mathbf{I},$$

$${\mathbf{I}}_h^T\mathbf{R}{\mathbf{I}}_k=0{\sigma }_h\ne {\sigma }_k.$$

$$\mathbf{I}=\sum _{h=1}^{N_{DoF}}\frac{1}{{\chi }^{-1}-{\sigma }_h}\frac{{\mathbf{I}}_h^T\mathbf{V}}{{\mathbf{I}}_h^T\mathbf{R}{\mathbf{I}}_h}{\mathbf{I}}_h$$

$$\left|\text{mbox}Re\left\{{\sigma }_h\right\}\right|>\xi \cdot \underset{k=1\dots N_{tot}}{\max }\left|\text{mbox}Re\left\{{\sigma }_k\right\}\right|$$

$$\left|\text{mbox}Im\left\{{\sigma }_h\right\}\right|>\xi \cdot \underset{k=1\dots N_{tot}}{\max }\left|\text{mbox}Im\left\{{\sigma }_k\right\}\right|$$

$${\mathbf{E}}_i=e^{i{k}_{0}z}\hat{\mathbf{x}}$$

$$C_{sca}=\frac{c_0}{\omega }{\text{varoiint}}_{S_c}{\hat{\mathbf{e}}}_r\cdot \text{mbox}Im\left\{{\left(\nabla \times {\mathbf{E}}_S^{-}\right)}^{\ast }\times {\mathbf{E}}_S^{-}\right\}\text{mbox}dS,$$

$${\mathbf{E}}_S^{\infty }\left(\theta ,\varphi \right)=\underset{r\to \infty }{\lim }\left[r{e}^{-i{k}_{0}r}{\mathbf{E}}_S^{-}\left(r,\theta ,\varphi \right)\right].$$

$${\mathbf{E}}_S^{-}=\sum _{n=1}^{n_{max}}{E}_n\left(i{a}_n{\mathbf{N}}_{e1n}^{\left(3\right)}\left(k_0\mathbf{r}\right)-b_n{\mathbf{M}}_{o1n}^{\left(3\right)}\left(k_0\mathbf{r}\right)\right)$$

Full text

This method can be used to investigate resonant scattering phenomena, including plasmonic and photonic resonances, and to design the permittivity of the object to pursue a prescribed tailoring of the scattered field. Moreover, the presented modal expansion is computationally advantageous compared to direct solution of the volume integral equation when the scattered field has to be computed for many different values of the dielectric permittivity, given the size and shape of the dielectric body.

I Introduction

The solution of the scattering problem obtained by the majority of analytical and numerical method is often expressed in a form where the contributions of the material and of the geometry are interwoven and cannot be separated. For instance, in the Mie theory for the electromagnetic scattering from spherical homogeneous objects, the electric size and the permittivity both appear in the argument of the vector spherical wave functions, and, as a consequence, the Mie expansion coefficients are complicated function of their combination. Similar considerations hold true also when the scattering problem is solved by using Volume Integral Equations (VIE) [1, 2, 3, 4, 5], Surface Integral Equations (SIE) [6] or the Finite Difference Time Domain (FDTD) [7] solvers. Thus, the complex dependence of the solution of the scattering problem on the permittivity makes the design of the material to achieve assigned constraints on the scattered field complicated.

Over the years some authors have investigated the possibility of separating the dependence on the material properties from the dependence on the geometry by representing the scattered field in terms of a set of modes independent of the permittivity of the scatter. In the following we denote these modes as Material Independent Modes (MIMs). MIMs have been calculated in the quasi-static limit [8, 9, 10, 11, 12, 13, 14, 4], for the scalar Mie scattering [15]. They have been also derived within the quasi-static [16] and retarded [17] single dipole approximation.

More recently, MIMs have been derived for the full-retarded vector scattering by a homogeneous sphere [18], by a coated sphere [19], and by a flat slab [20]. Furthermore, it has been also shown that the investigation of the eigenvalues associated to the MIMs unveils important structural properties of plasmonic and photonic resonances [21]. The expansion of the scattered field in terms of MIMs leads to a natural separation of the contributions of the permittivity and of the geometry because the corresponding expansion coefficients are simple rational functions of the permittivity, whereas the MIMs depend solely on the geometry. This fact has suggested a straightforward methodology to design the permittivity of the object to pursue a prescribed tailoring of the scattered field, including the cancellation of the backscattering, the suppression of a given multipolar order, and the maximization of the scattered field in the near-field zone [18, 19].

MIMs could also be useful in solving imaging problems. In particular, it has been already shown that in Eddy Current Tomography the eigenvalues related to the modes are monotonic with respect to the electrical resistivity [22, 23]. This property has enabled the use of a real-time imaging algorithm [24].

Deriving a general and accurate method to compute MIMs of arbitrarily shaped particle is a problem of great importance that has been addressed so far only in the electrostatic regime [12, 14, 4]. In this paper, we tackle this problem by introducing a general approach for the calculation of MIMs in arbitrary shaped homogeneous particles based on a volume integral formulation of the Maxwell’s equation. This approach reduces to the one presented in Ref. [4] in the electrostatic limit. The problem unknown is the sum of the conduction and polarization current densities. It is represented in terms of its loop and star components.

The MIMs are not the only modes that can be defined starting from the full-wave electromagnetic scattering problem in the frequency domain. Different choices are possible resulting in different kinds of modes. Here, we briefly outline two alternative mode definitions, the quasi normal modes and the characteristic modes, pointing out their profound differences with the MIMs.

The Quasi Normal Modes (QNMs), also known as morphology dependent modes and leaky modes [25, 26] are often used to investigate open systems. Reviews of QNMs can be found in different contexts, including light scattering [27] and nanophotonics [28]. The QNMs are not orthogonal in the usual sense, moreover they diverge exponentially at large distances [28]. Therefore, to be used in any practical application they need to be normalized [29]. QNMs of a homogeneous dielectric object do depend on the material, shape and size but are independent of the frequency.

Characteristic modes (CMs), also known as characteristic currents, were first introduced for perfect conductors by the pioneering works of R. J. Garbacz [30, 31], R. F. Harrington and J. Mautz [32]. Later, Harrington et al. have proposed generalization of CMs to structures with dielectric and magnetic materials [33, 34]. CMs have been extensively used in the analysis and design of radiation and scattering of two and three dimensional objects [35], providing guidance for tailoring the excitation of antennas or scatterers to selectively excite desirable radiation modes. CMs have been also used as a set of basis functions for the solution of the scattering problem from an array of small objects [36], resulting in significant reduction of the total number of unknowns compared to a standard MoM approach. CMs are real and satisfy a weighted orthogonality. Although CMs are independent of the excitation, they do depend on the shape, size, frequency, and material composition of the scattering object.

The layout of the paper is as follows. In Sec. II we derive an auxiliary eigenvalue problem which defines the MIMs starting from a volume integral formulation of the Maxwell’s equations. We also discuss the main properties of the eigenvalues and of the associated modes. Then, in Sec. III we introduce the numerical discretization of the unknowns in terms of the star and loop shape functions. We provide the explicit expression of the matrix equations. Thus, we use the MIMs to expand the scattered field in the presence of an arbitrary external excitation. In Sec. IV we apply the introduced numerical approach to compute the MIMs of a sphere and a cylinder of dimensions comparable to the incident wavelength. In particular, the eigenvalues associated to MIMs of a homogeneous sphere are validated against analogous quantities evaluated through an independent approach [18]. The scattering cross section and the radiation diagrams of both sphere and cylinder are obtained starting from MIMs, and compared against the Mie theory and a surface integral formulation of the Maxwell’s equations. It is also shown by numerical examples that only few modes are necessary to represent the far field scattered by a dimension comparable to the wavelength. A scenario, in which the presented method is computational advantageous compared to the direct solution of the volume integral equations is also illustrated.

II Mathematical Model

Let us consider an isotropic and homogeneous material occupying a volume $\Omega$, which is bounded by a closed surface $\partial\Omega$ with outward-pointing normal $\hat{\bf n}$. The medium is non-magnetic with relative permittivity $\varepsilon_{R}\left(\omega\right)=\varepsilon_{R}^{\prime}\left(\omega\right)+i\,\varepsilon_{R}^{\prime\prime}\left(\omega\right)$, surrounded by vacuum. The object is excited by a time harmonic electromagnetic field incoming from infinity $\mbox{Re}\left\{{\bf E}_{i}\left({\bf r}\right)e^{-i\omega t}\right\}$. Let ${\bf E}_{S}^{+}$ and ${\bf E}_{S}^{-}$ be the scattered electric fields in $\mathring{\Omega}$ and $\mathbb{R}^{3}\backslash\bar{\Omega}$, respectively, where $\mathring{\Omega}$ denotes the interior of $\Omega$.

Following Ref. [37], we now consider the generalized current density in the material region:

[TABLE]

where $\chi\left(\omega\right)=\varepsilon_{R}\left(\omega\right)-1$ is the dielectric susceptibility, and $\varepsilon_{0}$ is the free space dielectric constant. We choice $\mathbf{J}$ as the problem’s unknown. This fact allows us to limit the discretization only to the material regions and to recast the problem as a volume integral equation whose kernel is the free-space Green function [38, 39, 37].

We express the electric field produced by $\bf J$ in terms of the vector and scalar potentials, assuming the Lorenz gauge. The vector and scalar potentials can be directly obtained from the current density and the surface charge density through the scalar free space Green function $g\left(\mathbf{r}\right)\triangleq e^{ik_{0}r}/\left(4\pi r\right)$, $k_{0}=\omega/c$, and $c$ is the speed of light in vacuum. Therefore, the electric field ${\bf E}_{S}^{+}$ sustained by $\mathbf{J}$ is:

[TABLE]

where the quantity ${\bf J}$ occurring in the surface integral on the r.h.s. represents the limit of the volume current $\bf J$ along the normal to $\partial\Omega$ as the evaluation point approaches the surface from the internal face of $\partial\Omega$. By combining Eqs. 1 and 2, we obtain the volume integral equation in terms of the unknown $\mathbf{J}$:

[TABLE]

where we have introduced the operator $\mathcal{L}\left\{\cdot\right\}$

[TABLE]

Since the material is assumed to be homogeneous, ${\bf J}$ is divergence free in $\Omega$, whereas its normal component on the domain boundary is proportional to surface charge density. It is straightforward to show that the operator $\mathcal{L}$ is symmetric:

[TABLE]

where

[TABLE]

Aiming at the reduction of the scattering problem to an algebraic form, we introduce the following auxiliary eigenvalue problem

[TABLE]

where $\sigma$ is the eigenvalue. The adjoint of the operator of $\mathcal{L}$ is

[TABLE]

Therefore, the operator $\mathcal{L}$ is not self-adjoint but symmetric. Its eigenvalues are complex with $\mbox{Im}\left\{\sigma_{r}^{-1}\right\}<0$. The eigenmodes ${\bf u}_{r}$ and ${\bf u}_{s}$ corresponding to different eigenvalues $\sigma_{r}$ and $\sigma_{s}$ are not orthogonal in the usual sense, i.e. $\langle{\bf u}_{r}^{*},{\bf u}_{s}\rangle_{\Omega}\neq 0$ Nevertheless, it can be proved that

[TABLE]

We denote with $\mathbf{C}_{r}$ the electric field produced in the whole space by the current ${\bf u}_{r}$. By considering the volume integral of the quantity $\nabla\cdot\left[\mathbf{C}_{r}\times\nabla\times\mathbf{C}_{r}^{*}\right]$ over $\mathbb{R}^{3}\backslash\bar{\Omega}$, exploiting the divergence theorem and the properties of the eigenfunctions we readily obtain

[TABLE]

where $\left\|{\bf A}\right\|_{V}^{2}=\langle{\bf A}^{*},{\bf A}\rangle_{V}$ and $S_{\infty}$ is an auxiliary closed surface enclosing the scatterer and contained in the far zone. Equation 10 suggests that $\mbox{Re}\left\{{\sigma_{r}^{-1}}\right\}$ does not have a definite sign, while Equation 11 implies that $\mbox{Im}\left\{\sigma_{r}^{-1}\right\}$ is strictly negative. In particular, $\mbox{Im}\left\{\sigma_{r}^{-1}\right\}$ is proportional to the contribution of the corresponding mode to the power radiated to infinity, accounting for its radiative losses.

III Numerical Formulation

Following [37], we split the unknown $\mathbf{J}$ into the sum of its loop and star components, denoted by $\mathbf{J}_{L}$ and $\mathbf{J}_{S}$, respectively, namely:

[TABLE]

where:

[TABLE]

Thus, we introduce a finite-dimensional approximation of the currents $\mathbf{J}_{L}$, $\mathbf{J}_{S}$ in terms of linear combinations of suitable shape functions, namely $\mathbf{w}_{k}^{L}$’s and $\mathbf{w}_{k}^{S}$’s, respectively. The $k$-th loop shape function $\mathbf{w}_{k}^{L}$ is associated to the $k$-th edge of the finite element discretization of the volume $\Omega$, and it is defined as the curl of the $k$-th edge-element shape functions $\mathbf{N}_{k}$:

[TABLE]

The shape functions $\mathbf{w}_{k}^{S}$’s are used to discretize the star component $\mathbf{J}_{S}$ describing the effects due to the surface charges appearing on the boundary ${\partial\Omega}$. The function $\mathbf{w}_{k}^{S}$ is defined as the curl of the $k$-th edge-element shape function corresponding to an edge on ${\partial\Omega}$.

Then, the unknown current density distribution is represented, at the discrete level, as

[TABLE]

where the coefficients $I_{k}^{L}$’s and $I_{k}^{S}$’s are the degrees of freedom (DoFs) for the loop and star components, while $N_{L},\ N_{S}$ are the number of loop and star functions involved in the discretization. Their sum correspond to the total number of degrees of freedom $N_{DoF}=N_{L}+N_{S}$.

We obtain the discrete model, by first substituting the current representation of Eq. 13 into the integral equation 3, and by then applying the Galerkin method, projecting along the loop and the star shape functions. Following the outlined steps, we obtain the discretization of the non-homogeneous problem:

[TABLE]

where

[TABLE]

and the block matrices and vector are defined as:

[TABLE]

It is worth to point out that the use of the Galerkin method guarantees that the matrix $\mathbf{K}$, which represents the discretization of the operator $\mathcal{L}$, preserves its symmetry. The generic occurrences of the block matrices $\bf R$ and $\bf L$ are:

[TABLE]

$\forall\alpha,\beta\in\left\{L,S\right\}$. The non-vanishing elements of the matrix $\mathbf{D}$ are:

[TABLE]

The occurrences of the vector $\mathbf{V}$ are:

[TABLE]

In absence of external excitation from Eq. 14 we obtain the discrete eigenvalue problem

[TABLE]

where $\sigma$ is the eigenvalue, We denote the numerical eigenvalues with $\sigma_{h}$ and the corresponding eigenvectors with ${\bf I}_{h}$. The eigenvectors are not orthogonal in the usual sense but we have:

[TABLE]

By expanding the solution ${\bf I}$ of the discrete problem of Eq. 14 in terms of the eigenvalues $\sigma_{h}$ and the eigenvectors ${\bf I}_{h}$ of Eq. 16 we obtain:

[TABLE]

In Eq. 18 the dependence of ${\bf I}$ from the material and the geometry are naturally separated. The eigenvalues $\sigma_{h}$ and the eigenvectors ${\bf I}_{h}$ are permittivity independent, and they depend on the shape and size of the dielectric object, and on the frequency. The susceptibility appears in the multiplicative factors only as $1/\left(\chi^{-1}-\sigma_{h}\right)$.

For passive materials with non-negative imaginary part of the susceptibility we have $\mbox{Im}\left\{\chi\right\}\geq 0$, thus the quantity $\left|\chi^{-1}-\sigma_{h}\right|$ in Eq. 18 does not vanish because $\mbox{Im}\left\{\sigma_{h}^{-1}\right\}<0$. Nevertheless, the amplitude of the h-th mode increases as the distance between $\chi^{-1}$ and $\sigma_{h}$ is reduced. In other words, Eq. 18 exemplifies that, for a fixed frequency, when the scatterer’s material closely “matches” an eigenvalue, the corresponding mode undergoes a boosting, namely a “resonance” in a “material picture”. This picture is dual with respect to the usual “frequency picture”, where the material is instead fixed and the frequency plays the role of spectral parameter. The “material picture” is particularly relevant in light of the latest advances in Metamaterials’ design and fabrication techniques, which are enabling the effective value of material’s permittivity and permeability to be engineered with increasing precision.

Moreover, when $\mbox{Re}\left\{\sigma_{h}^{-1}\right\}<-1$ the amplitude of the h-th mode will be particularly strong for materials with negative permittivity, in this case we denote the corresponding mode as a “plasmonic” mode. On the contrary when $\mbox{Re}\left\{\sigma_{h}^{-1}\right\}>-1$ will be particularly strong for materials with positive permittivity, in this case we denote the corresponding mode as “photonic”.

In the next sections, we refer to the introduced method as MIM-VIE.

IV Numerical Results

We carried out the computation presented in this section on a cluster of 25 Intel Xeon CPU E5-2690 cores operating at 2.90GHz and equipped with 128 GB of RAM memory. The Fortran library ARPACK [40] has been used for the calculation of the generalized eigenvalues and eigenvectors of Eq. 16. These routines are based upon an algorithmic variant of the Arnoldi process called the Implicitly Restarted Arnoldi Method (IRAM) [40]. In particular, we compute $N_{SR}$ eigenvalues of smallest real part, $N_{LI}$ eigenvalues of largest imaginary part, $N_{SI}$ eigenvalues of smallest imaginary part, and we consider the union of these three sets, whose dimension is denoted with $N_{tot}$. We set $N_{SR}=N_{LI}=N_{SI}=200$ and a converge tolerance $\tau=10^{-2}$. It is apparent that $N_{tot}\leq N_{SR}+N_{LI}+N_{SI}$, because these three sets may have finite intersection. It is worth to point out that the set of eigenvalues with the largest real part have not been computed since they are associated to modes of very high order playing no role in the scattering process.

Subsequently, we filter the modes by retaining only those associated to the eigenvalues located outside a box centred in the origin of the complex plane, satisfying simultaneously the following criteria

[TABLE]

where we assumed different values of $\xi$. We denote with $h_{max}$ the total number of modes which pass this filter. In this section, we restrict the sum of Eq. 18 to only these $h_{max}$ MIMs.

In this section, we investigate the scattering cross section and the radiation patterns of a sphere and a cylinder by using the MIMs expansion. They are excited by a plane wave of unit intensity polarized along $\hat{x}$ and propagating along $\mathbf{\hat{z}}$,

[TABLE]

In particular, the scattering cross section is [41, 42]:

[TABLE]

where $S_{c}$ is an auxiliary surface enclosing $\Omega$. It is worth noting that interference among several MIMs may take place in the total scattered power, because the MIMs are not orthogonal[21]. The radiation pattern is [42]

[TABLE]

IV-A Sphere

Let us consider a homogeneous sphere of diameter $D=\lambda$, where $\lambda$ is the wavelength in vacuum. We considered an hexahedral volume mesh with $N_{p}=26075$ points and $N_{e}=24800$ elements, leading to $N_{L}=48401$ loop shape functions, and $N_{S}=2399$ star shape functions, corresponding to a total of $50800$ complex unknowns.

First, we validate the calculation of the eigenvalues $\sigma_{h}$ against the corresponding quantities calculated by using the approach introduced in Ref. [18]. The latter approach consists in numerically finding the the roots of a polynomial, whose coefficients are analytically known for a sphere. For simplicity of representation, we compare the quantities $\gamma_{h}=\sigma_{h}^{-1}+1$. In Fig. 1, we show with a black circle the $\gamma_{h}$’s associated to the eigenvalues calculated with the approach introduced in this paper, with a red down-pointing triangle the same quantity associated to magnetic (TE) modes, and with a blue up-pointing triangle the eigenvalues associated to electric (TM) modes, both calculated by using Ref. [18]. Only the $\gamma_{h}$’s belonging to the box of the complex plane $\left[-10,10\right]\times\left[-4,0\right]$ are shown. They are obtained by using the filter of Eq. 19 with $\xi=10^{-3}$ are shown. A small error is appreciable only for moderately positive $\text{Re}\left\{\gamma_{h}\right\}$, where the MIM-VIE method tends to overestimate the real part of $\gamma_{h}$. This error does not improve by increasing the tolerance $\tau$ up to $10^{16}$, while the eigenvalues’ positions remain the same. It is worth noting that the quantities shown in Fig. 1 depend neither on the excitation conditions nor on the permittivity $\varepsilon_{r}$, but they solely depend on the quantity $D/\lambda$.

Next, in Fig. 2 we show the scattering cross section $C_{sca}$ of the sphere as a function of the real part of its permittivity $\varepsilon_{r}^{\prime}\in\left[-10,10\right]$, whereas the imaginary part is fixed to the value $\varepsilon_{r}^{\prime\prime}=0.1$. The sphere is excited by the plane wave of Eq. 20 with $\lambda=1m$. We computed $C_{sca}$ using Eqs. 21 and 18. As a reference solution we use the standard Mie Theory [43], where the scattered electric field has been expanded in terms of vector spherical wave functions (VSWFs):

[TABLE]

where $a_{n}$ and $b_{n}$ are the Mie scattering coefficients, which can be found in Ref. [43] and $E_{n}=i^{n}(2n+1)/\left[n(n+1)\right]$ and the cross section is obtained by

[TABLE]

where we assumed $n_{max}=10$. The functions ${\bf N}_{emn}^{\left(3\right)}$ and ${\bf M}_{omn}^{\left(3\right)}$ are the VSWFs of the radiative kind and the subscripts $e$ and $o$ denote even and odd azimuthal dependence.

Specifically, in Fig. 2 (a) we investigate the convergence of $C_{sca}$ as we increase the number of considered MIMs. In particular, we performed the filtering of the modes by using Eq. 19, and assuming $\xi=5\cdot 10^{-2}$, $\xi=5\cdot 10^{-3}$, and $\xi=10^{-3}$, thus retaining $h_{max}=48$ (green line), $h_{max}=89$ (blue line), and $h_{max}=147$ (black line) eigenvalues, respectively. We assumed a tolerance of $\tau=10^{-2}$. We compare the solutions with the standard Mie theory by using $n_{max}=10$ (red line). It is apparent that $48$ modes are able to satisfactory describe the $C_{sca}$ only within a small interval centred at $\varepsilon_{r}^{\prime}=1$. Then, by considering $h_{max}=89$ MIMs we note an increase of the permittivities’interval in which the VIE solution agrees well with the MIE solution. Eventually, for $h_{max}=147$ we achieved a very good agreement within the whole investigated interval. We also appreciate a slight shift of the peaks of the MIM-VIE solution with respect to the reference Mie solution for moderate positive permittivities, due to the overestimation of the real part of $\gamma_{h}$, already apparent in Fig. 1.

In Fig. 2 (b) we investigate the accuracy of the computed $C_{sca}$ as a function of the tolerance $\tau$. In particular, we considered $\tau=10^{-2}$, $\tau=10^{-3}$, and $\tau=10^{-16}$. The number of modes was instead fixed to $h_{max}=147$ ($\xi=10^{-3}$). We compare the solutions with the standard Mie theory by using $n_{max}=10$ (red line). It is apparent that with tolerance of $\tau=10^{-2}$ we achieve a good agreement with the Mie theory. A further increase of the tolerance does not appreciably improves the solution.

Next, we show in Fig. 3 the radiation pattern intensity $\left|\mathbf{E}_{S}^{\infty}\left(\theta,\phi\right)\right|^{2}$ as a function of the inclination angle $\theta$ in the $zx$ plane (a), (c) and in the $yz$ plane (b), (d). In particular, we investigate the convergence of $\left|\mathbf{E}_{S}^{\infty}\left(\theta,\phi\right)\right|^{2}$ as a function of the number of considered MIMs. The sphere is excited by the plane wave of Eq. 20 with $\lambda=1$. The angles $\theta=0$ and $\theta=\pi$ correspond to the directions of forward- and back- scattering, respectively. We assumed a permittivity $\varepsilon_{r}=5+0.1i$, and we performed the calculation using Eq. 18 with $h_{max}=48$ (green line), $h_{max}=89$ (blue line), and $h_{max}=147$ (black line) MIMs and by the standard Mie theory (red line). We obtained good agreement only for $h_{max}=147$.

In conclusion, the introduced approach correctly evaluates the eigenvalues $\sigma_{h}$ of a sphere of diameter $D=\lambda$. Moreover, it is also capable to accurately compute the scattering cross section and the radiation diagrams by using only a moderate number of MIMs.

We now show that the MIM-VIE solution may be computational advantageous in a scenario relevant to practical applications if compared to a direct solution of the VIE formulation of Eq. 14, namely

[TABLE]

First, let us derive the time needed for the calculation of the scattering response at a given frequency of a number $n_{\varepsilon_{r}}$ of different spheres of given size and shape, but having $n_{\varepsilon_{r}}$ different values of $\varepsilon_{r}$, by using both the MIM-VIE and the direct VIE solvers. The cumulative time for assembling both matrices $\bf K$ and $\bf R$ on $25$ processors is the same for both approaches, namely $t_{\text{ASM}}=2h\,12m$. The total time needed for the calculation of the sets of eigenvalues and corresponding MIMs is $t_{\text{EIG}}=44m$. The total time needed for the MIM-VIE formulation is with good approximation independent of $n_{\varepsilon_{r}}$, being approximatively $t_{\text{MIM-VIE}}\approx t_{\text{ASM}}+t_{\text{EIG}}$. On the contrary, using the direct VIE solver, we have to find the solution of a different linear system for each of the $n_{\varepsilon_{r}}$ values of permittivity. We accomplished this task by LU factorization using the ScaLAPACK routine called PZGETRF [44]. Each inversion requires a time $t_{\text{FCT}}=11m$. Therefore, the total time is $t_{\text{VIE}}=t_{\text{ASM}}+n_{\varepsilon_{r}}\times t_{\text{FCT}}$. The break-even value of $n_{\varepsilon_{r}}$, in correspondence of which the MIM-VIE method becomes computationally advantageous compared to a direct VIE solution, is $n_{\varepsilon_{r}}\approx 4$. It is also worth to point out that if we are interested in finding the scattering response of the object for just one value of permittivity, i.e. $n_{\varepsilon_{r}}=1$, the MIM-VIE solution is roughly $1.23$ times slower than the direct VIE solution. This analysis clearly demonstrates that the VIE-MIM approach is particularly suited for optimizing the the material permittivity.

IV-B Cylinder

We now consider a cylinder of diameter $\lambda$ and height $H=0.5\lambda$. The finite elements mesh (hexahedral elements) is shown in Fig. 4, whose details are reported in the figure caption.

In Fig. 5 we show with black circles the quantity $\gamma_{h}=\sigma_{h}^{-1}+1$. We used a tolerance $\tau=10^{-2}$, and by using a filter with $\xi=10^{-3}$ we only retain $h_{max}=208$ eigenvalues. Only the $\gamma_{h}$’s belonging to the box of the complex plane $\left[-10,10\right]\times\left[-6,0\right]$ are shown.

Next, we investigate the convergence of the scattering cross section $C_{sca}$ by varying the number of considered MIMs. In particular, in Fig. 6 we plot the $C_{sca}$ of the investigated homogeneous cylinder as a function of the real part of its permittivity $\varepsilon_{r}^{\prime}\in\left[-10,10\right]$, whereas the imaginary part is fixed to the value $\varepsilon_{r}^{\prime\prime}=0.1$. The cylinder is excited by the plane wave of Eq. 20 with $\lambda=1m$. As reference solution we use the solution of the Poggio-Miller-Chang-Harrington-Wu-Tsai (PMCHWT) Surface Integral Formulation of the Maxwell’s equations [6] by RWG basis functions using a triangular mesh with $N_{p}=1373$ nodes, $N_{t}=2742$ triangular elements, $N_{e}=4113$ edges, corresponding to $N_{DoF}=8226$ degrees of freedom. First, we computed $C_{sca}$ using the VIE approach with $h_{max}=34$ (green line). As for the case of a sphere, we notice a good agreement only within a small interval centred at $\varepsilon_{r}^{\prime}=1$. The convergence improves when $h_{max}=107$ modes are considered. Finally, by using $h_{max}=208$ MIMs we can appreciate a very good agreement. Only for negative values of $\varepsilon_{r}^{\prime}$, the two plasmonic resonances appear to be slightly shifted with respect to the surface integral formulation. We plot in Fig. 7 the ten modes that play the most relevant role in correspondence of the peak at $\varepsilon_{r}=4.135+0.1j$. The values of $\gamma_{h}$ associated to these eigenvalues are tabulated in Tab. I, and also plotted with red stars in Fig. 5. By considering exclusively these 10 modes in Eq. 18, we would only commit an error of $9\%$ in the estimation of $C_{sca}$.

Then, we investigate the convergence of the radiation pattern intensity $\left|\mathbf{E}_{S}^{\infty}\left(\theta,\phi\right)\right|^{2}$ as a function of the number of considered MIMs. Specifically, we plot in Fig. 8 the quantity $\left|\mathbf{E}_{S}^{\infty}\left(\theta,\phi\right)\right|^{2}$ as a function of the inclination angle $\theta$ in the $zx$ plane (a) and in the $yz$ plane (b), for different values of $h_{max}$. The cylinder is excited by an $x$-polarized plane wave of unit intensity, incoming from the negative $z$- axis. Also in this case, the cylinder is excited by the plane wave of Eq. 20 with $\lambda=1m$. We assumed a permittivity $\varepsilon_{r}=5+0.1i$, and we performed the calculation using Eq. 18 with $\tau=10^{-2}$ and by the PMCHWT formulation. The green line corresponds to the MIM-VIE solution computed using $h_{max}=34$ modes. It is apparent that we get a very poor agreement with the SIE solution for both the investigated scattering planes. By increasing the number of modes to $h_{max}=107$ (blue line), we obtain convergence toward the SIE solution only for the $xz$ plane, as the other plane still shows lack of convergence. Finally, we obtain very good agreement for both the planes for $h_{max}=208$.

Similarly to what we have done for a sphere, we now calculate the break-even value of $n_{\varepsilon_{r}}$, in correspondence of which the MIM-VIE method becomes computationally advantageous compared to a direct VIE solution for the considered cylinder using $\tau=10^{-2}$. In this case, the cumulative time for assembling both matrices $\bf K$ and $\bf R$ on $25$ processors is the same for both approaches, namely $t_{\text{ASM}}=4h\,10m$. The total time needed for the calculation of the sets of eigenvalues and corresponding MIMs is $t_{\text{EIG}}=1h\,30m$. For the direct VIE solution, the factorization time $t_{\text{FCT}}=31m$. Therefore we get $n_{\varepsilon_{r}}\approx 3$. In conclusion, if we need to calculate the scattering from more than three cylinders having the geometry specified above but three different material compositions the MIM-VIE approach is computational advantageous compared to the direct VIE solver.

V Conclusions

We introduced a volume integral formulation for the calculation of a set of modes independent of the permittivity of the scatterer, which we have denoted as Material Independent Modes (MIMs). The solution of the scattering problem by a homogeneous object can be expanded in terms of MIMs, leading to a natural separation of the contributions of the permittivity and of the geometry. Specifically, the expansion coefficients are simple rational functions of the permittivity, whereas the MIMs depend solely on the geometry.

The calculation of the eigenvalues associated to the MIMs of a sphere is validated against the approach of Ref. [18]. We showed through numerical examples that for objects of linear dimension comparable to the incident wavelength only a moderate number of MIM is needed to accurately describe the scattering cross section for a wide range of permittivities and the radiation pattern. This fact has been shown by investigating the scattering from a sphere and a cylinder of size comparable to the incident wavelength and validating the results exploiting the standard Mie theory and the PMCHWT surface integral formulation.

In addition, we show that, exploiting the fact that the MIMs are independent of the permittivity of the scatterer, the representation of the solution in terms of them may be computationally advantageous compared to direct solution of the volume integral equation when the scattered field has to be computed for many different values of the dielectric permittivity, given the size and shape of the dielectric body.

Eventually, it is also worth to point out that the application of the present method is limited to objects of size comparable to the wavelength due to its memory requirements and computational burden, unless suitable parallel sparsification techniques are used.