# Nanophotonic resonance modes with the nanobem toolbox

Ulrich Hohenester  Nikita Reichelt  Gerhard Unger  Institute of Physics,
University of Graz, Universitätsplatz 5, 8010 Graz, Austria  Institute of
Advanced Mathemathics, Technical University of Graz, Steyrergasse 30, 8010
Graz, Austria

###### Abstract

nanobem is a matlab toolbox for the solution of Maxwell’s equations for
nanophotonic systems and the computation of resonance modes, sometimes also
referred to as quasinormal modes or resonance states. It is based on a
Galerkin scheme for the boundary element method, using Raviart-Thomas shape
elements for the representation of the tangential electromagnetic fields at
the particle boundary. The toolbox is written in an object-oriented manner
with the focus on clarity rather than speed, and has been developed and tested
for small to intermediate problems with a few thousand boundary elements. The
computation of the resonance modes uses the contour integral method of Beyn.

###### keywords:

quasinormal modes, boundary element method, Maxwell’s equations, nanophotonics

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: nanobem
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: (if available)
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions(please choose one): GNU General Public License 3
Programming language: Matlab
Supplementary material:
Nature of problem: Solve Maxwell’s equations and compute resonance modes for
optical resonators and nanophotonic systems with linear, homogeneous, and
local materials separated by abrupt interfaces.
Solution method: Galerkin implementation of boundary element method approach
using Raviart-Thomas shape elements, and contour integral method for the
computation of nanophotonic resonance modes.
Additional comments including restrictions and unusual features: Toolbox has
been developed and tested for small to intermediate problems with a few
thousand boundary elements.

##  1 Introduction

