ADVANCED ELECTROMAGNETICS, V�����, N����,�������� � ���� An Electric Field Volume Integral Equation Approach

In this paper we present an electric field volume integral equation approach to simulate surface plasmon propagation along metal/dielectric interfaces. Metallic objects embedded in homogeneous dielectric media are considered. Starting point is a so-called weak-form of the electric field integral equation. This form is discretized on a uniform tensor-product grid resulting in a system matrix whose action on a vector can be computed via the fast Fourier transform. The GMRES iterative solver is used to solve the discretized set of equations and numerical examples, illustrating surface plasmon propagation, are presented. The convergence rate of GMRES is discussed in terms of the spectrum of the system matrix and through numerical experiments we show how the eigenvalues of the discretized volume scattering operator are related to plasmon propagation and the medium parameters of a metallic object.


Introduction
While photons and phonons are quanta of energy for light and mechanical vibrations, plasmons result from the quantization of plasma oscillations in a conductor.Plasmons are not elementary particles, like photons, but quasiparticles, like phonons; nonetheless, they can be used to transfer energy and information, when coupled to a photon to create a polariton [1,2].Surface plasmon polaritons (SPPs) occur usually at the interface between a dielectric and a conductor [3,4].
In this paper SPPs are studied, in particular how they form, propagate, and dissipate in a variety of solid-state configurations mimicking the surface of a conventional CMOS structure.The paper outlines how the equations governing SPPs can be derived from Maxwell's equations and it describes the engine simulating SPP transients in details.The goal is to achieve the equivalent of a device simulator for SPPs for guiding researchers in the design process of data processing systems based on SPPs.
It is well known that SPPs can be excited by electromagnetic fields in three dimensions and by H-polarized fields in two dimensions (magnetic field strength parallel to the invariance direction) [4,5].Local and global solution methods are available to simulate these SPPs along a variety of metallic objects.In a local method, Maxwell's equations are discretized directly (as in the well known Finite-Difference Time-Domain method), whereas in a global method an integral form of these equations is solved.In this paper we use a global volume integral approach, since we are interested in electromagnetic fields operating in steadystate and electromagnetic field strength unknowns are defined on the scattering domain only.In addition, outward radiation is automatically taken into account via the Green's tensor of the background medium and there is no need to implement any absorbing boundary conditions as is the case in a local method.A global integral equation approach also allows us to analyze plasmon propagation in terms of the spectrum of the volume integral operator.
For three-dimensional electromagnetic fields and for H-polarized fields in two dimensions, the electric field strength satisfies a vectorial volume integral equation containing a gradient-divergence term that takes the effects of induced charges into account.It is well known that this term must be handled with care when discretizing the integral equation [4,6].In this paper, we follow the approach proposed, for example, in [6] and approximate the gradient-divergence term by second-order centered finitedifferences.Furthermore, the singular Green's function is weakened in a manner that is consistent with the finitedifference approximation error and we make use of the Kronecker product to show that the discretized set of equations has a similar structure as the continuous volume integral operator.Details of the discretization procedure are provided, but to keep the bookkeeping to a minimum, we restrict ourselves to two-dimensional configurations only.Threedimensional problems can be handled in a similar manner.
This paper is organized as follows.In Section 2, we briefly review the basic integral representation for the scattered electric field strength and formulate an integral equation for the total electric field strength inside the object of interest.This integral equation is discretized in Section 3 and a detailed description of the discretization procedure and the structure of the discretized set of equations is provided.In Section 4 we illustrate the performance of our integral equation solver through a number of numerical ex-periments in which we simulate surface plasmon propagation along a golden strip and a silver plasmonic waveguide.The discretized set of equations is solved using the Generalized Minimum Residuals (GMRES) iterative solver [7] and a numerical analysis of its convergence rate for plasmonic configurations is presented as well.

