Analysis of composite nanoparticles with surface integral equations and the multilevel fast multipole algorithm

Computational simulations of nanoparticles at optical frequencies can be extremely useful by enabling detailed analysis of novel designs before their actual realizations and preventing waste of sources and time required for building and testing prototypes. Compared to radio and microwave frequencies, however, simulations at optical frequencies are usually difficult. For example, highly conductive metals that are simply modeled as perfect electric conductors (PECs) at lower frequencies become penetrable at higher frequencies [1]. In fact, nanoparticles often involve different parts, such as lossy or lossless dielectric, magnetic, perfectly conducting, and plasmonic [2], depending on frequency. In some cases, a part of an overall structure may need to be modeled as a double-negative (DNG) medium with negative effective permittivity and permeability. This letter presents a rigorous solver for the solution of those structures with coexisting multiple regions with diverse electromagnetic properties.

In the literature, simulations of nanoparticles are dominated by volume formulations. Unfortunately, these formulations require volumetric discretizations of structures and host media, and they are usually limited to relatively small problems and/or periodic structures with infinite dimensions. A recently popular approach is extending well-developed surface-integral-equation techniques for radio and microwave frequencies to optical frequencies. Various surface formulations have already been used for the efficient analysis of left-handed materials [3–5]. In a recent work [6], composite structures involving DNG and plasmonic regions are analyzed rigorously with the electric and magnetic current combined-field integral equation (JMCFIE) [7] and the adaptive integral method. It is also shown [8] that single-body plasmonic structures can easily be investigated using JMCFIE and the multilevel fast multipole algorithm (MLFMA) [9]. In this work, we employ JMCFIE and a parallel implementation of MLFMA to analyze composite nanoparticles with diverse material properties. We also consider the accuracy of simulations in respect to the formulation, based on our previous experience at lower frequencies.

Consider time-harmonic electromagnetic fields in a three-dimensional space involving U + 1 different homogeneous regions D_u for u = 0,1,2,...,U. Let D₀ be a free space extending to infinity (host medium). Each region can be characterized by complex permittivity _u = _0u,r and complex permeability μ_u = μ₀μ_u,r. Such an electromagnetics problem can be formulated with JMCFIE and converted into a 2N × 2N matrix equation via discretization. As opposed to volume formulations, only the interfaces between different regions (surfaces) are discretized using edge-based basis and testing functions. In an iterative solution with MLFMA, a matrix–vector multiplication can be written as

$$\begin{eqnarray} &&\boldsymbol{y}=\bar {\boldsymbol{Z}}\boldsymbol{\cdot }\boldsymbol{x}=\sum _{u=0}^{U}\bar {\boldsymbol{Z}}_{\mathrm{NF}}^{u}\boldsymbol{\cdot }\boldsymbol{x}+\sum _{u=0}^{U}\bar {\boldsymbol{Z}}_{\mathrm{FF}}^{u}\boldsymbol{\cdot }\boldsymbol{x}, \end{eqnarray} \tag{ 1 }$$

where ${\bar {\boldsymbol{Z}}}_{\mathrm{NF}}^{u}$ and ${\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u}$ represent near-field and far-field interactions, respectively. Near-field interactions, which are between nearly located basis and testing functions, are calculated directly and combined in a single sparse matrix with (N) nonzero elements. On the other hand, far-field interactions are performed efficiently using the factorization and diagonalization of the Green's function [9]. Specifically, for each region [10],

$$\begin{eqnarray} {\boldsymbol{y}}_{\mathrm{FF}}^{u}={\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u}\boldsymbol{\cdot }\boldsymbol{x}=\left[\begin{array}{@{}cc@{}} \displaystyle {\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u(11)}&\displaystyle {\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u(12)}\\ \displaystyle {\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u(21)}&\displaystyle {\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u(22)} \end{array}\right]\boldsymbol{\cdot }\left[\begin{array}{@{}c@{}} \displaystyle {\boldsymbol{x}}^{(1)}\\ \displaystyle {\boldsymbol{x}}^{(2)} \end{array}\right], \end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray} &&{\bar {\boldsymbol{Z}}}_{\mathrm{FF}}^{u(a b)}[m,n]\approx \int {\mathrm{d}}^{2}\hat {\boldsymbol{k}}{\boldsymbol{F}}_{m}^{u(a b)}(\hat {\boldsymbol{k}})\boldsymbol{\cdot }\tau (\hat {\boldsymbol{k}}){\boldsymbol{R}}_{n}^{u(a b)}(\hat {\boldsymbol{k}}). \end{eqnarray} \tag{ 3 }$$