Resonance modes, sometimes also referred to as quasinormal modes or resonant
states, have recently received considerable interest in the field of
nanophotonics [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. In optical resonators and
cavities, light can be confined on a length scale comparable to the
wavelength. By binding light to coherent excitations of nanostructures, such
as surface plasmons or surface phonons, one can achieve light confinement down
to extreme sub-wavelength volumes [11, 12]. Such optical resonators play an
important role in science and engineering, for instance for microwave
resonators, semiconductor lasers, or optical microcavities and plasmonic nano-
cavities.

From a computational perspective, the simulation of optical resonators builds
on the solution of Maxwell’s equations. Let $u$ denote the degrees of freedom
for the representation of the (total or scattered) electromagnetic fields,
e.g. curl-conforming Nedelec elements for the finite element method approach
[13, 14] or Raviart-Thomas shape elements for the boundary element method
approach [15, 16, 12], and $q_{\rm inc}$ the incoming electromagnetic fields,
such as a planewave excitation of the optical resonator. For the linear
material properties of our present concern, $u$ is related to an excitation
$q_{\rm inc}$ oscillating with angular frequency $\omega$ through

$A(\omega)u=q_{\rm inc}\,,$ (1)

where the transmission matrix $A(\omega)$ has to be derived from Maxwell’s
equations subject to the proper boundary conditions between different
materials. Solving Eq. (1) for the unknowns $u$ provides us with the solutions
of Maxwell’s equations. Instead of solving Eq. (1) directly, we can also seek
for the resonance frequencies $\omega_{k}$ and modes $v_{k}$ of the
transmission problem, which are defined such that they can persist in absence
of an external excitation,

$A(\omega_{k})v_{k}=0\,.$ (2)

For absorbing media or open cavities, Eq. (2) can only be fulfilled for
complex frequencies
$\omega_{k}=\omega_{k}^{\prime}+i\omega_{k}^{\prime\prime}$, where the
material and optical losses are compensated by the imaginary part
$\omega_{k}^{\prime\prime}$. In a second step, the electromagnetic response
$u$ to the incoming fields $q_{\rm inc}$ is expressed in terms of these
resonance modes

$u\approx\sum_{k}C_{k}(\omega)v_{k}\,,$ (3)

Figure 1: Schematics of the boundary integral method (bim) and the resonance
modes. (a) We consider a setup where a nanoparticle with homogeneous material
parameters mat1 is embedded in a medium with material parameters mat2. The
boundary separating the two media is denoted with $\partial\Omega$, the outer
surface normal of the boundary is $\hat{\bm{n}}$. (b) In bim, the response of
the nanoparticle to some external perturbation (e.g. planewave excitation or
oscillating dipole) can be described in terms of the tangential
electromagnetic fields $\bm{u}_{E,H}$ at the particle boundary, fields away
from the boundary can be computed from the representation formula. (c) In the
resonance mode approach we don’t solve Maxwell’s equation directly, but first
determine the resonance frequencies $\omega_{k}$ and modes $v_{k}$ of the
nanoparticle, and then expand the tangential fields $\bm{u}_{E,H}$ in terms of
these modes. The second step is usually referred to as the “reconstruction
problem”.

where the determination of the expansion coefficients $C_{k}(\omega)$ is
usually referred to as the reconstruction problem [9], see also Fig. 1 for a
schematic representation. The advantage of the resonance mode solution (3) is
often the speedup in comparison to the direct solution of Eq. (1), provided
that the resonance modes are at hand. Probably even more important is the
interpretation of the solutions in terms of intuitive and physically
transparent resonance modes, which also govern the physics of optical
resonators. These modes can be additionally submitted to a quantization
procedure [17], which is beneficial for quantum optical investigations.

The consideration of complex frequencies in Eq. (2) for the evaluation of the
resonance modes can be complicated for volume-based schemes, such as the
finite element method approach, because of the increase of the electromagnetic
fields away from the resonators. Also the damping of such fields in perfectly
matched layers surrounding the simulation volume can be critical. Although
these problems can be overcome in principle, approaches such as the boundary
element method (bem) avoid these difficulties because the solution vector $u$
only involves the tangential electromagnetic fields at the boundaries
separating different materials. As has been demonstrated in [18], the
determination of resonance modes using a bem approach is relatively simple and
free from any conceptional or computational difficulties.

In this paper we present the nanobem toolbox for the computation of resonance
modes using a bem approach based on a Galerkin scheme [15, 16, 12]. The
toolbox is implemented in matlab in an object-oriented manner, and comes along
with a number of demo files together with detailed help pages. The main
purpose of the toolbox is to provide a fast, efficient, and flexible toolkit
for computing and analyzing resonance modes for simple nanophotonic
structures, such as dielectric or metallic nanospheres or nanorods, and to
allow for easy modifications and extensions in case of future improvements. At
the same time, the toolbox supersedes the mnpbem toolbox developed by one of
the authors of this paper [19, 20, 21], which is based on a collocation
approach for the electrodynamic potentials [22]. The user interface and the
program structures are similar for the two toolboxes, but all internal classes
have been completely rewritten for the nanobem toolbox to accommodate for the
Galerkin approach using the electromagnetic fields as basic quantities.

We have organized this work as follows. In Sec. 2 we show how to set up the
toolbox and how to run simple simulations for nanophotonic systems. The theory
underlying the bem approach and the contour integral method (cim) for
computing the resonance modes is presented in Sec. 3, implementation details
of the toolbox are discussed in Sec. 4. Finally, in Sec. 5 we list a few dos
and don’ts for the cim approach and the computation of resonance modes. Some
theoretical and computational details are postponed to the Appendices.

##  2 Getting started

In this section, we briefly describe how to set up the toolbox, and how to run
simple boundary element method (bem) simulations and simulations using
resonance modes. To install the toolbox, one must add the path of the main
directory nanobemdir of the nanobem toolbox as well as the paths of all
subdirectories to the matlab search path. This can be done, for instance,
through

>> addpath(genpath(nanobemdir));To set up the help pages, one must once change
to the main directory of the nanobem toolbox and run makehelp

>> cd nanobemdir;>> makehelp;Once this is done, the help pages, which provide
detailed information about the toolbox, are available in the matlab help
browser.

###  2.1 A simple BEM example

Figure 1(a) shows a typical bem example, where a gold nanosphere with a
diameter of 50 nm is embedded in a dielectric medium representative for water,
with a refractive index of $1.33$. In the following we briefly discuss a
simulation where the optical cross section is computed. The corresponding
matlab files are available as `demoretspec01.m` and `democim01.m` in the help
pages of the toolbox. The nanophotonic environment is set up through

mat1 = Material( 1.33 ^ 2, 1 ); % material properties for watermat2 =
Material( epsdrude( ’Au’ ), 1 ); % material properties for goldmat = [ mat1,
mat2 ]; % unique material vectordiameter = 50; % diameter of sphere in
nanometersp = trisphere( 144, diameter ); % triangulated nanosphere
boundarytau = BoundaryEdge( mat, p, [ 2, 1 ] ); % boundary elements including
shape elementsIn the first two lines we set up material objects
`Material(eps,mu)`, where `eps` is either a permittivity function or constant
and `mu` a permeability function or constant. The `epsdrude` object provides a
Drude dielectric function representative for gold. We finally specify a
material vector `mat` for the nanophotonic environment, where the first entry
specifies the material properties of the background medium, and the other
entries specify the material properties of the particles. From here on, the
same material vector should be passed to all toolbox functions. See also Sec.
3.2, B, and the help pages for a more detailed description. We next generate a
triangulation of the nanosphere boundary with 144 vertices. With

>> plot( p, ’EdgeColor’, ’b’ );>> plot( p, ’EdgeColor’, ’b’, ’nvec’, 1 ); %
show outer surface normalswe can plot the triangulated boundary and the outer
surface normals. Finally, we create a `BoundaryEdge` object with the material
vector, the triangulated boundary, and the information that the material at
the boundary inside (measured with respect to the outer surface normal) is
`mat(2)` and the material at the boundary outside is `mat(1)`. The boundary
elements `tau` are additionally equipped with shape elements for the
representation of the tangential electromagnetic fields, see also Fig. 5, as
will be discussed in more detail below.

Once the nanophotonic environment is set up, we can run a simulation by
specifying the external excitation, for instance a planewave excitation,
initializing the bem solver, and solving the bem equations for a number of
light wavelengths $\lambda$. In the toolbox we use the free-space wavenumber
$k_{0}=\nicefrac{{2\pi}}{{\lambda}}$ instead.

pol = [ 1, 0, 0 ]; % light polarizationdir = [ 0, 0, 1 ]; % light propagation
directionexc = galerkin.planewave( pol, dir ); % planewave excitationbem =
galerkin.bemsolver( tau ); % BEM solverlambda = linspace( 400, 800, 20 );k0 =
2 * pi ./ lambda; % wavenumber of light in vacuumcext = zeros( size( k0 ) ); %
initialize array for extinction cross sectionfor i = 1 : numel( k0 ) % loop
over wavenumbers sol = bem \ exc( tau, k0( i ) ); % solve BEM equations cext(
i ) = extinction( exc, sol ); % compute exctinction cross sectionendThe above
code is quite self-explanatory and we will discuss the various commands in
more detail in Sec. 3.2. When solving the BEM equations, we obtain a solution
object `sol` that stores the tangential electromagnetic fields at the sphere
boundary. From these fields all other electromagnetic quantities can be
determined, such as the extinction cross section. We can finally plot the
extinction cross section

>> plot( lambda, cext );to obtain the extinction cross section shown in Fig.
2(a). See also Fig. 3 for the example of a dielectric nanosphere with a Drude-
Lorentz dielectric function representative for silicon.

Figure 2: (a) Extinction cross sections for gold nanosphere (diameter 50 nm)
embedded in water (black line) and air (gray line), as obtained from the
solutions of the bem equations (circles) and the resonance mode expansion
($+$, $\times$). (b) Same as (a) but using for the abscissa the photon energy
in eV rather than the light wavelength in nanometers. (c) Extinction cross
sections for a gold nanoellipsoid with a long axis of 100 nm and short axis of
50 nm. (d,e) cim contour in the complex frequency plane (solid lines) and
resonance mode energies $\omega_{k}$ (symbols). Note that for small
nanoparticles only modes with a non-vanishing dipole moment can be excited
optically, with a similar behavior for larger nanoparticles where typically
only a few modes contribute noticeably to the optical spectra. Figure 3: Same
as Fig. 2 but for a dielectric nanosphere with a diameter of 150 nm. For the
nanosphere we consider a dielectric function representative for silicon, with
material parameters taken from Ref. [23], which is embedded in a background
medium with a refractive index of one. (a,b) The bem and cim results using the
42 resonance modes shown in panel (d) agree well for smaller photon energies,
say below 2.2 eV, but start to deviate consistently for larger energies. (c,e)
By considering in cim the same resonance modes as in panel (d) together with a
non-resonant background (dashed line), one obtains almost perfect agreement
for all photon energies. For details see discussion in Sec. 4.7.

###  2.2 Resonance modes and the cim solver

We next show how to compute the optical cross sections using the resonance
mode expansion. Details of our theoretical approach and of the working
principle for the contour integral method (cim) approach will be presented in
Secs. 3 and 4, respectively. In brief, we are seeking for a description of the
nanophotonic environment in terms of (complex) resonance frequencies
$\omega_{k}$ and modes, as schematically shown in Fig. 1(c). To find these
resonances, the user must define a contour in the complex frequency space, see
ellipse in Figs. 2(d,e), and the cim algorithm then returns the resonances
located inside this region in an (almost) black-box manner. We will discuss
later in Sec. 5 that the approach only works for permittivity and permeability
functions that can be extended to the complex frequency space, such as the
Drude dielectric function used in the example of Sec. 2.1, and the contour
must not cross branch cuts of the square root of these functions, and must
also exclude accumulation points of resonances. Within the toolbox, the cim
solver is set up with

cim = galerkin.cimsolver( tau );cim.contour = cimbase.ellipse( [ 1.5, 2.6 ],
0.2, 60 );

Figure 4: Graphical user interface of the tolselect toolbox function for the
computation of the resonance modes using the approach of Beyn [24]. The left
panel shows the singular values of the Beyn approach, see Sec. 4.7 for
details, and the user can interactively select those resonance modes that have
a sufficiently large singular value. This step is somewhat critical, because
several things can go wrong in the cim approach, as briefly discussed in Sec.
5. In the panel of the left matlab figure three special icons have been added:
when the ellipse is pressed, the figure shown in the right panel opens up and
displays the contour in the complex plane together with the resonance
frequencies of the modes selected by the user. When the green triangle
$\triangleright$ is pressed, the cim calculation continues, when the red
circle $\bullet$ is pressed, the cim calculation stops.

In the second line we define an ellipse in the complex frequency plane
$\omega=\omega^{\prime}+i\omega^{\prime\prime}$, where the axis along
$\omega^{\prime}$ goes from 1.5 eV to 2.6 eV, and the axis for
$\omega^{\prime\prime}$ goes from $-0.2$ eV to $0.2$ eV. For the evaluation of
the resonance modes we use 60 quadrature points, following the approach of
Beyn [24] to be discussed in Sec. 4.7. We then solve the bem equations in the
complex frequency plane, and the user must select the resonance modes to be
kept in the simulations from a singular value decomposition [24]. This is
accomplished with

data = eval1( cim ); % solve BEM equations along contourtol = tolselect( cim,
data, ’tol’, 1e-2 ); % select modes according to singular valuesUpon execution
of the `tolselect` command, the window shown in the left panel of Fig. 4 opens
up. It displays the singular values returned from the cim algorithm. For
typical problems, the resonance modes have the largest singular values and are
separated by a sufficiently large gap from those modes which should be
discarded because they do not noticeably contribute to the resonance mode
expansion. We will comment on this point in more detail in Sec. 5. The user
can select the relevant modes interactively with the mouse cursor (red line)
and by clicking the left mouse button. Once the modes (red dots in left panel)
are selected, we can visualize the resonance frequencies in the complex plane
by pressing the ellipse icon in the matlab figure panel. This opens the figure
shown in the right panel of Fig. 4. To continue from here on with the cim
simulation, the user should press the green triangle icon $\triangleright$. If
something has gone wrong with the computation of the resonance modes, the user
should press the red circle $\bullet$. In the calling program the progress of
the cim calculation is then handled with

if isempty( tol ), return; endcim = eval2( cim, data, ’tol’, tol );The first
line makes sure that the user has not pressed the red circle icon $\bullet$,
otherwise the calculation is stopped. In the second line we keep the modes
with a sufficiently large singular value and normalize the modes, using the
approach discussed in Sec. 3. Once the cim solver has been initialized, it can
be used in the same manner as the previously described bem solver. For
instance, using the excitation object `exc` for plane waves we can compute the
bem solution through

sol = bem \ exc( tau, k0 );This step is usually fast because all relevant
information about the electromagnetic response to an external excitation is
already contained in the resonance modes. Yet, the description is often
extremely accurate, as can be seen for instance by comparing the results of
the bem and cim simulations in Figs. 2 and 3, which are almost
indistinguishable.

##  3 Theory

In this section, we first introduce to the theory of the boundary integral
method approach and its numerical implementation using boundary elements, and
then discuss the basic ideas behind resonance modes and resonance mode
expansions.

###  3.1 Boundary integral method

Within the boundary integral method one considers one or several particles
embedded in a background medium, as depicted in Fig. 1(a). In the following we
discuss for simplicity the case of a single particle only, but our approach
can be easily generalized and the nanobem toolbox implements general
geometries with possibly coupled or coated particles. Our main assumption
about the system under study is that the electromagnetic response of the
particle and the background medium is linear and can be described in terms of
homogeneous, isotropic, and local material parameters. We denote the
permittivity and permeability of the particle with $\varepsilon_{1}$ and
$\mu_{1}$, respectively, and the material properties of the background medium
with $\varepsilon_{2}$, $\mu_{2}$. The particle volume is $\Omega$, the
particle boundary is $\partial\Omega$, and the outer surface normal pointing
away from the particle is $\hat{\bm{n}}$. The boundary integral method, to be
discussed below, is a solution scheme for Maxwell’s equations that only
involves the tangential electromagnetic fields at the particle boundary
$\bm{s}\in\partial\Omega$,

$\bm{u}_{E}(\bm{s})=\hat{\bm{n}}\times\bm{E}(\bm{s})\,,\qquad\bm{u}_{H}(\bm{s})=\hat{\bm{n}}\times\bm{H}(\bm{s})\,.$
(4)

We consider the case of incoming electromagnetic fields $\bm{E}_{2}^{\rm
inc}$, $\bm{H}_{2}^{\rm inc}$ oscillating with a given angular frequency
$\omega$, which are either solutions of the homogeneous Maxwell’s equations,
e.g., plane waves, or are produced by current sources located within the
background medium, e.g., oscillating dipoles. These incoming fields would be
the solutions of Maxwell’s equations in absence of the particle. In presence
of the particle we additionally have to consider the scattered fields
$\bm{E}_{1}^{\rm sca}$, $\bm{H}_{1}^{\rm sca}$ and $\bm{E}_{2}^{\rm sca}$,
$\bm{H}_{2}^{\rm sca}$ inside and outside the particle, respectively, which
have to be chosen such that the boundary conditions of Maxwell’s equations
hold at the particle boundary. The total electromagnetic fields are the sum of
incoming and scattered fields. As discussed in Refs. [15, 12] and briefly
outlined in A, for positions $\bm{r}$ inside the particle the electromagnetic
fields can be obtained from the representation formulas (we use SI units
throughout)

$\displaystyle\bm{E}(\bm{r})$ $\displaystyle=$
$\displaystyle-i\omega\mu_{1}\bigl{[}\mathbb{S}_{1}\,\bm{u}_{H}\bigr{]}(\bm{r})+\bigl{[}\mathbb{D}_{1}\,\bm{u}_{E}\bigr{]}(\bm{r})$
$\displaystyle\bm{H}(\bm{r})$ $\displaystyle=$
$\displaystyle+i\omega\varepsilon_{1}\bigl{[}\mathbb{S}_{1}\,\bm{u}_{E}\bigr{]}(\bm{r})+\bigl{[}\mathbb{D}_{1}\,\bm{u}_{H}\bigr{]}(\bm{r})\,.\qquad$
(5a) Similarly, for positions $\bm{r}$ outside the particle the
electromagnetic fields can be obtained from $\displaystyle\bm{E}(\bm{r})$
$\displaystyle=$ $\displaystyle\bm{E}_{2}^{\rm
inc}(\bm{r})+i\omega\mu_{2}\bigl{[}\mathbb{S}_{2}\,\bm{u}_{H}\bigr{]}(\bm{r})-\bigl{[}\mathbb{D}_{2}\,\bm{u}_{E}\bigr{]}(\bm{r})$
$\displaystyle\bm{H}(\bm{r})$ $\displaystyle=$ $\displaystyle\bm{H}_{2}^{\rm
inc}(\bm{r})-i\omega\varepsilon_{2}\bigl{[}\mathbb{S}_{2}\,\bm{u}_{E}\bigr{]}(\bm{r})-\bigl{[}\mathbb{D}_{2}\,\bm{u}_{H}\bigr{]}(\bm{r})\,.\qquad$
(5b)

Here $\bigl{[}\mathbb{S}\,\bm{u}_{E,H}\bigr{]}(\bm{r})$ and
$\bigl{[}\mathbb{D}\,\bm{u}_{E,H}\bigr{]}(\bm{r})$ are the single and double
layer potentials, for details see A, which propagate the tangential
electromagnetic fields $\bm{u}_{E,H}$ from the particle boundary to position
$\bm{r}$. Note that owing to the boundary conditions of Maxwell’s equations
$\bm{u}_{E,H}$ are the same inside and outside the particle. Eq. (5) is
reminiscent of Huygen’s principle which propagates the fields at the wavefront
(here $\bm{u}_{E,H}$) to another position in space (here $\bm{r}$). Thus, once
$\bm{u}_{E,H}$ are known at the boundary, we can compute the electromagnetic
fields everywhere else using the representation formulas of Eq. (5).

The tangential fields $\bm{u}_{E,H}$ can be obtained by letting $\bm{r}$ in
Eq. (5) approach the boundary of the particle, and by exploiting the boundary
conditions of Maxwell’s equations. While this limit can be performed safely
for the single layer potential,
$\lim_{\bm{r}\to\bm{s}}\hat{\bm{n}}\times\bigl{[}\mathbb{S}\,\bm{u}_{E,H}\bigr{]}(\bm{r})=\hat{\bm{n}}\times\bigl{[}\mathbb{S}\,\bm{u}_{E,H}\bigr{]}(\bm{s})$,
in the evaluation of the double layer potential we have to be careful on
whether we approach the boundary from the inside or outside [12]

$\lim_{\bm{r}\to\bm{s}}\hat{\bm{n}}\times\bigl{[}\mathbb{D}\,\bm{u}_{E,H}\bigr{]}(\bm{r})=\pm\mbox{$\frac{1}{2}$}\bm{u}_{E,H}(\bm{s})+\hat{\bm{n}}\times\bigl{[}\mathbb{D}\,\bm{u}_{E,H}\bigr{]}(\bm{s})\,.$
(6)

Here the positive sign has to be taken for the limit from the inside, and the
negative sign for the limit from the outside. Thus, if we consider in Eq. (5a)
the tangential fields $\hat{\bm{n}}\times\bm{E}$, $\hat{\bm{n}}\times\bm{H}$
and approach the boundary from the inside, we get

$\displaystyle\mbox{$\frac{1}{2}$}\bm{u}_{E}(\bm{s})$ $\displaystyle=$
$\displaystyle-i\omega\mu_{1}\hat{\bm{n}}\times\bigl{[}\mathbb{S}_{1}\,\bm{u}_{H}\bigr{]}(\bm{s})+\hat{\bm{n}}\times\bigl{[}\mathbb{D}_{1}\,\bm{u}_{E}\bigr{]}(\bm{s})$
$\displaystyle\mbox{$\frac{1}{2}$}\bm{u}_{H}(\bm{s})$ $\displaystyle=$
$\displaystyle+i\omega\varepsilon_{1}\hat{\bm{n}}\times\bigl{[}\mathbb{S}_{1}\,\bm{u}_{E}\bigr{]}(\bm{s})+\hat{\bm{n}}\times\bigl{[}\mathbb{D}_{1}\,\bm{u}_{H}\bigr{]}(\bm{s})\,.\qquad$
(7a) Similarly, we obtain from Eq. (5b)
$\displaystyle\mbox{$\frac{1}{2}$}\bm{u}_{E}(\bm{s})$ $\displaystyle=$
$\displaystyle\hat{\bm{n}}\times\bm{E}_{2}^{\rm
inc}(\bm{s})+i\omega\mu_{2}\hat{\bm{n}}\times\bigl{[}\mathbb{S}_{2}\,\bm{u}_{H}\bigr{]}(\bm{s})-\hat{\bm{n}}\times\bigl{[}\mathbb{D}_{2}\,\bm{u}_{E}\bigr{]}(\bm{s})$
$\displaystyle\mbox{$\frac{1}{2}$}\bm{u}_{H}(\bm{s})$ $\displaystyle=$
$\displaystyle\hat{\bm{n}}\times\bm{H}_{2}^{\rm
inc}(\bm{s})-i\omega\varepsilon_{2}\hat{\bm{n}}\times\bigl{[}\mathbb{S}_{2}\,\bm{u}_{E}\bigr{]}(\bm{s})-\hat{\bm{n}}\times\bigl{[}\mathbb{D}_{2}\,\bm{u}_{H}\bigr{]}(\bm{s})\,.\qquad$
(7b)

For given incoming fields $\bm{E}_{2}^{\rm inc}$, $\bm{H}_{2}^{\rm inc}$, we
must extract from the two sets of equations (7) the unknown $\bm{u}_{E,H}$. In
the Poggio-Miller-Chang-Harrington-Wu-Tsai formulation [25, 26, 27] one
subtracts the two sets of equations to arrive at [12]

$\displaystyle\hat{\bm{n}}\times\Bigl{(}{\mathbb{D}}_{1}[\bm{u}_{E}](\bm{s})+{\mathbb{D}}_{2}[\bm{u}_{E}](\bm{s})\Bigr{)}-\hat{\bm{n}}\times\Bigl{(}i\omega\mu_{1}{\mathbb{S}}_{1}[\bm{u}_{H}](\bm{s})+i\omega\mu_{2}{\mathbb{S}}_{2}[\bm{u}_{H}](\bm{s})\Bigr{)}$
$\displaystyle=$ $\displaystyle\hat{\bm{n}}\times\bm{E}_{2}^{\rm inc}(\bm{s})$
$\displaystyle\hat{\bm{n}}\times\Bigl{(}{\mathbb{D}}_{1}[\bm{u}_{H}](\bm{s})+{\mathbb{D}}_{2}[\bm{u}_{H}](\bm{s})\Bigr{)}+\hat{\bm{n}}\times\Bigl{(}i\omega\varepsilon_{1}{\mathbb{S}}_{1}[\bm{u}_{E}](\bm{s})+i\omega\varepsilon_{2}{\mathbb{S}}_{2}[\bm{u}_{E}](\bm{s})\Bigr{)}$
$\displaystyle=$ $\displaystyle\hat{\bm{n}}\times\bm{H}_{2}^{\rm
inc}(\bm{s})\,.$ (8)

These are the working equations for the boundary integral method approach.

###  3.2 Boundary element method

Figure 5: Discretized particle boundary and Raviart-Thomas shape elements for
boundary element method (bem) approach. (a) In bem, the nanoparticle boundary
is discretized by boundary elements, here of triangular shape. (b) In each
triangle $\tau_{i}$ we introduce Raviart-Thomas shape elements
$f_{ia}(\bm{s})$ for the representation of the tangential electromagnetic
fields. The shape elements are constructed such that the inflow to an edge
equals the outflow to the adjacent triangle. For the representation of the
tangential fields we attach a value
$\normalsize\llbracket{u_{E,H}^{e}}\normalsize\rrbracket$ to each edge and
interpolate the fields in the triangles using the shape elements. Each edge is
labeled by a unique global degree of freedom $\nu$, within each triangle we
use the local degrees of freedom $i,a$, where $a$ denotes the different edges
of triangle $\tau_{i}$.

The boundary element method (bem) approach implements Eq. (3.1) for a
discretized particle boundary

$\partial\Omega=\bigcup_{i=1}^{m}\tau_{i}\,,$ (9)

where $\tau_{i}$ denote the different boundary elements and $m$ is the total
number of elements. See also Fig. 5. In the nanobem toolbox we currently
provide boundary elements of triangular shape only and consider regular mesh
triangulations where each edge of a given triangle is shared by a single
adjacent triangle. In addition to the discretized particle boundary, we also
need a discretization for the tangential electromagnetic field $\bm{u}_{E,H}$.
A convenient choice are the Raviart-Thomas or Rao-Wilton-Glison basis
elements, which guarantee the continuity of the tangential fields when going
from one triangle to an adjacent triangle [15, 12]. Technically, this is done
by assigning to each edge $\nu$ of the discretized particle a value for
$\normalsize\llbracket{u_{E,H}^{e}}\normalsize\rrbracket_{\nu}$ and by using
basis functions $\bm{f}_{\nu}^{e}(\bm{s})$ that are nonzero in the two
adjacent triangles only (a so-called local support), and which are constructed
such that the outflow from one triangle equals the inflow to the other
triangle. For details see [15, 12]. The discretized tangential fields can then
be approximated through

$\bm{u}_{E,H}(\bm{s})=\sum_{\nu=1}^{n}\big{\llbracket}{u_{E,H}^{e}}\big{\rrbracket}_{\nu}\,\bm{f}^{e}_{\nu}(\bm{s})\,,$
(10)

where $n$ is the total number of individual edges that equals the number of
degrees of freedom for the bem approach. Finally, Eq. (3.1) is converted into
a matrix equation using a Galerkin scheme. For the single layer potential we
obtain [12]

$\big{\llbracket}{\mathbb{S}}\big{\rrbracket}_{\nu\nu^{\prime}}=\oint\bm{f}^{e}_{\nu}(\bm{s})\cdot\big{[}\mathbb{S}\,\bm{f}^{e}_{\nu^{\prime}}\big{]}(\bm{s})\,dS\,,$
(11)

where $dS$ denotes the integration over the boundary, with a similar
expression for the double layer potential. Note that because of the local
support of the basis function $\bm{f}_{\nu}^{e}$ the above integral can be
broken down into integrals over triangle pairs. Finally, the working equations
(3.1) of the boundary integral method approach can be approximated through a
matrix equation

$\begin{pmatrix}\normalsize\llbracket{{\mathbb{D}}_{1}+{\mathbb{D}}_{2}}\normalsize\rrbracket&-i\omega\normalsize\llbracket{\mu_{1}{\mathbb{S}}_{1}+\mu_{2}{\mathbb{S}}_{2}}\normalsize\rrbracket\\\
i\omega\normalsize\llbracket{\varepsilon_{1}{\mathbb{S}}_{1}+\varepsilon_{2}{\mathbb{S}}_{2}}\normalsize\rrbracket&\normalsize\llbracket{{\mathbb{D}}_{1}+{\mathbb{D}}_{2}}\normalsize\rrbracket\end{pmatrix}\left(\begin{array}[]{c}\normalsize\llbracket{u_{E}^{e}}\normalsize\rrbracket\\\
\normalsize\llbracket{u_{H}^{e}}\normalsize\rrbracket\end{array}\right)=\left(\begin{array}[]{c}\normalsize\llbracket{q_{E}^{\rm
inc}}\normalsize\rrbracket\\\ \normalsize\llbracket{q_{H}^{\rm
inc}}\normalsize\rrbracket\end{array}\right)\,,$ (12)

with the inhomogeneity of the incoming fields

$\llbracket{q_{E,H}^{\rm
inc}}\rrbracket_{\nu}=\oint\bm{f}_{\nu}^{e}(\bm{s})\cdot\left\\{\begin{array}[]{c}\bm{E}^{\rm
inc}(\bm{s})\\\ \bm{H}^{\rm inc}(\bm{s})\end{array}\right\\}\,dS\,.$ (13)

Details of the implementation of Eq. (12) within the nanobem toolbox will be
given below.

###  3.3 Resonance modes and contour integral method

In many cases of interest, the working equation (12) can be approximately
solved by using resonance modes. We start by rewriting the working equations
in the more compact form

$\begin{pmatrix}A_{11}&A_{12}\\\ A_{21}&A_{22}\\\
\end{pmatrix}\begin{pmatrix}u_{E}\\\ u_{H}\\\
\end{pmatrix}=\begin{pmatrix}q_{E}\\\ q_{H}\\\ \end{pmatrix}\,,$ (14)

or $A(\omega)u=q$ in short, where we have explictly indicated the dependence
of the transmission matrix $A$ on the angular frequency $\omega$ of the
external excitation. A resonance mode is characterized through a (in general
complex) resonance frequency $\omega_{k}$ for which the following expression
holds

$A(\omega_{k})v_{k}=0\,,\qquad v_{k}=\begin{pmatrix}v_{E}\\\ v_{H}\\\
\end{pmatrix}\,.$ (15)

Here $v_{k}$ is the resonance mode associated with $\omega_{k}$. Due to our
formulation of the bem equations, see also A, it is guaranteed that these
modes have the proper boundary conditions of outgoing waves at infinity. The
computation of the resonance modes is particularly simple within the bem
approach, contrary to many other computational schemes such as the finite
element method approach, where the consideration of complex frequencies
requires a careful treatment of the perfectly matched layers in order to
properly damp the outgoing waves. In contrast, within bem the single and
double layer potentials only connect positions on the boundary of the
particle, and $\mathbb{S}$, $\mathbb{D}$ can be safely evaluated for complex
frequencies. In the nanobem toolbox we employ the contour integral method
approach of Beyn [24], which allows computing the eigenmodes in a complete
black box manner: we solve Eq. (14) for random inhomogeneities along a contour
in the complex frequency plane, and the algorithm provides the resonances
located inside the contour. Details of our implementation of the contour
integral method (cim) and its working principles will be presented in Sec.
4.7.

In addition to the right eigenmodes of Eq. (15) there also exist left
eigenmodes $\tilde{\nu}_{k}$ defined through

$\tilde{v}_{k}^{T}A(\omega_{k})=0\,,\qquad\tilde{v}_{k}=\begin{pmatrix}v_{H}\\\
v_{E}\\\ \end{pmatrix}\,,$ (16)

where $\tilde{v}_{k}$ is obtained from $v_{k}$ by exchanging the electric and
magnetic components. This can be easily proven by taking the transpose of Eq.
(15) and noting that all of the submatrices $A_{11}$, $A_{12}$, $A_{21}$,
$A_{22}$ in Eq. (14) are symmetric, and additonally $A_{11}=A_{22}$ holds, as
can be inferred from the matrix in Eq. (12).

In a second step, we approach the reconstruction problem and seek for a
description of the electromagnetic response to some incoming fields in terms
of resonance modes via

$u_{\rm CIM}(\omega)\approx\sum_{k=1}^{K}C_{k}(\omega)\,v_{k}\,,$ (17)

where $K$ is the number of modes under consideration, e.g., those located
inside a given contour in the complex frequency plane, and $C_{k}(\omega)$ are
the expansion coefficients to be determined. In case of a single mode $v_{k}$,
we can expand close to the resonance $\omega\approx\omega_{k}$ the
transmission matrix into a Taylor series and rewrite the bem equation
$A(\omega)u=q$ approximately as

$\Big{(}A(\omega_{k})+A^{\prime}(\omega_{k})\big{(}\omega-\omega_{k})\Big{)}C_{k}(\omega)v_{k}\approx
q\,,$ (18)

where $A^{\prime}(\omega_{k})$ denotes the derivative of the transmission
matrix at the resonance frequency. Multiplying the equation from the left-hand
side with $\tilde{v}_{k}^{T}$ and using Eq. (16), allows us to express the
expansion coefficient as

$C_{k}(\omega)\approx\Big{(}\tilde{v}_{k}^{T}A^{\prime}(\omega_{k})v_{k}\Big{)}^{-1}\frac{\tilde{v}_{k}^{T}q}{\omega-\omega_{k}}\,.$
(19)

The term in parentheses is a normalization constant, and the fraction is a
Lorentzian lineshape function with the resonance frequency $\omega_{k}$ and
the oscillator strength $\tilde{v}_{k}^{T}q$. In case of many resonance modes,
one can show that the expansion of Eq. (17) can be written in the form

$u_{\rm
CIM}(\omega)=\sum_{k=1}^{K}\Big{(}\tilde{v}_{k}^{T}A^{\prime}(\omega_{k})v_{k}\Big{)}^{-1}\frac{\tilde{v}_{k}^{T}q}{\omega-\omega_{k}}v_{k}+h(\omega)\,,$
(20)

where $h(\omega)$ is a holomorphic function that contains contributions from
resonance modes located outside the contour and terms originating from higher-
order contributions of the Taylor expansion of the transmission matrix
$A(\omega)$. Eq. (20) follows from Keldysh’s theorem [29, 30], a central
result of the theory of eigenvalue problems for holomorphic Fredholm operator-
valued functions, which can be directly applied to our bem formulation. In
many cases of interest, the contribution of $h(\omega)$ can be neglected, but
there is no guarantee that such an approach always works. See Sec. 4.7 for a
discussion of this issue and the possibility to consider a non-resonant
background.

##  4 The nanobem toolbox

The nanobem toolbox is a matlab implementation of the Galerkin scheme for the
boundary element method approach using Raviart-Thomas basis functions. It is
written in an object oriented manner with the main emphasis on clarity rather
than speed. We have currently tested the code for relatively small problems
only, say a few thousand boundary elements. In the following we discuss the
various steps of a typical bem simulation in more detail, which require
setting up the particle boundaries and the dielectric environment, specifying
the external perturbation, and solving the bem equations, with a possible post
processing of the final solution vector. Note that the toolbox comes along
with matlab help pages that provide further information about the
implementation details. Table 4 lists the most important classes of the
toolbox.

Table 1: List of the most important classes of the nanobem toolbox. A detailed
description of the various classes is given in Sec. 4 and in B.

Class name |  Description
---|---
Material |  Material properties of media
particle |  Discretization of particle boundary in terms of vertices and faces
BoundaryEdge |  Boundary elements with Raviart-Thomas shape elements
Point |  Point embedded in photonic environment
galerkin.planewave |  Planewave excitation
galerkin.dipole |  Oscillating dipole excitation
galerkin.bemsolver |  bem solver for solution of working equation (12)
galerkin.cimsolver |  Resonance mode solver based on the contour integral method (cim)
galerkin.solution |  Solution vector for $\normalsize\llbracket{u_{E,H}}\normalsize\rrbracket$ of Eq. (10)

###  4.1 Units

With exception of the Material class and the cim solver, in the nanobem
toolbox lengths are measured in nanometers and frequencies are converted to
the corresponding wavenumbers of light in vacuum

$k_{0}=\frac{2\pi}{\lambda}=\frac{\omega}{c}\,.$ (21)

$\lambda$ is the wavelength of light in vacuum, $\omega$ the angular
frequency, and $c$ the speed of light. Occasionally one would like to convert
from eV to wavenumbers. This can be done through

ene = 2; % photon energy in eVunits; % load conversion factor from eV to
wavelength in nmlambda = eV2nm / ene; % convert to wavelength in nmk0 = 2 * pi
/ lambda; % convert wavelength to wavenumberElectric fields are given in units
of the field strength of the incoming field, which can be also set to one,
such that the simulations give the field enhancement only. Instead of the
magnetic field we internally use $Z_{0}H$, with $Z_{0}$ being the free space
impedance, such that the electric field and the scaled magnetic field have the
same units.

###  4.2 Material class

The material properties of a medium are implemented in the nanobem toolbox
through the Material class

mat = Material( eps, mu ); % initialize Material objectHere eps is either a
constant or a permittivity function `eps(ene)` with the photon energy in eV as
input, with a corresponding functionality for the permeability `mu`. Both
`eps` and `mu` should return dimensionless numbers. We provide functions for
tabulated and Drude dielectric functions

eps = epstable( ’gold.dat’ ); % use tabulated dielectric functioneps =
epsdrude( ’Au’ ); % Drude dielectric functionOnce the material object is
initialized, we can compute for a given free space wavenumber `k0` the
permittivity `mat.eps(k0)`, the permeability `mat.mu(k0)`, the effective
wavenumber `mat.k(k0)`, and the impedance `mat.Z(k0)`. In all nanobem
simulations we first have to set up a table of material properties, e.g., for
air and silver

mat1 = Material( 1, 1 ); % airmat2 = Material( epstable( ’silver.dat’ ), 1 );
% silvermat = [ mat1, mat2 ]; % unique material vectorThe material vector mat
should be specified once, and then the same vector should be passed to all
objects and functions of the toolbox. We recommend setting the first entry of
the material vector to the material of the background medium.

###  4.3 Boundary

In a typical bem simulation we deal with several connected or unconnected
boundaries

$\partial\Omega=\bigcup_{j=1}^{n}\partial\Omega_{j}\,,$ (22)

where each boundary $\partial\Omega_{j}$ has an orientation, this is an inside
and outside specified with respect to the outer surface normal $\hat{\bm{n}}$,
and each boundary $\partial\Omega_{j}$ must have a single material at the
inside and outside. The boundaries $\partial\Omega_{j}$ are discretized using
triangles $\tau_{i}$ through

$\partial\Omega_{j}=\bigcup_{i=1}^{m}\tau_{i}\,,$ (23)

where $m$ is the total number of triangles. In the nanobem toolbox these
boundaries are specified using the `particle` class

p = particle( verts, faces ); % set up particle boundaryHere `verts` are the
triangle corners or vertices, and `faces` is a $m\times 3$ matrix describing
how these vertices are connected within the individual triangles. Note that
this implementation is in complete analogy to the face-vertex structures
provided by matlab. The toolbox provides a few implementations for simple
shapes, such as spheres, rods, or cubes, together with a plotting function for
the particle boundary. Details can be found in the help pages of the toolbox.
For more complicated setups the user must provide a discretization of the
boundaries of the desired structure in terms of vertices and triangular faces,
which often constitutes the most difficult part in setting up a bem
simulation.

###  4.4 Shape elements

Once we have set up the material vector and the boundaries, we can initialize
the Raviart-Thomas shape elements which play a central role in our bem
approach,

tau = BoundaryEdge( mat, p, inout );tau = BoundaryEdge( mat, p1, inout1, p2,
inout2, ... );Here `mat` is a unique material vector, `p` is a particle object
characterizing the discretized particle boundary, and `inout` is a vector
`[i1,i2]` where `mat(i1)` is the material at the particle inside and `mat(i2)`
the material at the particle outside. See also Fig. 5. For example, to set up
a silver sphere embedded in glass (refractive index 1.5) we proceed as follows

mat = [ Material( 1.5 ^ 2, 1 ), Material( epstable( ’silver.dat’ ), 1 )
];diameter = 50;p = trisphere( 144, diameter );tau = BoundaryEdge( mat, p, [
2, 1 ] );In the first line we set up the material vector, in the third line we
define a sphere triangulation with 144 vertices (diameter 50 nm), using the
toolbox function `trisphere`, and in the last line we set up the shape
elements for the nanosphere boundary with silver at the inside and glass at
the outside. `tau` is a BoundaryEdge vector specifying the shape elements for
the discretized particle boundary, where each vector element has the following
properties

tau(i).nu; % global degrees of freedom for triangle edgestau(i).val; %
auxiliary information for Thomas-Raviart basis elementstau(i).inout; % index
to material at inside and outsidetau(i).verts; % triangle vertices`nu` are the
three global degrees of freedom for the shape elements assigned to the three
edges of the triangle, see Eq. (10) and Fig. 5, and `val` is a value for the
Raviart-Thomas elements that is chosen such that the flow of $u_{E,H}$ in and
out of the triangle edges is conserved. `inout` specifies the material at the
inside and outside of the triangle boundary, and `verts` are the triangle
corners. Note that in our present implementation we internally assign the
shape elements to the boundary elements, and evaluate the single and double
layer potentials for triangle pairs rather than shape element pairs, which
extend over two triangles. This allows for some speedup for discretizations
with a moderate number of boundary elements, say a few thousands, but may
become inefficient for larger discretizations. From here on all information
about the discretized boundary is stored in the `BoundaryEdge` vector, the
information of the `particle` object is no longer needed. We note in passing
that the toolbox also implements linear shape elements `BoundaryVert` for
simulations using the quasistatic approach, for details see the help pages of
the toolbox.

###  4.5 Plane wave and dipole excitation

At present, the nanobem toolbox provides two types of external excitations,
namely plane waves and oscillating dipoles. Other excitations can be easily
implemented, as described in more detail in B, all a user must provide is a
structure of the form

qinc = struct( ’e’, e, ’h’, h, ’tau’, tau, ’k0’, k0 );Here `e`, `h` are the
inhomogeneities given in Eq. (13), `tau` is a vector for the boundary
elements, and `k0` is the free space wavenumber of the incoming fields. The
dimensions of `e`, `h` must be $n\times 1$ for a single excitation, where $n$
is the total number of edges or degrees of freedom, see Eq. (10). It is also
possible to specify multiple excitations, e.g. different polarization vectors,
by providing array dimensions of $n\times m$ or any other multidimensional
shape with leading dimension $n$. For plane wave excitations, we first set up
a `galerkin.planewave` object which is then evaluated for a given bounday
`tau` and wavenumber `k0`

exc = galerkin.planewave( [ 1, 0, 0; 0, 1, 0 ], [ 0, 0, 1; 0, 0, 1 ] );qinc =
exc( tau, k0 );The first array in the call to `galerkin.planewave` gives the
polarization vectors, here polarization along $x$ and $y$, and the second
array gives the light propagation direction, here along $z$. To set up an
excitation of oscillating dipoles, we must first specify the positions of the
dipoles and in which medium they are located

pt = Point( mat, imat, pos );`mat` is the same material vector that is passed
to `BoundaryEdge`, `imat` specifies the medium (scalar) or media (vector) in
which the dipole positions `pos` are located. Alternatively, the positions can
be also placed automatically in the photonic environment specified by `tau`

pt = Point( tau, pos );Once the `Point` object is initialized, we can set up
and evaluate the `galerkin.dipole` object in the same way as for plane waves

dip = galerkin.dipole( pt );qinc = dip( tau, k0 );`e`, `h` are arrays of
dimension $n\times m_{\rm pt}\times 3$, where $m_{\rm pt}$ is the number of
dipole positions and the last dimension accounts for the dipole orientations
along $x$, $y$, and $z$.

###  4.6 The boundary element method solver

To solve the bem equations (12), one has to set up a `galerkin.bemsolver`
object and invert the matrix equation for the excitation structure `qinc`

bem = galerkin.bemsolver( tau );sol = bem \ qinc;Both the initialization and
the solution of the bem solver are typically the most time consuming parts of
a simulation. The performance and the working principle of the bem solver can
be controlled through additional arguments passed in the initialization call
to the solver, as described in more detail in B.3 and the help pages of the
toolbox. The bem solution `sol` is an object of type

sol = galerkin.solution( tau, k0, e, h );where `tau` is the discretized
boundary, `k0` is the free space wavenumber for the external excitation, and
`e`, `h` hold the values
$\normalsize\llbracket{u_{E,H}^{e}}\normalsize\rrbracket$ characterizing the
tangential electromagnetic fields of Eq. (10). The bem solution object `sol`
can be used for various kinds of post processing. For a planewave excitation
we can pass the solution to the planewave object in order to compute the
optical cross sections

csca = scattering( exc, sol ); % scattering cross section in (nm^2)cabs =
absorption( exc, sol ); % absorption cross section in (nm^2)cext = extinction(
exc, sol ); % extinction cross section in (nm^2)For a dipole excitation, we
can use the dipole object in order to compute the total and radiative decay
rates in units of the free-space decay rate [12, chapter 10]

[ tot, rad ] = decayrate( dip, sol );For all types of excitations we can
compute the electromagnetic fields at the particle boundary and the surface
charge distribution

[ e1, h1 ] = interp( sol, ’inout’, 1 ); % fields at boundary inside[ e2, h2 ]
= interp( sol, ’inout’, 2 ); % fields at boundary outsidesig = surfc( sol ); %
surface charge distributionThe fields and the surface charge distribution are
given at the centroids of the boundary elements. The toolbox provides a number
of plotting routines for the fields and the surface charge distribution, as
described in more detail in the help pages. In order to evaluate the fields
away from the boundary, we can use the representation formulas (5). To this
end we set up a `Point` object and pass it together with the solution vector
to the `fields` function

[ x, y ] = ndgrid( len * linspace( -0.5, 0.5, n ) );pt = Point( tau, [ x( : ),
y( : ), 0 * x( : ) ] );[ e, h ] = fields( sol, pt ); % evaluate scattered
fields[ ei, hi ] = fields( exc, pt, k0 ); % evaluate incoming fieldsIn the
first line we set up a regular grid of size `len` with $\mbox{\tt
n}\times\mbox{\tt n}$ points centered around the origin, and place in the
second line the points into the photonic environment of `tau`. In the third
and fourth line we compute the scattered and incoming fields, which have the
dimension of $\mbox{\tt npt}\times 3\times\mbox{\tt m}$, where `npt` is the
number of points and `m` the number of excitations. Finally, we can also
compute the optical far fields through

[ efar, hfar ] = farfields( sol, pt );The points `pt` should be placed in the
embedding medium and give the direction (measured with respect to the origin)
along which the fields propagate far away from the nanoparticle.

###  4.7 The contour integral method solver

The contour integral method of Beyn [24] allows to determine the resonance
frequencies $\omega_{k}$ and modes $v_{k}$ of Eq. (15) that are located inside
a given closed contour $\Gamma$ in the complex frequency plane (see also Ref.
[31] for a related scheme). We start by briefly summarizing the method of Beyn
[24] for the determination of the modes fulfilling $A(\omega_{k})v_{k}=0$.
First, we set up a random matrix $V_{r}$ of dimension $n\times n_{r}$, where
$n$ is the matrix size of $A(\omega)$ (this is the total number of global
degrees of freedom for representing the electric and magnetic fields) and
$n_{r}$ is a value that has to be larger than the number of resonance modes in
the contour (which can only be determined at the end of the algorithm). We
usually set $n_{r}$ to a value of a few hundreds. In a second step, we compute
the matrices

$A_{j}=\frac{1}{2\pi i}\oint_{\Gamma}z^{j}A^{-1}(z)V_{r}\,dz$ (24)

for different values of $j$. We use an elliptic contour $\Gamma$ and
approximate the contour integration through a simple trapezoidal quadrature.
The solution of $A(z)u=V_{r}$ for $u$ in Eq. (24) requires the solution of the
bem equation, which is typically the time consuming part of cim (together with
the computation of the normalization constants to be described at the end). In
the following we describe the procedure for obtaining the resonance modes from
the knowledge of the matrices $A_{0}$, $A_{1}$, but we have also implemented
the algorithm using higher-order matrices along the lines discussed in the
paper of Beyn. Quite generally, in all our computer experiments the cim
approach worked perfectly when using $A_{0}$, $A_{1}$ only. In a third step,
we perform a singular value decomposition of
$A_{0}=V_{0}\Sigma_{0}W_{0}^{\dagger}$ and keep only singular values that are
larger than some cutoff value. In the toolbox, this step is implemented
through the graphical user interface described in Sec. 2, and we will discuss
this step in more detail below in Sec. 5. In the fourth and last step, we
compute the resonance frequencies $\omega_{k}$ and modes $v_{k}$ from the
solutions of the eigenvalue problem

$\left(V_{0}^{\dagger}A_{1}W_{0}\Sigma_{0}^{-1}\right)s_{k}=\omega_{k}s_{k}\,,\quad
v_{k}=V_{0}s_{k}\,.$ (25)

Once the resonance modes have been computed, we can check the accuracy of the
numerical scheme by computing the residuals

$r_{k}=\left\|A(\omega_{k})v_{k}\right\|\,.$ (26)

Finally, we compute the left modes $\tilde{v}_{k}$ from Eq. (16) and the
normalization factors in Eq. (20). Inside the toolbox these normalization
factors are absorbed into $\tilde{v}_{k}$. The resonance modes and the contour
in the complex frequency plane are given in units of the photon energy
$\hbar\omega$ in eV. The cim solver is set up with

cim = galerkin.cimsolver( bem ); % use previously initialized BEM solvercim =
galerkin.cimsolver( tau ); % or set up CIM solver using boundary elementsOnce
the cim solver is initialized, we can set the size $n_{r}$ for the random
matrix, the order `nz` for cim, and the contour $\Gamma$ through

cim.nr = 200; % size of random matrixcim.nz = 1; % order for CIM, default
value is 1cim.contour = cimbase.ellipse( [ zmin, zmax ], irad, nt );In the
initialization of the contour, `zmin` and `zmax` are the minimum and maximum
of the contour on the real axis, `irad` is the extension of the contour along
the imagainary axis, and `nt` is the number of quadrature points along the
contour. If needed, users can also shift the ellipse into the complex plane by
adding an imaginary part to the first argument `[zmin,zmax]+1i*zshift`. We
next call

data = eval1( cim, PropertyPairs );In this call the bem solver is initialized
if we have previously set up the cim solver with the boundary elements, where
the performance of the bem solver can be controlled with the same property
pairs as in the initialization of the bem solver, see also B, and additionally
the matrices $A_{j}$ of Eq. (24) are evaluated along the discretized contour.
This evaluation can be rather time consuming depending on the number of
boundary elements `tau` and quadrature points. The cutoff value for the
singular values can be set interactively by the user with the `tolselect`
function described in Sec. 2,

tol = tolselect( cim, data );if isempty( tol ), return; endcim = eval2( cim,
data, ’tol’, tol );The second line is needed for the case that the user stops
the simulation in the call to `tolselect`. In the call to `eval2` we compute
the left resonance modes and normalization constants, which can be again time
consuming because of the evaluation of $A^{\prime}(\omega_{k})$ for the
normalization constant. The `cim` object now stores the following properties

cim.nev % number of resonance modescim.ene % complex resonance
frequenciescim.sol % right resonance modescim.sol2 % left resonance
modescim.res % residualsFrom here on the cim solver can be used in the same
manner as the bem solver. For instance, for given inhomogeneities `qinc`, such
as those of a planewave or dipole excitation, the approximate solution $u_{\rm
CIM}$ of Eq. (20) can be obtained through

sol = cim \ qinc;The outcome of this call is a `galerkin.solution` vector that
can be used for post processing purposes in exactly the same way as previously
described in Sec. 3.2.

For larger nanoparticles one occasionally observes a consistent deviation
between the bem and cim results, see for instance the extinction spectra in
Fig. 3(a,b) or Fig. 6(b,c). As discussed for instance in [32, 33, 34], a
spectrally broad, featureless and non-resonant background can account for this
discrepancy. We will not enter into this topic here, but present a simple
approach for considering such a non-resonant background. In the spirit of Eq.
(20), we associate the difference between the inverse of the bem transmission
matrix $A^{-1}(\omega)$ and its cim approximate as the holomorphic matrix

$H(\omega)=A^{-1}(\omega)-\sum_{k=1}^{K}\Big{(}\tilde{v}_{k}^{T}A^{\prime}(\omega_{k})v_{k}\Big{)}^{-1}\frac{v_{k}\tilde{v}_{k}^{T}}{\omega-\omega_{k}}\,,$
(27)

which accounts for the background. As $H(\omega)$ is expected to have a weak
frequency dependence, it usually suffices to compute $H(\omega)$ at a few
sampling points only and to interpolate the function in between using a
quadratic or cubic fit. For instance, in Fig. 3(c) the three sampling energies
are indicated by arrows in the bottom of the figure. The spectral contribution
of this background is shown by the dashed line, and by adding it to the cim
results one obtains perfect agreement with the full bem results (compare
circles and cross symbols). To add a non-resonant background in the nanobem
toolbox, one first provides the sampling wavenumbers `ktab` and then adds the
contributions of Eq. (27) to the cim solver via

% sampling wavenumbers for interpolation of nonresonant matrixktab = 2 * pi ./
[ 400, 600, 700 ];% add nonresonant background to CIM solver, use BEM solver
stored e.g. in data.bemcim.nonresonant = diffcalderon( cim, bem, ktab );Once
the non-resonant part has been added, the approximate cim solution can again
be obtained through `sol=cim\qinc`.

##  5 Frequently asked questions about resonance modes

Figure 6: Comparison of results obtained with the nanobem toolbox and the
bem++ software [35]. We consider gold nanoparticles (Drude dielectric function
taken from [8, 18]) embedded in glass (refractive index $n=1.5$), and use
identical boundary discretizations and cim contours in the nanobem and bem++
simulations. The black lines and the cross symbols ($+$) give the results of
the bem and cim simulations obtained with the nanobem toolbox, the circle
symbols ($\circ$) give the results of the cim simulations obtained with the
bem++ software. The different panels show results for a nanosphere with a
diameter of (a) 50 nm and (b) 100 nm (305 vertices), and a nanocube with a
side length of (c) 50 nm (153 vertices) and (d) 100 nm (611 vertices). The
upper parts of the panels show the extinction cross sections $C_{\rm ext}$,
the lower parts display the contour for the determination of the resonance
modes and the resonance energies $\hbar\omega_{k}$. The nanobem and bem++
results are in perfect agreement throughout. The slight but systematic
disagreement between the bem and cim results for the larger structures in
panels (b,d) is due to the approximate nature of cim, as discussed in more
detail in [18], which can be approximately accounted for through a non-
resonant background (see Sec. 4.7). Figure 7: Convergence properties and
computer time. (a) Relative error between resonance energies for dipole and
quadrupole modes as computed within cim and Mie theory for the gold nanosphere
with a diameter of 100 nm investigated in Fig. 6(b) and for different boundary
discretizations. $h$ is the mean edge length of the triangulated sphere
boundary, for better visibility we use a quadratic scale. We use sphere
discretizations with 100, 200, 400, 600, 800, 1000, and 2000 vertices. (b)
Mean computer times for the solution of the bem equations at a single
frequency for the different sphere discretizations withd an increasing number
of global degrees of freedom, as obtained with a standard Intel Core processor
i7-7700 CPU at 3.60GHz. Typical cim simulations take from less than a minute
for coarse discretizations to several hours for fine discretizations.

The main purpose of this paper is to provide a toolbox for the simple and
efficient computation of resonance modes using the bem approach. Users
interested in such an approach are advised to set up their own problems, and
to critically examine them using computer experiments. From our own experience
we found that several things can go wrong in the computation of resonance
modes. In the following we list a few dos and don’ts.

Accumulation points. For nanostructures there usually exist accumulation
points where an infinite number of resonance frequencies converges at a single
frequency, e.g. the edge modes of a cube. The cim algorithm may fail when such
an accumulation point falls into the contour $\Gamma$. In such situations one
often finds that all singular values have the same order of magnitude, rather
than the standard situation where the large singular values are separated from
the smaller ones through a noticeable gap (see Fig. 4).

Accuracy for boundary integration. In particular for boundary discretizations
where the boundary elements have different sizes, it is important to select
high accuracy for triangle integrations. When the results of bem and cim
simulations differ, it might be a good idea to increase the accuracy for
boundary integrations, see table 3 for the control of the quadrature
performance using an option structure.

Accuracy of toolbox. We have tried our best to make the implementation of the
bem approach within the nanobem toolbox as transparent and accurate as
possible, where most of the quadrature rules can be directly controlled by the
user. Comparison of our toolbox results with Mie theory shows excellent
agreement throughout, demo programs can be obtained in the help pages. Fig. 6
shows for selected geometries a comparison of our results with those obtained
with the independent bem++ software [35], giving excellent agreement
throughout.

Analytic continuation. While a Drude dielectric function of the form

$\varepsilon(\omega)=1-\frac{\omega_{p}^{2}}{\omega(\omega+i\gamma)}$ (28)

can be easily extended to the complex frequency plane, a function tabulated
for real frequencies (e.g. `epstable`) cannot. This doesn’t necessarily mean
that there will be an error message in the computation of the resonance modes,
but the results will be meaningless. For this reason, users must assure that
their permittivity and permeability functions are suitable for a continuation
into the complex plane. For instance, both functions should be causal and thus
fulfill the Kramers-Kronig relation.

BEM and CIM results differ. The resonance mode expansion via a finite set of
modes is only approximate, and it may happen that the bem and cim results
differ systematically, for example if not all relevant modes are considered in
the cim. Also for larger structures a systematic discrepancy between bem and
cim might show up, see for instance panels (b,d) in Fig. 6. See also the
discussion about a non-resonant background given below.

Branch cuts. The contour $\Gamma$ of Eq. (24) should avoid crossing branch
cuts of the square root of the permittivity and permeability functions. For
Drude dielectric functions similar to Eq. (28) the toolbox has implemented a
suitable `zsqrt` function that avoids crossing branch cuts, for other material
functions one should be at least aware of this problem.

Convergence. Quite generally, the bem results should converge to the exact
solutions of the transmission problem of Eq. (1) for sufficiently small mesh
sizes $h$ and for sufficiently regular solutions. Beyond a certain threshold
value, which is hard to determine for systems where no analytic solutions are
available, the error due to the boundary discretization scales with
$h^{\nicefrac{{3}}{{2}}}$. As a rule of thumb, the ratio between $h$ and the
effective wavelength in the medium should be at least of the order of
$1\\!:\\!5$. In any case, convergence is a critical issue and should be
considered carefully in all simulations. Fig. 7 shows the convergence
properties for selected resonance modes of a sphere towards the exact Mie
results, as well as typical simulation runtimes.

Contour. We recommend to always use an elliptic contour for the cim approach,
as it can be shown to give numerically accurate results.

Cutoff for singular values. Singular values with a small magnitude might be
due to resonances that lie outside the contour, resonances that lie inside the
contour but cannot be resolved well enough due to a too small number of
quadrature points along the contour, and numerical artifacts. For a viable
choice of the cutoff parameter, the results of the reconstruction problem
should not not depend decisevely on its precise value because modes with a
small singular value contribute only little to the expansion of Eq. (20). For
this reason, the choice of this parameter is usually less critical than one
could expect.

Non-resonant background. For larger nanoparticles one occasionally observes a
consistent difference between the bem and cim results, which can be attributed
to a non-resonant background whose origin is not fully understood [32, 33,
34]. For a brief discussion and a possibility to account for such a background
within the nanobem toolbox see Sec. 4.7.

Number of resonance modes. For the optical spectra it usually suffices to keep
only a small number of modes in the simulation. For excitations of oscillating
dipoles, see e.g. democim02.m, the number of simulation modes should be
larger, in particular for dipoles situated close to the nanoparticle.

Residuals. To check the accuracy of the cim approach, it might be a good idea
to inspect the residuals for the resonance modes. See Sec. 4.7 for details.

Resonance frequency tracking. When bem and cim simulations agree for small
structures but differ for larger structures, it might be a good idea to track
the modes when continuously increasing the system size. It can happen that
specific resonances move out of the contour in the complex plane and are no
longer detected by the cim algorithm. In such cases, the contour should be
enlarged to encompass all relevant resonance frequencies.

Spectra versus nearfields. It often happens that the optical spectra computed
with bem and cim are very similar, but the optical nearfields differ. This is
no surprise because the spectra are only governed by the farfields features of
the resonance modes. For this reason, it may be beneficial to compute the
scattering spectra from the optical far fields, see for instance
`farscattering` command in democim01.m.

## Acknowledgment

We thank Marko Šimić for valuable feedback on the toolbox and Thomas Weiss for
a careful reading of the manuscript and most helpful suggestions. This work
has been supported in part by the Austrian Science Fund FWF under project P
31264 and by NAWI Graz.

##  Appendix A Representation formula

In this appendix we derive the representation formula (5) and its
implementation within the boundary element method approach. The derivation
closely follows the original work of Stratton and Chu [36], which is discussed
in more detail in [15, 12], and is repeated here for the sake of completeness.
Our starting point is the wave equation for the electric field

$-\nabla\times\nabla\times\bm{E}(\bm{r})+k^{2}\bm{E}(\bm{r})=-i\mu\omega\bm{J}(\bm{r})\,,$
(29)

where $k$ is the wavenumber inside a given medium and $\bm{J}(\bm{r})$ is an
external current distribution. It turns out to be convenient introducing a
Green’s function defined through

$-\nabla\times\nabla\times\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{r}^{\prime})+k^{2}\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{r}^{\prime})=-\delta(\bm{r}-\bm{r}^{\prime})\overset{\leftrightarrow}{\bm{I}},,$
(30)