Basic equations
We consider steady-state H-polarized fields in a twodimensional configuration that is invariant in the zdirection.A (metallic) object occupies a bounded domain D obj in the transverse xy-plane (see Figure 1) and is characterized by a position dependent admittance per unit length η s (x) = σ s (x) + jωε s (x) and constant impedance per unit length ζ = jωµ 0 .In these expressions, σ s (x) and ε s (x) are the conductivity and the permittivity of the object, respectively, and µ 0 is the permeability of vacuum.
The object is embedded in a homogeneous background medium that is characterized by a constant conductivity σ, permittivity ε, and permeability µ 0 .We write the corresponding admittance and impedance per unit length as η = σ + jωε and ζ = jωµ 0 , respectively.Furthermore, the wave number of the background medium is introduced as with Im(k b ) ≤ 0. Now let the external sources occupy a bounded domain D src with D src ∩ D obj = ∅ (see Figure 1).To find the electric field strength in and around the object, we make use of the linearity of Maxwell's equations and set up a scattering formalism.Specifically, we first introduce the incident electric field as the field that would be present in our domain of interest if the object was absent.We denote this field by E inc and assume that it is known.Second, we take the presence of the object into account by introducing the scattered electric field as E sc = E − E inc , where E is the total electric field strength.The scatterer acts as a source for the scattered electric field and it is well known that this field can be found at any point in the transverse plane from the integral representation [8] (2) In the above equation, G is the Green's function of the homogeneous background medium and is given by where H (2) 0 is the Hankel function of the second kind and order zero.This function has a logarithmic singularity at the origin [9].Furthermore, is the contrast function of the object and it obviously vanishes if η s = η, that is, if there is no object present.From the integral representation (2) it is clear that the scattered electric field strength is known as soon as the total electric field strength inside the object has been found.To find this total field, we restrict the observation vector x to the object domain D obj and use the definition of the scattered electric field strength to obtain where we have introduced the vector potential A as Equation ( 5) is an integral equation for the total electric field strength inside the scattering object.Written out in components, this equation becomes and with x ∈ D obj .As mentioned above, by solving the above integral equation for the electric field strength, we have essentially solved the complete scattering problem, since the scattered field at any point in the transverse plane can then be found using the integral representation of Eq. (2).

The contrast function for metallic objects
For metallic objects, it is customary to use the relative complex permittivity ε r as a constitutive parameter.Writing this parameter in terms of its real and imaginary part as where ε 1 and ε 2 are real-valued, we have Now if the background medium is vacuum, we have η = jωε 0 and the contrast function becomes showing that the contrast is simply a shifted version of the relative permittivity of the metal.
As a first step in the discretization procedure, we require that Eqs. ( 7) and ( 8) hold for all grid points within the scattering domain with position vectors x m,n , where m = 1, 2, ..., M and n = 1, 2, ..., N .Consequently, we have and for m = 1, 2, ..., M and n = 1, 2, ..., N .Notice that the discretization cells located at the outer boundary are not included.These cells are used as dummy cells to properly handle the discretization of the gradient-divergence term (see Subsection 3.1).
As a second step, we approximate the spatial derivatives in Eqs. ( 15) and ( 16) by second-order centered finite-differences.For example, the mixed derivative term and for the double derivative term (one-dimensional Laplacian) in the x-direction, we have Similar finite-difference expressions are used for ∂ y ∂ x A x and ∂ 2 y A y at x = x m,n .All these local finite-difference approximations can be written in a more compact global form by introducing a number of field and differentiation matrices.Specifically, let us introduce the (M + 2)-by-(N + 2) electric field matrices E x and E y with elements and for m = 1, 2, ..., M + 2 and n = 1, 2, ..., N + 2. Incident electric field matrices E inc x and E inc y are defined similarly.In addition, we introduce the vector potential matrices A x and A y with elements and for m = 1, 2, ..., M + 2 and n = 1, 2, ..., N + 2. Differentiation in the x-direction is carried out by the M -by-(M + 2) differentiation matrices and while the N -by-(N + 2) matrices and take care of differentiation in the y-direction.The first entry between the brackets in the above equations corresponds to a diagonal element.Finally, we introduce the N -by-(N + 2) restriction matrix where I N is the identity matrix of order N .
With the introduction of all these matrices, we can write the finite-difference approximation of Eqs. ( 15) and ( 16) in global form as and respectively.
As a final step, we turn these equations into matrixvector form by applying the vec-operation to both equations.Recall that for a given matrix A, vec(A) stacks the columns of matrix A from left to right into a single column vector [10].Using this operation, we first introduce the vectors and a x = vec(A x ) and a y = vec(A y ).
Applying now the vec-operation to Eqs. ( 23) and (24) separately, using the linearity of this operator, and the property [10] vec where ⊗ denotes the Kronecker (tensor) product, we arrive at and where R = R N ⊗ R M .The only thing left to do now, is to relate the vector potentials a x and a y to the electric field strengths e x and e y using the definition of the vector potential.

