Computing electromagnetic modes with scuffspectrum
scuffspectrum is a tool within the scuffem code suite for computing frequencies and field patterns of resonant modes in arbitrary material geometries. More specifically, scuffspectrum implements Beyn's contourintegral algorithm for nonlinear eigenproblems, described in more detail below, to pinpoint frequencies at which the BEM system matrix is singular, allowing nonzero surface currents to flow in the absence of external fields. scuffspectrum reports mode frequencies and, optionally, visualization diagrams and various types of information on the current and field patterns of the corresponding eigenvectors.
Of course, the information reported by scuffspectrum could also be obtained indirectly from other codes in the scuffem application suitefor example, the resonant modes of the spherical dielectric cavity shown below would show up as blips in (say) the scattering cross section as the frequency is scanned past the vicinity of the resonance. However, this would be an inefficient and imprecise way to determine resonance frequencieswe would probably wind up computing scattering cross sections at hundreds of frequencies on a dense grid just to approximate the resonance with one or twodigit accuracyand it would not yield direct information on the current and field patterns of the corresponding modes. In contrast, the method used by scuffspectrum is optimized to pinpoint mode frequencies with manydigit accuracyusing only a small number of samplesand, moreover, directly furnishes all information we may need about the current a and field patterns of the modes themselves.
1. Beyn's algorithm: the computational engine behind scuffspectrum
Mode frequencies as nonlinear eigenvalues
The basic problem solved by scuffspectrum is to find values of the angular frequency at which the BEM system matrix is singular (has nonempty nullspace), in which case we can find a nonzero vector that solves the equation
This is just the usual linear system of equations solved in scuffem scattering problems, but now with zero righthand side, i.e. there is no incident field or other external stimulus exciting the structure; instead, the eigenvector describes a configuration of surface currents that can exist in the absence of external fieldsa resonance mode.
Equation (1) resembles an eigenvalue problem, but with the complication that the matrix depends in a complicated nonlinear way on the frequency ; this means we can't simply use standard solvers for linear eigenproblems such as those implemented in lapack.
Beyn's contourintegral algorithm: A lightning overview
Instead, scuffspectrum implements the contourintegral approach to nonlinear eigenproblems proposed by W. Beyn in this 2012 paper:

WolfJürgen Beyn, "An integral method for solving nonlinear eigenvalue problems." Linear Algebra and its Applications 436 3839 (May 2012).
In a nutshell, Beyn's method locates eigenvalues of (1) by evaluating a certain contour integral in the complex plane; the evaluation proceeds via numerical quadrature, sampling the integrand at a total of quadrature points (see figure below). The number of quadrature points, and the shape and location of the contour, are usertweakable parameters that you will specify as inputs to scuffspectrum; the basic idea is that Beyn's method will identify all eigenvalues located inside the contour, and will ignore any eigenvalues lying outside the contour, so the challenge from the user's perspective is to have a rough idea of where your eigenvalues will beand to design contours that enclose the ones you want while omitting the others.
We are being deliberately vague here about precisely what contour integral we are evaluating, as this is not crucial knowledge for using scuffspectrum (for details, see the original paper linked above). Pretty much all you need to know is this: Evaluating the integrand at each quadrature point involves (a) forming and factorizing the matrix , and then (b) doing some simple linearalgebra calculations involving the solution of a small number (see below) of linear systems of equations of the form for certain given RHS vectors Of course, assembling the system matrix and solving linear systems is exactly what scuffem does to solve scattering problems, and so mechanically Beyn's algorithm boils down simply to solving integralequation scattering problems, at some number of frequencies, with some number of RHS vectors at each frequencynot different, in principle, from what we would wind up doing in the bruteforce (frequencyscanning) approach to computing mode frequencies discussed above. The difference is that Beyn's algorithm chooses the frequencies at which we calculate, and combines the results of all those calculations, in clever ways to ensure that maximal benefit is extracted from the computations, allowing us to converge quickly to highly accurate values for the eigenvalues in (1). Of course, in addition to the eigenvalue , Beyn's method also gives us the eigenvectors , which we can use in postprocessing to do things like visualizing the distribution of currents and fields in the various eigenmodes we find (see examples below).
Beyn's algorithm: Mechanics
Having briefly outlined the idea of Beyn's methodand referring interested readers to Beyn's paper cited above for more detail on the theorylet's now turn to a discussion of the mechanics of running Beynmethod calculations in scuffspectrum.
The following figure illustrates a typical elliptical contour in the complex plane over which we might wish to execute Beyn's algorithm. (Contours in scuffspectrum are always elliptical.)
In this figure, blue crosses indicate eigenvalues of (1); note that there are 3 eigenvalues contained inside the contour labeled as well as several other eigenvalues not enclosed by the contour, which we do not label as they will be ignored by our calculation. The black dots on the contour indicate quadrature points; the values represented by these points are the complex frequencies at which scuffscatter will assemble the matrix and do some linear algebra to compute the integrand of the Beynmethod contour integrals.
Items in red in the figure above indicate userspecified inputs to scuffscatter: these are