which gives the electromagnetic response to a unit source located at position
$\bm{r}^{\prime}$ subject to the boundary condition of outgoing waves at
infinity (Silver-Müller radiation condition). It can be shown that the Green’s
function has the form of an outgoing spherical wave [15, 12]

$\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{r}^{\prime})=\left(\overset{\leftrightarrow}{\bm{I}}+\frac{\nabla\nabla}{k^{2}}\right)g(\bm{r},\bm{r}^{\prime})\,,\qquad
g(\bm{r},\bm{r}^{\prime})=\frac{e^{ik|\bm{r}-\bm{r}^{\prime}|}}{4\pi|\bm{r}-\bm{r}^{\prime}|}\,,$
(31)

where we have used the dyadic product
$(\nabla\nabla)_{ij}=\partial_{i}\partial_{j}$ and
$\overset{\leftrightarrow}{\bm{I}}$ is the unit matrix. The matrix
$\overset{\leftrightarrow}{\bm{G}}$ is often referred to as the dyadic Green’s
function or Green’s dyadics, in short. Apart from a prefactor,
$G_{ij}(\bm{r},\bm{r}^{\prime})$ relates a current source at position
$\bm{r}^{\prime}$ oriented along the Cartesian direction $j$ to the $i$’th
component of the electric field at position $\bm{r}$. We next multiply the
wave equation (29) for the electric field from the right-hand side with the
Green’s dyadics