Weak form of the vector potential and its discretization
With x ∈ D obj , the expression for the vector potential cannot be discretized using standard Newton-Cotes quadrature formulas, because of the logarithmic singularity of the Green's function.Our approach is, therefore, to first weaken the Green's function such that this weakening is consistent with the second-order finite-difference discretization procedure discussed above [6].In particular, we introduce the weakened Green's function G w that satisfies and the radiation condition at infinity.The right-hand side f (x) is given by where D circ is a circular disk with radius a = 1 2 min{δx, δy} and centered at the origin.Clearly, f satisfies and G w approaches the original Green's function G as a ↓ 0.
Replacing now the original Green's function by its weakened counterpart, we have and the expression on the right-hand side of the above equation is used to approximate the vector potential at the grid nodes with position vectors x = x m,n .
To this end, we consider for m = 0, 1, ..., M + 1 and n = 0, 1, ..., N + 1. Subsequently, we restrict ourselves to piecewise constant contrast functions for which χ is constant within each discretization cell.Writing where χ i,j is position independent, we obtain for m = 0, 1, ..., M +1 and n = 0, 1, ..., N +1.Finally, we approximate the remaining integral over the discretization cell using the midpoint rule and arrive at for m = 0, 1, ..., M + 1 and n = 0, 1, ..., N + 1.Notice that the vector potential approximations are required at grid nodes with indices running from m = 0 to m = M + 1 and n = 0 to n = N + 1.The summations in Eq. ( 35) run from i = 0 to M + 1 and from j = 0 to j = N + 1 as well.The electric field strength unknowns E(x m,n ), however, are defined at grid node coordinates x = x m,n with m = 1, 2, ..., M and n = 1, 2, ..., N .To relate the vector potential to these electric field strength unknowns only, we always set the contrast to zero at the outer cells of the computational domain, that is, we set χ 0,n = 0 and χ M +1,n = 0 (36) for n = 0, 1, ..., N + 1, and χ m,0 = 0 and χ m,N +1 = 0 (37) for m = 1, ..., M .Equation ( 35) can then equivalently be written as for m = 0, 1, ..., M + 1 and n = 0, 1, ..., N + 1.However, to evaluate the vector potential values at the grid nodes we do not use Eq. ( 38), but follow an embedding procedure and implement Eq. ( 35) with the contrast at the boundaries set to zero.The electric field strength quantities at the outer cells then become redundant dummy variables, but the advantage of keeping these extra variables is that all approximate vector potential values can be computed very efficiently by evaluating Eq. ( 35) via two-dimensional FFTs.This is not possible if we use Eq. ( 38) to evaluate the vector potential at the required grid nodes.
From Eq. (35), we observe that the weakened Green's function is required at grid nodes outside the circular disk D circ and at the origin.To find these function values, we return to Eq. (30).The solution of this equation is given by and evaluating the integral for x / ∈ D circ , we find (see Appendix A) where J 1 is the first-order Bessel function of the first kind.Notice that the above expression for the weakened Green's function is just a scaled version of the original Green's function.To show that this scaling is consistent with the second-order finite-difference approximations of Eqs. ( 17) and (18), we use the series expansion of J 1 (k b a) around zero [9] and obtain showing that the weakening procedure is indeed consistent with the second-order finite-difference discretization procedure.Finally, at the origin the integral evaluates to (see Appendix A) and we call this weakened Green's function value the self patch element.With this result, we have completely determined the discretized vector potential operator of Eq. (35).To write this operator in matrix-vector notation, we introduce the contrast vector c and contrast matrix C as c = vec(χ i,j ) and C = diag(c). (42) In addition, we introduce the matrices Gj,n , j, n = 0, 1, ..., N + 1, of order M + 2 with elements for i, m = 0, 1, ..., M + 1.With the introduction of these matrices, we can write where matrix G is given by Using Eq. ( 43), it is easily verified that matrix G is block Toeplitz and each block is Toeplitz as well (matrix G is a so-called Block-Toeplitz Toeplitz-Block or BTTB matrix).This implies that its action on a vector can be computed via two-dimensional FFTs.Furthermore, matrix G is symmetric (but not Hermitian) and all self patch elements are located on the diagonal.Substituting Eq. ( 44) in Eqs. ( 28) and (29), we arrive at the system of equations where the system matrix is given by Since the action of matrix K on a vector can be computed efficiently via FFTs, we solve Eq. ( 46) using an iterative method such as GMRES or BiCGStab [7].