In (3), τ represents a combination of shift and translation operators to transform the radiated field of the nth basis function ${\boldsymbol{R}}_{n}^{u(a b)}$ into an incoming field, which is multiplied with the receiving pattern of the mth testing function ${\boldsymbol{F}}_{m}^{u(a b)}$. Note that a pair of basis and testing functions interact through a region D_u only if they are both located on the surface of that region.

The complexity of MLFMA is (NlogN); hence it enables fast and accurate solutions of large problems (in terms of unknowns) that may not be solved directly. For faster solutions of larger problems, MLFMA can be parallelized efficiently using the hierarchical partitioning strategy [11, 12]. This strategy is based on simultaneous partitioning of discretization elements and their field samples among processors. This way, load-balancing is improved and communications between processors are reduced, as detailed in [11, 12].

JMCFIE is a mixed formulation involving directly (tangentially) and rotationally (normally) tested electric and magnetic fields [7]. Specifically, this formulation can be written as a convex combination of the combined tangential formulation (CTF) and the combined normal formulation (CNF) [13], i.e. 

$$\begin{eqnarray} &&\mathrm{JMCFIE}=\alpha \mathrm{CTF}+(1-\alpha )\mathrm{CNF}, \end{eqnarray} \tag{ 4 }$$

where α∈[0,1] is a combination parameter. Depending on the value of α, well-tested identity operators in CNF make contributions in the matrix equations derived from JMCFIE. These operators improve the conditioning of matrix equations and lead to faster iterative convergences, but unfortunately they behave like an operator with a highly singular kernel such that they become major error sources for low-order discretizations. Increasing the value of α towards unity improves the accuracy, but the efficiency drops due to an increasing number of iterations. This difficult tradeoff can be considered as an inherited disadvantage of using surface formulations.

For numerical comparisons, we consider solutions of scattering problems involving spherical objects. A sphere of radius 500 nm (core) is placed inside another sphere of radius 1 μm (shell). The relative permittivity of the shell is selected as 2.0, whereas the core has diverse electromagnetic properties. Specifically, we consider a dielectric core with a relative permittivity of 4.0, a PEC core with infinite conductivity, a lossy core with a relative permittivity of 4.0 + 3.0i, a plasmonic core with a relative permittivity of  − 4.0 + 3.0i, and a DNG core with a relative permittivity of  − 4.0 + 3.0i and a relative permeability of  − 1.0. The objects are located at the origin in free space and illuminated by plane waves propagating in the z direction with the electric field polarized in the x direction at 600 THz. The scattering problems are formulated with JMCFIE using α = 0.5, α = 0.7, and α = 1.0, which are denoted by JMCFIE (0.5), JMCFIE (0.7), and CTF, respectively. Discretization of the surfaces with the Rao–Wilton–Glisson functions on λ₀/10 triangles, where λ₀ is the wavelength in the host medium, leads to matrix equations involving 51 222 unknowns. Iterative solutions are performed by the biconjugate-gradient-stabilized (BiCGStab) method employing MLFMA with two digits of accuracy. Iterative convergences are also accelerated with block-diagonal preconditioners for JMCFIE (0.5) and JMCFIE (0.7). The number of iterations to reach 10⁻³ residual error is 25–29 for JMCFIE (0.5), 41–54 for JMCFIE (0.7), and 246–360 for CTF. Each solution is performed in 2–15 min on 16 cores of Intel Xeon Nehalem processors with a 2.80 GHz clock rate.

Figure 1 depicts the total electric field obtained with CTF on the x–y plane in the vicinities of the objects. The field distribution inside the core significantly depends on its electromagnetic properties. The effect of the core in the scattered fields is also remarkable in the far zone, as illustrated in figure 2. In this figure, the scattering cross section (SCS) is plotted as a function of the bistatic angle θ on the z–x plane, where θ = 0° and θ = 180° correspond to the forward-scattering and back-scattering directions, respectively. SCS values obtained with different formulations are compared with those obtained via analytical Mie-series solutions. We observe that all solutions agree well with each other.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For more quantitative comparisons, table 1 lists the root-mean-square (RMS) error in the computational SCS values with respect to the analytical SCS values, along with the number of iterations required for the solutions. The tradeoff between the efficiency and accuracy is clearly observed; increasing the value of α improves the accuracy at the cost of more iterations.

Numerical examples show that composite nanoparticles can be analyzed rigorously with surface integral equations and MLFMA. Compared to volume formulations, surface formulations provide more efficient simulations since only interfaces between different regions are discretized. In addition, solutions can easily be accelerated with fast solvers, such as MLFMA. At the same time, similar to simulations at lower frequencies, the type of the formulation and the discretization scheme are critical factors for the efficiency and accuracy of solutions.

This work was supported by the Centre for Numerical Algorithms and Intelligent Software (EPSRC-EP/G036136/1) and by the Engineering and Physical Sciences Research Council (EPSRC) under research grant EP/J007471/1.