$\Big{[}-\nabla^{\prime}\times\nabla^{\prime}\times\bm{E}(\bm{r}^{\prime})\Big{]}\cdot\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})+k^{2}\,\bm{E}(\bm{r}^{\prime})\cdot\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})=-i\omega\mu\bm{J}(\bm{r}^{\prime})\cdot\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})\,.$

For later convenience we have written the wave equation for the primed
positions $\bm{r}^{\prime}$ rather than the unprimed ones. Similarly, we
multiply Eq. (30) from the left-hand side with $\bm{E}(\bm{r}^{\prime})$,

$-\bm{E}(\bm{r}^{\prime})\cdot\nabla^{\prime}\times\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})+k^{2}\,\bm{E}(\bm{r}^{\prime})\cdot\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})=-\bm{E}(\bm{r}^{\prime})\delta(\bm{r}^{\prime}-\bm{r})\,.$

We next integrate the two equations over $\bm{r}^{\prime}$, which is assumed
to be located in volume $\Omega$, and subtract them to arrive at [15, Eq.
(8.1.18)]

$\int_{\Omega}\left\\{\Big{[}-\nabla^{\prime}\times\nabla^{\prime}\times\bm{E}(\bm{r}^{\prime})\Big{]}\cdot\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})+\bm{E}(\bm{r}^{\prime})\cdot\nabla^{\prime}\times\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})\right\\}d^{3}r^{\prime}=-\bm{E}_{\rm
inc}(\bm{r})+\bm{E}(\bm{r})\,,$ (32)