Simulations
To illustrate the performance of our integral equation approach, we present two numerical experiments in which we simulate surface plasmon propagation.In both examples, the electromagnetic field is computed by iteratively solving Eq. ( 46) using the GMRES algorithm [7].This algorithm takes an initial guess u 0 as input and generates field approximations u 1 , u 2 , ..., such that the Euclidean norm of the corresponding residual is minimized at every iteration.In our numerical experiments, we take u 0 = 0 as an initial guess leading to an initial residual r 0 = u inc .Iterations are terminated as soon as the normalized Euclidean norm of the residual r n / r 0 falls below a user specified tolerance.As a first example, we compute the electromagnetic field in and around a golden strip (see Figure 2).The strip is illuminated from the left by electromagnetic waves that are generated by a line source operating in steady-state at a frequency of 4.73 • 10 14 Hz (λ = 633 nm).The complex relative permittivity of gold is ε r = −11.6 − 1.2j at this frequency and the iteration process is terminated as soon as the normalized residual falls below 10 −6 .A snapshot of the instantaneous magnitude of the electric field strength is shown in Figure 3 and the convergence history of GMRES is shown in Figure 4. We observe that initially there is a fairly sharp drop in the norm of the residual, but convergence slows down to a steady rate of decrease as GMRES    proceeds (after about 50 iterations in the above example).
To study this effect, we first recall that GMRES constructs field approximations u k such that [7] where P k is the set of polynomials of degree ≤ k.In other words, of all polynomials p k belonging to P k and normalized such that p k (0) = 1, GMRES constructs the polynomial for which p k (K)r 0 is minimum.The roots of this optimal polynomial are called harmonic Ritz values [11] and approximate the eigenvalues of the system matrix K.
We have computed these eigenvalue estimates after 50 and 100 iterations of the GMRES algorithm.We also computed the first 30 eigenvalues of the system matrix with the smallest imaginary part.These eigenvalues are shown as circles in Figures 5 and 6.In Figure 5, the crosses show the first 50 harmonic Ritz values obtained after 50 GMRES iterations, while the first 100 harmonic Ritz values are shown as crosses in Figure 6.From these figures, we observe that the harmonic Ritz values cluster along curves in the complex plane.Eigenvalues of banded Toeplitz matrices are known to have this property [11], but here we have a full BTTB-type matrix.Also note that the leftmost point of the upperleft spectral line corresponds to the complex relative permittivity ε r = −11.6 − 1.2j of the golden strip (see Figures 5 and 6).Finally, we observe that GMRES quickly approximates the eigenvalues located at the outer boundary of the spectrum as is typical for Krylov subspace methods [11].The inner eigenvalues are approximated as GM-RES proceeds and the harmonic Ritz values start to cluster along curves in the complex plane.Convergence slows down to a steady rate of decrease as soon as GMRES starts "filling" these curves.
In our second set of experiments, we consider a plasmonic waveguide consisting of two silver strips and an air gap (see Figure 7).The same source as in the previous example is placed to the left of the waveguide and it again operates at a wavelength of 633 nm.The complex relative permittivity of silver at this frequency is ε r = −18.2− 0.5j [4].Since the total waveguiding structure consists of air and silver and since the contrast function of silver has a larger absolute real part and a smaller absolute imaginary part (smaller losses) compared with gold, we expect that our solver converges at a slower rate than in the previous example.Figure 8 shows that this is indeed the case.After 750    GMRES iterations, the normalized residual has dropped to 1.73 • 10 −3 , while for the golden strip problem 500 iterations are sufficient to obtain a normalized residual that falls below 10 −6 .
To illustrate how the field approximations improve as iteration proceeds, we show in Figures 9 and 10 snapshots of the magnitude of the electric field strength in and around the waveguide.Figure 9 shows the approximate electric field strength obtained after 100 iterations, while Figure 10 shows the approximate field obtained after 500 iterations.Comparing both figures, we observe that as the number of iterations increases, improvements go from left to right, that is, from the source (left-hand side of Figures 9 and 10) to the far end of the waveguide (right-hand side of Figures 9  and 10).The plasmonic wave is built up along the waveguide as iteration proceeds.
The first 30 eigenvalues of the system matrix with the smallest imaginary part have also been computed for this waveguide problem (circles in Figures 11 and 12) along with the harmonic Ritz values obtained after 50 and 100 iterations of GMRES (crosses in Figures 11 and 12, respectively).Again, we observe that eigenvalues at the outer boundary of the spectrum are approximated first, although after 50 iterations GMRES misses an extremal eigenvalue (see Figure 11).Furthermore, the harmonic Ritz values again cluster along curves in the complex plane.As opposed to the previous example, it is now difficult to recognize any material parameter values from the spectral curves in the complex plane, since the total scattering object is now inhomogeneous.The spectral analyses of this section do seem to indicate, however, that the clustering of harmonic Ritz values along curves may provide a spectral explanation of how a plasmonic wave along a metal/dielectric interface is approximated by an iterative solver such as GMRES.
Figure 12: The first 30 eigenvalues of the system matrix with the smallest imaginary part (circles) and the first 100 harmonic Ritz values (crosses) for the plasmonic waveguide problem.