the complexvalued frequency at which the contour is centered

the horizontal and vertical radii (halfminor axes) and of the elliptical contour

the number of quadrature points (12 in this case)

an integer which should be greater than or equal to the number of eigenvalues you expect to be found within the contour; if scuffspectrum finds more than this number of eigenvalues, Beyn's method breaks down and must be restarted with a larger value of . (You will get a console warning in this case.)
On the other hand, the algorithm works fine if winds up being larger than the number of eigenvalues found inside the contour, so you should give yourself some leeway by choosing generous values for (typical values might be 10 or 20); higher values of do result in slightly greater computation time, but not much. (More specifically, is the number of linearsystem solves of the form that must be done at each quadrature point ; increasing from 10 to 20 or 30 or so does require more solves, but that cost is generally negligible compared to the cost of assembling and factorizing the matrix , so don't bother with excessive parsimony in choosing .)
Running scuffspectrum with the above inputs as illustrated in the figure above would yield accurate values for the three eigenvalues inside the contourwith the accuracy increasing with the number of quadrature points as well as values for the corresponding eigenvectors.
2. scuffspectrum tutorial: Modes of a spherical dielectric cavity
In this example we'll use the Beyn algorithm as implemented by scuffspectrum to compute the modes of a spherical dielectric cavity, i.e. a simple dielectric sphere.
Exact (sphericalwave) calculation
For this simple geometry, the mode frequencies can be calculated to any desired numerical precision by looking for poles of the Miescattering coefficients for plane waves impinging on a dielectric sphere; the calculation is discussed on page 15 of this memo and implemented by this simple mathematica code. For a sphere of relative permittivity (refractive index ), some typical results are:
 For type (electric multipole) waves with there is a resonance at
where with the sphere radius.
 For type (magnetic multipole) waves with there is a resonance at
Next let's ask how well we can reproduce these results in scuffspectrum.
scuffspectrum calculation
scuffem geometry files for spherical cavity
For scuffspectrum calculations I use a simple .scuffgeo
file describing an isolated sphere with ;
this file is called E4Sphere_501.scuffgeo
:
OBJECT Sphere
MESHFILE Sphere_501.msh
MATERIAL CONST_EPS_4
ENDOBJECT
The file Sphere_501.msh
to which this refers is a
gmsh
mesh file produced from a gmsh
geometry file named Sphere.geo
by running gmsh 2 Sphere.geo o Sphere.msh
; to investigate
the effect of meshing resolution I will also create a series
of more finely discretized spheres by saying e.g.
gmsh 2 Sphere.geo clscale 0.75 o Sphere_Finer.msh
(where the clscale
option, short for "characteristic
length scale", may be used to refine the discretization length
everywhere in the geometry.)
The finestresolution sphere I will use has interior
edges; for comparison , here are images of the coarsest
and finest sphere meshes I will use.
Running scuffspectrum to pinpoint mode frequencies
Based on the discussion above, I expect to find modes at frequencies near and . Thus, I will use scuffspectrum to execute Beyn's method for contours centered near these points and enclosing them.
For example, here's a run in which I specify a contour centered at , with horizontal and vertical radii , using =14 quadrature points, and allowing a budget of up to =5 eigenvalues to be found within the contour:
% scuffspectrum geometry E4Sphere_501.scuffgeo \
omega0 1.10.63i \
Rx 0.2 Ry 0.2 \
N 14 L 5
This produces the file E4Sphere_501.ModeFrequencies
:
# For contour w0=1.1+0.63i, Rx=2.000000e01, Ry=2.000000e01, N=14, L=5# w0=1.1+0.63i, Rx=2.000000e01, Ry=2.000000e01, N=14, L=5:
# re(w) im(w) estimated error in re(w), im(w)
+1.149426e+00 6.378139e01 +3.180420e05 +7.130373e05
As the header says, the 4 numbers reported here are the real and imaginary parts of the mode frequency (which is being estimated at ) followed by estimates of the integration error, i.e. the error incurred by approximating the contour integral via numerical quadrature; in this case the error is evidently tiny, and we can be confident that Beyn's method has converged even with just the 14 quadrature points we used. (For larger contours we would probably need more points.)
Also, in this case Beyn's method only found a single eigenvalue within the contour. If we had considered a larger contour, there would most probably have been several lines of data each with the 4number format shown above, corresponding to multiple eigenvalues identified within the contour.
Mesh convergence
Of course, even if Beyn's method converges to manydigit precision on an eigenmode, we don't necessarily expect the resulting mode frequency to agree with Mietheoretic predictions to high accuracy, because we are here studying a slightly different structurenamely, a sphere approximated by a discretized version in the form of a 501sided polygonbut we may expect the mode frequencies to converge toward exact results for spheres as the discretization is refined.
To test this, I repeated the scuffspectrum calculations above using more finelymeshed spheres. Here are results for mode frequency vs. mesh discretization, characterized by the number N of interior edges in the discretized sphere:
N  Re()  Im() 

501  1.14943  0.63781 
1482  1.14076  0.63318 
2604  1.13877  0.63204 
4107  1.13787  0.63157 
inf  1.13627  0.63071 
exact  1.13622  0.63063 
In this table, the first four data rows show
results of scuffspectrum calculations
for discretized meshes with interior edges.
The column marked inf
corresponds to an extrapolation
of the finite data to the limit
(basically, just fit the data to and
retain only the coefficient.)
The row labeled exact
shows
numerically exact predictions of Mie theory.
Here are the analogous results for the type spherical cavity mode:
N  Re()  Im() 

501  2.09595  0.147858 
1482  2.07996  0.146926 
2604  2.07615  0.146679 
4107  2.07442  0.146564 
inf  2.07151  0.146403 
exact  2.07141  0.146361 
These results demonstrate that the Beyn method implemented by scuffspectrum easily achieves 4digit precision in determining mode frequencies, despite the modest computational burden (equivalent to solving scattering problems at 14 different frequencies). This is much better accuracy than we could hope to get for the same computational cost from a bruteforce frequencyscan method.
Running scuffspectrum to analyze eigenmode field patterns
Moreover, having identified the frequencies (eigenvalues) of our structure with high precision, we can now look at the surface currents (eigenvectors) and the field patterns they produce. scuffspectrum offers several commandline options to facilitate this analysis, many of which are similar to the postprocessing outputs available in scuffscatter and other scuffem codes:

EPFile MyEPFile
requests computation of E and H field components at each evaluation point in a userspecified list of evaluation points. 
FVMesh MyFVMesh.msh
requests generation of gmsh visualization files showing fields and fluxes on the userspecified visualization mesh (screen)MyFVMesh.msh.

CartesianMomentFile MyFile
requests that values of the Cartesian multipole moments corresponding to the eigencurrent distribution be written to fileMyFile.

SphericalMomentFile MyFile
is like the previous option, but for spherical multipole moments. 
PlotSurfaceCurrents
requests visualization files showing the distribution of electric and magnetic surface currents on the geometry.
For example, here are images showing the distribution of surface currents for the two eigenmodes captured by the above tables.
The quadrupole structure of the resonance is clearly distinguishable from the dipole structure of the mode.