with the incoming electric field

$\bm{E}_{\rm
inc}(\bm{r})=i\mu\omega\int_{\Omega}\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{r}^{\prime})\cdot\bm{J}(\bm{r}^{\prime})\,d^{3}r^{\prime}$
(33)

that would be the solution of Maxwell’s equations for an unbounded medium in
absence of the boundary $\partial\Omega$. The term in curly brackets of Eq.
(32) can be rewritten as [15, Eq. (8.1.19)]

$-\nabla^{\prime}\cdot\Bigl{\\{}[\nabla^{\prime}\times\bm{E}(\bm{r}^{\prime})]\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})+\bm{E}(\bm{r}^{\prime})\times[\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})]\Bigr{\\}}\,.$
(34)

We next use this modified expression together with Gauss’ theorem to obtain a
boundary integral of the form

$\bm{E}(\bm{r})=\bm{E}_{\rm
inc}(\bm{r})-\oint_{\partial\Omega}\Bigl{\\{}[\nabla^{\prime}\times\bm{E}(\bm{r}^{\prime})]\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})+\bm{E}(\bm{r}^{\prime})\times\bigl{[}\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r}^{\prime},\bm{r})\bigr{]}\Bigr{\\}}\cdot\hat{\bm{n}}^{\prime}dS^{\prime}\,,$
(35)