Conclusions
In this paper, we have discussed a volume integral equation method to simulate propagation of surface plasmon polaritons in media that mimic the silicon-to-silicon dioxide interface on standard CMOS integrated circuits.We started from Maxwell's equations and formulated an integral equation for the electric field strength which contains a gradient-divergence term that takes the induces charges into account.The discretization of the integral equation was discussed in detail and by formulating the discretized set of equations in terms of Kronecker products, it was shown that the discretized system matrix has a similar structure as the continuous volume scattering operator.Moreover, the action of the system matrix on a vector can be computed via FFTs and we therefore used GMRES to iteratively solve the discretized integral equation.The convergence rate of the algorithm was discussed in detail, while numerical examples illustrated surface plasmon generation, propagation, and decay in different media.
Even though the simulator efficiently solves practical surface plasmon problems, we can still further improve its performance by preconditioning the discretized integral equation.Future work will focus on the development of an effective preconditioner.This may be beneficial in particular for surface plasmon simulations in electrically large geometries (geometries with a spatial extent in the order of, say, a hundred wavelengths).The eigenvalue analysis presented in this paper may then be helpful, since knowledge about the spectrum or pseudospectrum of the system matrix may provide us a means to construct an effective preconditioner [11].
We believe that the simulator presented in this paper is a practical and efficient tool to simulate surface plasmons and is especially useful to guide optimization of CMOS processes and geometries for plasmonics based data and energy transfer and miniaturization.

A. Explicit expression for the weak Green's function
To evaluate the integral in Eq. (39), we make use of the addition theorem for the Hankel function H (2) 0 (see, for example [9]): for x / ∈ D circ .To evaluate G w (0), we take x ∈ D circ with |x| = < a and subsequently take the limit ↓ 0. Integration over D circ is split in two parts: in one part we integrate from |x | = to |x | = a (domain D 1; ), while in the second part we integrate over a disk centered at the origin and having a radius (domain D 2; ).Notice that when integrating over Using now the power series expansions of J 0 (x) and H (2) 1 (x) around x = 0, namely, and

Figure 1 :
Figure 1: A metallic object illuminated by electromagnetic waves.The object occupies the object domain D obj .The sources are located in the source domain D src .
47) and the field vectors are given by u = e x e y and u inc = e inc x e inc y .

Figure 2 :
Figure 2: Line source located next to a golden strip.

Figure 3 :
Figure 3: Snapshot of the instantaneous magnitude of the electric field strength in and around the golden strip.

Figure 4 :
Figure 4: Convergence history of GMRES for the golden strip problem.

Figure 5 :
Figure 5: The first 30 eigenvalues of the system matrix with the smallest imaginary part (circles) and the first 50 harmonic Ritz values (crosses) for the golden strip problem.

Figure 6 :
Figure 6: The first 30 eigenvalues of the system matrix with the smallest imaginary part (circles) and the first 100 harmonic Ritz values (crosses) for the golden strip problem.

Figure 7 :
Figure 7: Line source located next to a plasmonic waveguide consisting of two silver strips.

Figure 8 :
Figure 8: Convergence history of GMRES for the plasmonic waveguide problem.

Figure 9 :
Figure 9: Snapshot of the instantaneous magnitude of the electric field strength in and around the plasmonic waveguide after 100 GMRES iterations.

Figure 10 :
Figure 10: Snapshot of the instantaneous magnitude of the electric field strength in and around the plasmonic waveguide after 500 GMRES iterations.

Figure 11 :
Figure 11: The first 30 eigenvalues of the system matrix with the smallest imaginary part (circles) and the first 50 harmonic Ritz values (crosses) for the plasmonic waveguide problem.

( 2 )
k (k b |x|) exp(jkϕ), |x | ≤ |x|, (51)where J k is the kth order Bessel function of the first kind, H(2) kis the kth order Hankel function of the second kind, and ϕ is the angle between x and −x .Let us first evaluate the integral for observation points outside the circular disk D circ .In this case, we have |x| > |x | and substituting the second line of Eq. (51) in Eq. (39), we obtainG w (x) = − j G w (x) = − j 2k b a J 1 (k b a)H