where we have explicitly indicated the surface normal $\hat{\bm{n}}^{\prime}$
of the boundary. The first term in curly brackets can be simplified by
performing a cyclic permutation of the triple product and using Faraday’s law
for the curl of $\bm{E}$, and the second one by performing a cyclic
permutation of the triple product. This leads us to our final expression [12,
Eq. (5.26)]

$\bm{E}(\bm{r})=\bm{E}_{\rm
inc}(\bm{r})=-\oint_{\partial\Omega}\left\\{i\omega\mu\,\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{s}^{\prime})\cdot\hat{\bm{n}}^{\prime}\times\bm{H}(\bm{s}^{\prime})-\left[\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{s}^{\prime})\right]\cdot\hat{\bm{n}}^{\prime}\times\bm{E}(\bm{s}^{\prime})\right\\}\,dS^{\prime}\,.$
(36)

In deriving this expression we have exploited that
$\overset{\leftrightarrow}{\bm{G}}$ is a symmetric tensor. We next introduce
the single layer potential [12, Eq. (5.34)]

$\bigl{[}\mathbb{S}\,\bm{u}\bigr{]}(\bm{r})=\oint_{\partial\Omega}\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{s}^{\prime})\cdot\bm{u}(\bm{s}^{\prime})\,dS^{\prime}=\oint_{\partial\Omega}\left[g(\bm{r},\bm{s}^{\prime})\bm{u}(\bm{s}^{\prime})+\frac{1}{k^{2}}\nabla
g(\bm{r},\bm{s}^{\prime})\,\nabla^{\prime}\cdot\bm{u}(\bm{s}^{\prime})\right]\,dS^{\prime}\,,$
(37)

where we have used integration by parts to arrive at the second expression,
together with the double layer potential

$\bigl{[}\mathbb{D}\,\bm{u}\bigr{]}(\bm{r})=\oint_{\partial\Omega}\nabla^{\prime}\times\overset{\leftrightarrow}{\bm{G}}(\bm{r},\bm{s}^{\prime})\cdot\bm{u}(\bm{s}^{\prime})\,dS^{\prime}=-\oint_{\partial\Omega}\nabla\times
g(\bm{r},\bm{s}^{\prime})\bm{u}(\bm{s}^{\prime})\,dS^{\prime}\,$ (38)

With these potentials we can cast Eq. (36) to the representation formula (5)
given in the main text. The expressions outside of volume $\Omega$ can be
obtained in complete analogy to the above derivation but for the integration
over volume $\mathbb{R}^{3}\backslash\Omega$ and by observing that the
outgoing fields at infinity go to zero.

In the numerical implementation for the bem approach we have to evaluate the
single and double layer potentials for the Raviart-Thomas shape functions, see
Eq. (11), to arrive at the expressions

(39a) (39b)

where we have performed for $\mathbb{S}$ an integration by parts in order to
bring the derivative from the Green’s function to the Raviart-Thomas shape
element. The potentials of Eq. (39) are at the heart of our bem approach. More
details about the numerical evaluation will be presented in B.

##  Appendix B Details about internal toolbox functions

Table 2: List of the internal classes of the nanobem toolbox. The tensor class is described in B.1, the quadrature classes in B.2, and the quadrature engine in B.3. Class name |  Description
---|---
tensor |  Tensor class for simple manipulation of multidimensional arrays
quadboundary |  Quadrature rules for integration over boundary elemens
quadduffy |  Specialized quadrature rules for touching triangle pairs [37, 38]
potbase3 |  Auxiliary class for analytic integration over triangle pairs [39]
galerkin.pot1.base |  Base class for evaluation of single and double layer potentials, Eq. (12)
galerkin.pot1.std |  Standard evaluation over triangle pairs with low accuracy
galerkin.pot1.refine1 |  Refined integration over close triangle pairs w/o series expansion
galerkin.pot1.refine2 |  Refined integration over close triangle pairs with series expansion
galerkin.pot1.duffy1 |  Duffy-type integration over touching triangle pairs w/o series expansion
galerkin.pot1.duffy2 |  Duffy-type integration over touching triangle pairs with series expansion
galerkin.pot1.smooth1 |  Partially analytic integration over triangle pairs w/o series expansion
galerkin.pot1.smooth2 |  Partially analytic integration over triangle pairs with series expansion
galerkin.pot2.base |  Base class for evaluation of representation formula, Eq. (5)
galerkin.pot2.std |  Standard evaluation over triangle with low accuracy
galerkin.pot2.smooth |  Partially analytic integration over triangle [39]

In this appendix we provide information about the internal toolbox functions
and how their performance can be controlled through property pairs. Table 2
provides a list for the most important toolbox classes.

###  B.1 Tensor class

The nanobem toolbox provides a simple tensor class for the evaluation of
expressions like

$\quad a_{ij}b_{j}\,,\quad\sum_{j}a_{ij}b_{j}\,,\quad
A_{ikj}b_{ij}\,,\dots\,.$

The tensor class avoids the time consuming `for` loops of matlab and uses the
faster `bsxfun` command instead. We first introduce dummy expressions for the
indices

[ i, j, k ] = deal( 1, 2, 3 );Inside the toolbox functions a value of 1 then
represents index $i$, a value of 2 represents index $j$, and so on, with the
only important restriction that all indices must have different values. Tensor
objects are initialized via

a = tensor( rand( 3, 4 ), [ i, j ] );b = tensor( rand( 4, 1 ), j );A = tensor(
rand( 3, 5, 4 ), [ i, k, j ] );We here use random arrays for demonstration
purposes only. Once the tensor objects have been initialized, they can be
manipulated through

c1 = a * b;c2 = sum( a * b, j );c3 = A * b;Tensor objects can be manipulated
with other tensor objects and scalars, but not with regular matlab arrays. The
usual operators such as `+`, `-`, `*`, `./`, `sin`, `cos`, `exp` have been
overloaded for the tensor class and make use of the builtin `bsxfun` command
of matlab. The tensor class is particularly useful for clear and compact
coding, when speed is an issue the use of tensor objects should be examined
critically. Note that internally the consistency of the dimensions for the
various tensor arrays is not checked. Finally, conversion of the tensor object
back to a normal matlab array is done through

c1 = double( c1, [ i, j ] );c2 = double( c2, i );c3 = double( c3, [ i, j, k ]
);As a representative example we consider the case where the electromagnetic
fields for a planewave excitation are evaluated at the positions `pos` on the
boundary of the nanoparticle. The light polarizations `pol` and directions
`dir` are arrays of dimension $\mbox{\tt npol}\times 3$, and `pos` is an array
of dimension $\mbox{\tt ntau}\times\mbox{\tt nq}\times 3$, where `ntau` and
`nq` are the number of boundary elements and quadrature points per element,
respectively. We also provide the material object `mat` for the background
medium and the wavenumber `k0`. The electromagnetic fields can then be
computed with the following function:

function [ e, h ] = fields( pol, dir, mat, k0, pos ) [ i, q, ipol, k ] = deal(
1, 2, 3, 4 ); % dummy indices for tensor class pol = tensor( pol, [ ipol, k ]
); % convert polarizations to tensor class dir = tensor( dir, [ ipol, k ] ); %
convert directions to tensor class pos = tensor( pos, [ i, q, k ] ); % convert
positions to tensor class e = exp( 1i * mat.k( k0 ) * dot( pos, dir, k ) ) *
pol; h = mat.Z( k0 ) * cross( dir, e, k ); e = double( e, [ i, q, k, ipol ] );
% convert E to numeric array h = double( h, [ i, q, k, ipol ] ); % convert H
to numeric arrayend

###  B.2 Boundary integration

At present the toolbox can only deal with triangular boundary elements, but it
is written such that in the future also quadrilateral elements might be
considered. For triangle integrations, we first set up a structure

quad3 = struct( ’x’, x, ’y’, y, ’w’, w );Here, `x`, `y` are the quadrature
points for a unit triangle, and `w` are the integration weights. A number of
triangle quadrature rules is provided through

quad3 = triquad( rule ); % integration rules 1-19quad3 = triquad( rule,
’refine’, 2 ); % split triangleOnce the integration points and weights for the
unit triangle are specified, they can be used for the initialization of the
quadrature rules for boundary integration

pt = quadboundary( tau, ’quad3’, quad3 );`pt` is an `quadboundary` object or
an array thereof, with the elements

pt.tau % boundary elements with same materials at inside and outsidept.quad %
quadrature rules for boundary elementspt.npoly % number of edges of boundary
element, 3 for trianglespt.mat % material vectorpt.inout % index [i1,i2] to
material vector, i1 for inside and i2 for outsideEach `quadboundary` object
always groups boundary elements `tau` with the same shape (currently only
triangles) and with the same materials at the inside and outside. This can
considerably simplify the evaluation of boundary integrals, as will be
demonstrated below. With

[ pos, w, f ] = eval( pt );[ pos, w, f, fp ] = eval( pt );we obtain the
quadrature positions of dimension $\mbox{\tt ntau}\times\mbox{\tt nq}\times
3$, where `ntau` is the number of boundary elements `pt.tau` and `nq` the
number of quadrature points per boundary element, the quadrature weights `w`
of dimension $\mbox{\tt ntau}\times\mbox{\tt nq}$, and the Raviart-Thomas
shape elements `f` of Eq. (10) with dimension $\mbox{\tt ntau}\times\mbox{\tt
nq}\times\mbox{\tt npoly}\times 3$, where `npoly` labels the three shape
elements associated with the triangle edges. The additional output argument
`fp` gives the divergence of the shape elements $\nabla\cdot\bm{f}_{\nu}^{e}$
of size $\mbox{\tt ntau}\times\mbox{\tt nq}\times\mbox{\tt npoly}$, which is
needed in some cases.

As a representative example, in the following we show how to compute the
inhomogeneities of Eq. (13) for a planewave excitation

$\llbracket{q_{E}^{\rm
inc}}\rrbracket_{ia}=\oint_{\tau_{i}}\bm{f}_{ia}^{e}(\bm{s})\cdot\bm{E}^{\rm
inc}(\bm{s})\,dS\,,$ (39bd)

with a similar expression for the magnetic field. Note that the shape function
$\bm{f}_{ia}^{e}$ is here defined for a specific triangle $\tau_{i}$ and for
edge $a$, the translation to the global degrees of freedom $\nu$ will be done
at the end. Using the field evaluation function of B.1 and a single
`quadboundary` object `pt`, the inhomogeneities can be computed through

function [ qe, qh ] = qfields( pol, dir, mat, k0, pt ) [ pos, w, f ] = eval(
pt ); [ e, h ] = fields( pol, dir, mat, k0, pos ); [ i, ipol, q, a, k ] =
deal( 1, 2, 3, 4, 5 ); e = tensor( e, [ i, q, k, ipol ] ); h = tensor( h, [ i,
q, k, ipol ] ); f = tensor( f, [ i, q, a, k ] ); w = tensor( w, [ i, q ] ); qe
= sum( dot( f, e, k ) * w, q ); qe = double( qe, [ i, a, ipol ] ); qh = sum(
dot( f, h, k ) * w, q ); qh = double( qh, [ i, a, ipol ] );endWe can finally
use this function for the total evaluation of the inhomogeneities of Eq. (13).
In a first step we allocate the arrays for the inhomogeneities, using the
number of global degrees of freedom `ndof(tau)` for the boundary elements. We
then loop over the boundary elements and evaluate the inhomogeneities of Eq.
(39bd) provided that the boundary is connected to the embedding medium.
Finally, we translate from the local degrees of freedom $i,a$ to the global
degrees $\nu$.

[ qinc.e, qinc.h ] = deal( zeros( ndof( tau ), size( pol, 1 ) ); %
allocationfor pt = quadboundary( tau, ’quad3’, triquad( 3 ) ) % loop over
boundary elements if pt.inout( 2 ) == 1 [ qe, qh ] = qfields( pol, dir, mat,
k0, pt ); % field evalulation qe = reshape( qe, [], size( pol, 1 ) ); qh =
reshape( qh, [], size( pol, 1 ) ); nu = vertcat( pt.tau.nu ); % translation to
global degrees of freedom for it = 1 : numel( nu ) qinc.e( nu( it ), : ) =
qinc.e( nu( it ), : ) + qe( it, : ); qinc.h( nu( it ), : ) = qinc.h( nu( it ),
: ) + qh( it, : ); end endend

###  B.3 Quadrature engine

For the evaluation of the single and double layer potentials we have to
evaluate integrals of the form [12]

$\mathcal{I}_{ia,i^{\prime}a^{\prime}}=\oint_{\tau_{i}}\oint_{\tau_{i^{\prime}}}\bm{f}_{ia}^{e}(\bm{s})\cdot\bm{f}_{i^{\prime}a^{\prime}}^{e}(\bm{s}^{\prime})\frac{e^{ik_{j}|\bm{s}-\bm{s}^{\prime}|}}{4\pi|\bm{s}-\bm{s}^{\prime}|}\,dSdS^{\prime}\,.$
(39be)

The integrand diverges for $\bm{s}\approx\bm{s}^{\prime}$ and one has to be
careful in the evaluation of the integral for identical or touching triangle
pairs, as well as for triangles that are located close to each other. The
evaluation of integrals of the form of Eq. (39be) is performed in the toolbox
through a quadrature engine, which allows to use different quadrature rules
(coarse or refined) for different triangle pairs. We also provide the
possibility that the exponential in Eq. (39be) is expanded in a Taylor series,
such that the time consuming evaluations of the double integrals have to be
performed only once. The order of the Taylor expansion is controlled by the
parameter `order` (with a default value of three). In the following we
describe the quadrature engine and the parameters for its control in slightly
more detail.

Figure 8: Quadrature engine for the evaluation of integrals over triangle
pairs $\mathcal{I}(\tau_{1},\tau_{2})$. In the first step std, we evaluate the
integrals using quadrature rules of low order with only few quadrature points
per triangle. In the second step refine, we use a refined quadrature rules for
triangle pairs that are sufficiently close to each other. In the final step
duffy, we use specialized quadrature rules for identical or touching triangle
pairs to deal with the singular integrands of the single and double layer
potentials.

We first introduce a measure for the distance between two triangles
$\tau_{1}$, $\tau_{2}$. First, we enclose each boundary element $\tau$ in a
sphere with radius $\mbox{rad}(\tau)$ that is centered around
$\mbox{pos}(\tau)$. The mean distance scaled by the bounding box radius is
then defined as the relative distance

$d(\tau_{1},\tau_{2})=\frac{\|\mbox{pos}(\tau_{1})-\mbox{pos}(\tau_{2})\|}{\mbox{rad}(\tau_{1})+\mbox{rad}(\tau_{2})}\,.$
(39bf)

When this distance is smaller than some parameter `relcutoff` we usually
perform a refined integration. The working principle of our quadrature engine
is depicted in Fig. 8. We start with a coarse integration for all triangle
pairs, and then successively employ more refined quadrature rules. In each
refinement step the previous results are overwritten. To set up the quadrature
engine, we first define the quadrature rules for the standard and refined
integration

rules = quadboundary.rules( ’quad3’, triquad( 3 ) ); % standardrules1 =
quadboundary.rules( ’quad3’, triquad( 11 ) ); % refinedThe quadrature engine
is then set up for instance through

pot( 1 ) = galerkin.pot1.std( ’rules’, rules );pot( 2 ) =
galerkin.pot1.refine2( ’relcutoff’, 2, ’order’, 3, ’rules1’, rules1 );pot( 3 )
= galerkin.pot1.duffy2( ’nduffy’, 3 );Here `order` controls the Taylor series
expansion, `relcutoff` controls the triangle distance for refined integration,
and `nduffy` controls the quadrature rules for touching or identical triangle
pairs. The quadrature engine can be passed to the bem solver through

bem = galerkin.bemsolver( tau, ’engine’, pot );

Table 3: The performance of the galerkin.bem and galerkin.cim solvers can be controlled through a number of options that can be passed in the form of property names and values to the solvers. Alternatively one can also use a option structure, in complete analogy to the inputParser object ofmatlab. Property name | Default value |  Description
---|---|---
’relcutoff’ | 2 |  Cutoff for refined triangle integrations
’rules’ | triquad(3) |  Default integration rules, see quadboundary.rules
’rules1’ | triquad(11) |  Quadrature rules for refined integration
’order’ | 3 |  Order of series expansion, [] if none
’nduffy’ | 3 |  Parameter controlling Duffy integration rules for touching triangles
’waitbar’ | 0 |  Show waitbar during initialization

In the evaluation of the single and double layer potentials we then progress
in `bem` through the `pot` array, where each element provides the quadrature
rules for the coarse or refined triangle integrations. Alternatively, we can
also set up an option array and pass the arguments controlling the engine
directly to the bem solver (see also Table 3)

op = struct( ’relcutoff’, 2, ’order’, 3, ... );bem = galerkin.bemsolver( tau,
op );Note that with `order=[]` no series expansion is performed in the
exponential of Eq. (39be), which can lead to more accurate results at the
expense of a significant increase of simulation time. A similar quadrature
engine approach is used when computing the fields away from the boundary,
using the representation formula (5). The relative distance between position
$\bm{r}$ and the boundary element $\tau$ is defined through the distance
between $\bm{r}$ and $\mbox{pos}(\tau)$ scaled by the bounding box radius
$\mbox{rad}(\tau)$. The quadrature engine is set up in the same way as before

pot( 1 ) = galerkin.pot2.std( ’rules’, rules );pot( 2 ) =
galerkin.pot2.smooth( ’relcutoff’, 2, ’rules’, rules1 );This engine is used
for instance in the `fields` function for the bem solution vector, where the
performance can be controlled again by passing arguments to the function. The
`smooth` integration function uses the singularity subtraction method of [39]
for the evaluation of the integrals with divergent integrands. More details
about the performance and functionality of the engine functions can be found
in the help pages of the toolbox.

## References

  * [1] P. T. Leung, S. Y. Liu, K. Young, Completeness and orthogonality of quasinormal modes in leaky optical cavities, Phys. Rev. A 49 (1994) 3057.
  * [2] P. T. Kristensen, C. Van Vlack, S. Hughes, Generalized effective mode volume for leaky optical cavities, Opt. Lett. 37 (2012) 1649.
  * [3] Q. Bai, M. Perrin, C. Sauvan, J.-P. Hugonin, P. Lalanne, Efficient and intuitive method for the analysis of light scattering by a resonant nanostructure, Opt. Express 21 (22) (2013) 27371–27382.
  * [4] C. Sauvan, J. P. Hugonin, I. S. Maksymov, P. Lalanne, Theory of the spontaneous optical emission of nanosize photonic and plasmon resonators, Phys. Rev. Lett. 110 (2013) 237401.
  * [5] R. C. Ge, P. T. Kristensen, J. F. Young, S. Hughes, Quasinormal mode approach to modelling light-emission and propagation in nanoplasmonics, New J. Phys. 16 (2014) 113048.
  * [6] M. B. Doost, W. Langbein, E. A. Muljarov, Resonant-state expansion applied to three-dimensional open optical systems, Phys. Rev. A 90 (2014) 013834.
  * [7] E. A. Muljarov, W. Langbein, Resonant-state expansion of dispersive open optical systems: Creating gold from sand, Phys. Rev. B 93 (2016) 075417.
  * [8] M. Perrin, Eigen-energy effects and non-orthogonality in the quasi-normal mode expansion of Maxwell equations, Opt. Express 24 (24) (2016) 27137–27151.
  * [9] P. Lalanne, W. Yan, A. Gras, C. Sauvan, J.-P. Hugonin, M. Besbes, G. Demesy, M. D. Truong, B. Gralak, F. Zolla, A. Nicole, F. Binkowski, L. Zschiedrich, S. Burger, J. Zimmerling, R. Remis, P. Urbach, H. T. Liu, T. Weiss, Quasinormal mode solvers for resonators with dispersive materials, J. Opt. Soc. Am. A 36 (2019) 686.
  * [10] P. T. Kristensen, K. Herrmann, F. Intravaia, K. Busch, Modeling electromagnetic resonators using quasinormal modes, Adv. Opt. Photon. 12 (2020) 612–708.
  * [11] J. A. Schuller, E. S. Barnard, W. Cai, Y. C. Jun, J. S. White, M. L. Brongersma, Plasmonics for extreme light concentration and manipulation, Nature Mat. 9 (2010) 193.
  * [12] U. Hohenester, Nano and Quantum Optics, Springer, 2020.
  * [13] J. S. Hesthaven, T. Warburton, High-order/spectral methods on unstructured grids I. Time-domain solution of Maxwell’s equations, J. Comput. Phys. 181 (2002) 186.
  * [14] J. S. Hesthaven, High-order accurate methods in time-domain computational electromagnetics: A review, Advances in Imaging and Electron Physics 127 (2003) 59.
  * [15] W. C. Chew, Waves and fields in inhomogeneous media, IEEE Press, Picsatoway, 1995\.
  * [16] A. M. Kern, O. J. F. Martin, Surface integral formulation for 3d simulations of plasmonics and high permittivity nanostructures, J. Opt. Soc. Am. A 26 (2009) 732\.
  * [17] S. Franke, S. Hughes, M. K. Dezfouli, P. T. Kristensen, K. Busch, A. Knorr, M. Richter, Quantization of quasinormal modes for open cavities and plasmonic cavity quantum electrodynamics, Phys. Rev. Lett. 122 (2019) 213901.
  * [18] G. Unger, A. Trügler, U. Hohenester, Novel modal approximation scheme for plasmonic transmission problems, Phys. Rev. Lett. 121 (2018) 24680.
  * [19] U. Hohenester, A. Trügler, MNPBEM - A Matlab Toolbox for the simulation of plasmonic nanoparticles, Comp. Phys. Commun. 183 (2012) 370.
  * [20] U. Hohenester, Simulating electron energy loss spectoscopy with the MNPBEM toolbox, Comp. Phys. Commun. 185 (2014) 1177.
  * [21] U. Hohenester, Making simulations with the MNPBEM toolbox big: Hierarchical matrices and iterative solvers, Comp. Phys. Commun. 222 (2018) 209.
  * [22] F. J. García de Abajo, A. Howie, Retarded field calculation of electron energy loss in inhomogeneous dielectrics, Phys. Rev. B 65 (2002) 115418.
  * [23] C. Sauvan, Quasinormal mode expansion for nanoresonators made of absorbing dielectric materials: study of the role of static modes, Opt. Express 29 (2021) 8268.
  * [24] W.-J-Beyn, An integral method for solving nonlinear eigenvalue problems, Linear Algebra and its Applications 436 (2012) 3839.
  * [25] Y. Chang, R. Harrington, A surface formulation for characteristic modes of material bodies, IEEE Transactions on Antennas and Propagation 25 (6) (1977) 789–795.
  * [26] A. J. Poggio, E. K. Miller, Chapter 4: Integral equation solutions of three-dimensional scattering problems, in: R. Mittra (Ed.), Computer Techniques for Electromagnetics, International Series of Monographs in Electrical Engineering, Pergamon, 1973, pp. 159 – 264\.
  * [27] T. K. Wu, L. L. Tsai, Scattering from arbitrarily-shaped lossy dielectric bodies of revolution, Radio Science 12 (5) (1977) 709–718.
  * [28] F. Zolla, A. Nicolet, G. Demésy, Photonics in highly dispersive media: the exact modal expansion, Opt. Lett. 43 (2018) 5813–5816.
  * [29] V. Kozlov, V. Mazčya, Differential Equations with Operator Coefficients with Applications to Boundary Value Problems for Partial Differential Equations, Springer, Berlin, 1999.
  * [30] R. Mennicken, M. Möller, Non-Self-Adjoint Boundary Eigenvalue Problems, North-Holland Publishing Co., Amsterdam, 2003.
  * [31] D. A. Bykov, L. L. Doskolovich, Numerical methods for calculating poles of the scattering matrix with applications in grating theory, Journal of Lightwave Technology 31 (2013) 793–801.
  * [32] L. Zschiedrich, F. Binkowski, N. Nikolay, O. Benson, G. Kewes, and S. Burger, Riesz-projection-based theory of light-matter interaction in dispersive nanoresonators, Phys. Rev. A 98 (2018) 043806.
  * [33] T. Weiss and E. A. Muljarov, How to calculate the pole expansion of the optical scattering matrix from the resonant states, Phys. Rev. B 98 (2018) 085433.
  * [34] R. Colom, R. McPhedran, B. Stout, and N. Bonod, Modal expansion of the scattered field: Causality, nondivergence, and nonresonant contribution, Phys. Rev. B 98 (2018) 085418.
  * [35] W. Śmigaj, S. Arridge, T. Betcke, J. Phillips, J. Schweiger, Solving Boundary Integral Problems with BEM++, ACM Trans. Math. Software 41 (6) (2015) 1–40.
  * [36] J. A. Stratton, L. J. Chu, Diffraction theory of electromagnetic waves, Phys. Rev. 56 (1939) 99–107.
  * [37] D. J. Taylor, Accurate and efficient numerical integration of weakly singulars integrals in Galerkin EFIE solutions, IEEE Trans. Antennas Propag. 51 (2003) 2543.
  * [38] S. Sarraf, E. Lopez, G. Rios Rodriguez, J. D’Elia, Validation of a Galerkin technique on a boundary integral equation for creeping flow around a torus, Comp. Appl. Math. 33 (2014) 63.
  * [39] I. Hänninen, M. Taskinen, J. Sarvas, Singularity subtraction integral formulae for surface integral equations with RWG, rooftop and hybrid basis functions, Progress in Electromagnetic Research, PIERS 63 (2006) 243.