1 Introduction
Surface plasmon polaritons (SPPs) are charge density waves that are coupled to electromagnetic (EM) waves at the interface between a metal and a dielectric substrate. Exhibiting strong confinement and relatively low propagation losses, they are thought to be a novel way to confine and control light on the subwavelength scale in the field of nanophotonic technology. Such SPPs can be excited in graphene, a twodimensional carbon allotrope with a single atom layer that is arranged in a honeycomb lattice structure geim04 . It is characterized by strong confinement, low losses, and extreme tunability geim04 ; bludov13 ; In the infrared regime, the electric surface conductivity of such a 2D material is characterized by being complexvalued with a dominant positive imaginary part. This allows for the propagation of SPPs. Plasmons on graphene offer not only a lower ohmic loss than conventional plasmonic materials, but also a strong subwavelength confinement of the EM field liu16 ; oulton09 . The tunability of graphene by electrical gating or chemical doping, makes graphene a promising candidate for future compact plasmon devices vakil11 .
A conventional approach of analyzing a waveguide problem is to first reduce Maxwell’s equations to a Helmholtz eigenvalue problem. For a homogeneously filled waveguide, the EM fields decouple from one another, making the numerical simulation straightforward. However, if spatially dependent material parameters (gradientindex materials) are introduced, the field components are no longer independent from each other, and we are left with a coupled nonlinear eigenvalue problem.
Computational approaches for solving nonlinear eigenvalue problems have been studied in the literature gavin18 ; jay19 ; jay20 ; bai18 . They typically require specialized solvers not readily available in finite element dealII91 and numerical linear algebra tookits petsc04 ; slepc05 . In this paper we pursue a different approach that allows to use existing linear algebra techniques. To that end, we investigate a general class of waveguide configurations that consist of spatially dependent material parameters and contain (arbitrarily shaped) interior conducting 2D material interfaces. In detail, our contributions are as follows:

We derive a variational, nonlinear quartic eigenvalue problem for a waveguiding problem incorporating spatially dependent material parameters and interior conducting interfaces (see Section 2.2). The nonlinear quartic character of the eigenvalue problem stems from the fact that the spatially dependent material parameters cause a strong coupling between the otherwise decoupled transverse magnetic (TM) and transverse electric (TE) modes (as would normally be the case for the Helmholtz equation).

We solve the quartic eigenvalue problem numerically by transforming it to a spectrally equivalent linear system using a quadratification teran14 approach. Additional numerical tools, such as the Möbius transform and a perfectly matched layer (PML) are employed to assist with solving the eigenvalue problem. We verify our numerical method against analytical solutions for prototypical geometries with internal conducting interfaces.

As a practical example, we demonstrate how an improved quality factor (defined by the ratio of the real and the imaginary part of the computed eigenvalues) can be obtained (a) for a family of gradientindex host materials, and (b) by deformation of the geometry (see Section 5.3). Finally, we outline how our framework lays the groundwork for solving related shape optimization problems.
1.1 Related works
Optical properties of cylindrical waveguides with graphene interfaces have been extensively studied in the engineering community liu16 ; xu18 ; gao14a ; gao14b . Recently, focus has also shifted to gradientindex structures that couple with graphene xu16 ; moharrami20 . These structures are based on planar and cylindrical graphenedielectric multilayer metamaterials, and have shown potential applications in terahertz imaging, sensing, detecting, and communication areas xu16 ; moharrami20 . Numerical methods that compute eigenvalues of inhomogeneously loaded domains have been studied before dai13a ; dai13b ; chew99 ; white02
. A finite difference frequencydomain method is used to analyze eigenmodes of inhomogeneously loaded rectangular waveguides in
dai13a . Another study white02presents a method for computing solenoidal eigenmodes and the corresponding eigenvalues of the vector Helmholtz equation. However, neither of these references take into account possible conducting interfaces.
There exist a number of methods for numerically solving nonlinear eigenvalue problems. The FEAST algorithm polizzi09 ; gavin18 uses complex contour integration to compute a cluster of eigenvalues within some userdefined region in the complex plane. Studies have also used the algorithm to simulate propagation of light through optical fibers jay19 ; jay20 . Another technique is within the context of Rayleigh quotient optimization problem bai18 , where a nonlinear eigenvalue problem is solved using a spectral transformation based on nonlinear shifting and a reformulation using secondorder derivatives. In this paper, however, a different method based on quadratification is proposed that is both pragmatic and applicable to a wide range of optical waveguiding problems.
1.2 Paper organization
The remainder of the paper is organized as follows. In Section 2, we derive a variational quartic eigenvalue problem for the waveguiding problem based on timeharmonic Maxwell’s equations. In Section 3, we describe our numerical approach for solving the quartic eigenvalue problem, including a linearization based on quadratification, the use of a Möbius transform to shift the spectrum, and a PML. Section 4 discusses and derives analytical solutions for prototypical configurations, which will be used to verify our numerical method in the subsequent section. Section 5 presents numerical results in domains with and without azimuthal symmetry. We demonstrate how our numerical method can be extended to arbitrary computational domains. Section 6 concludes the paper with a summary of our results and an outlook.
2 Variational formulation
We introduce a variational formulation for a relevant eigenvalue problem prescribed with a gradientindex host material with (arbitrarily shaped) conducting interfaces in the context of a waveguide configuration. A convenient rescaling of the equations to dimensionless forms is used maier17a .
2.1 Preliminaries
The sourcefree timeharmonic Maxwell’s equations are given by
(1) 
where and denote the electric and magnetic field, respectively, and is the temporal frequency. and are complexvalued functions of the transverse coordinates, where denotes the magnetic permeability and denotes the electric permittivity. In order to study guided modes we make the additional ansatz
and decompose the fields and the gradient operator, , into their longitudinal and transverse parts, whence we obtain
where the subscript s denotes the transverse direction and denotes the unit vector in direction. In the strong sense, (1) holds true everywhere except on the points comprising the conducting interface, . The surface conductivity on the conducting interface gives rise to a jump condition on the tangential part of the magnetic field maier17a . In summary, we obtain
(2) 
where is a chosen unit normal vector field on , and denotes the jump over with respect to , viz.,
We also fix the notation
Next we introduce a convenient rescaling of the system to dimensionless form by setting the characteristic wavenumber of the ambient space to 1 maier17a :
To lighten the notation, we omit the breve sign in the remainder of this paper. Applying the rescaling to (2) and rewriting into tangential and normal parts leads to the following interface conditions:
(3) 
where and denote the derivative in the tangential and the normal direction, respectively; is a function in the transverse direction; and where we have used the identities (see Appendix A)
(4) 
The firstorder system (1) can be manipulated in a similar fashion (see Appendix A) to obtain
(5) 
2.2 Variational Statement
Let , where , be a simply connected and bounded domain with Lipschitzcontinous piecewise smooth boundary, . Assume, in addition, that is a Lipschitzcontinuous, piecewise smooth hypersurface. Let and denote the outer normal and the tangential vector on (see Figure 4). Some algebraic manipulation shows that and , which can be used in conjunction with (27) and (5) to obtain:
(6) 
We observe that if the domain is homogeneously filled and isotropic, the curl terms vanish, yielding the familiar decoupled Helmholtz equation for and . Now assume that
(7) 
For the sake of brevity, we summarize the derivation here and refer the reader to Appendix A for details. By multiplying (6) by and testing the first equation by and the second equation by , we obtain
(8) 
This shows the following statement.
3 Numerical approach
In this section, we outline our numerical approach for solving the quartic eigenvalue problem (9). In particular we discuss a quadratification approach transforming the quartic eigenvalue problem into a larger, linear eigenvalue problem with equivalent spectrum. A perfectly matched layer (PML), an artificial sponge layer placed near the boundary such that all outgoing waves decay exponentially, is introduced. The variational formulation (9) is discretized on a nonuniform quadrilateral mesh.
Proposition 2.
Let be a finite element subspace spanned by Lagrange finite elements .
(Q) 
Our goal is to translate (Q) into a finite dimensional linear problem, which then allows the use of a standard linear algebra solver.
3.1 Quadratification of the Quartic Eigenvalue Problem
We build upon the algebraic tool of quadratification introduced in teran14 , which allows us to reduce any even power matrix polynomial eigenvalue problem to a spectrally equivalent linear eigenvalue problem. Prop. 3 summarizes the main result (for a more general discussion of the ideas behind this reduction procedure, we refer the reader to teran14 ).
Proposition 3 (Theorem 5.3 and 5.4 of teran14 ).
Consider a quartic eigenvalue problem, to find , and , s. t.
where are given matrices. Then, the linearization stated below is spectrally equivalent to the original problem (c.f. teran14 Theorem 5.3 and 5.4): Find , and s. t.
(11) 
Here, denotes the identity matrix.
With this result at hand, we rewrite (Q) as a linear eigenvalue problem:
(LQ) 
where
Here, by some abuse of notation denotes the corresponding matrix formed by the bilinear form given in (10) and by fixing a basis of
. A quick computation shows that the eigenvectors of the original problem (
Q) and of the final linearized problem (LQ) are related as follows.Proposition 4.
Let and be an eigenvalue with corresponding eigenvector of (Q). Then, and given by
(12) 
is an eigenvalue with corresponding eigenvector of (LQ). Conversely, if and is an eigenvalue and eigenvector pair of (LQ), then provided that and , the vector is characterized by (12) and and are an eigenvalue and eigenvector pair of (Q).
3.2 Perfectly Matched Layer
A perfectly matched layer (PML) is a truncation procedure motivated from electromagnetic scattering problems in the time domain. The idea is to surround the computational domain with a socalled sponge layer, an artificial boundary wherein all outgoing electromagnetic waves decay exponentially with minimal artificial reflection. As outlined in collino98 ; monk03 ; maier17a , we carry out a change of coordinates from the computational domain with realvalued coordinates to a domain with complexvalued coordinates. We refer the reader to collino98 for details. For a spherical absorption layer, we define the transformation , where
and is an appropriately chosen, nonnegative scaling function. Applying the above transformation, the quartic eigenvalue problem takes the following rescaled form within the PML:
This can be rewritten as
where , and is the rotation matrix that rotates onto . We enforce the condition that the material parameters are constant outside the PML, i.e., and do not undergo a change of coordinate. Additionally, because the eigenmodes of our interest are confined to the conducting interface, which is situated inside the PML, no coordinate change is needed for the jump condition. The modified bilinear forms, , are
(13) 
3.3 Möbius Transform
The Möbius transformation is an efficient tool that transforms the given eigenvalue problem into a related one, for which eigenvalues of interest in the transformed spectrum are close to (or maximal), so that they can be selectively computed with conventional Krylovspace iteration techniques. It is a conformal mapping, defined as follows.
where
. Over arbitrary fields, the Möbius transformation preserves a number of spectral features of matrix polynomials, such as regularity, rank, minimal indicies, the location of zero entries, symmetry, and skewsymmetry
mackey14 . In particular, every Möbius transformation preserves the relation of spectral equivalence mackey14 . The computational eigenvalue problem, after introducing a PML, truncating the domain, and applying finite elements, can be written as(14) 
for a complexvalued vector and appropriate complexvalued matrices and . The implementation of the PML discussed in 3.2 necessitates changes to the definition of , and . The idea is to use the Möbius transform to map points near the origin to target values
(15) 
where are the Möbius transform parameters. The original eigenvalue can be retrieved via the inverse Möbius transform .
4 Validation of weak formulation
Eigenvalues from (L)  Eigenvalues from (9)  Eigenvalues from (19)  

Mode  
1  
2  
3 
In this section, the analytical solution for constant material parameters is derived and discussed. We use the analytic result to validate our numerical approach.
By assuming constant material parameters, the quartic eigenvalue problem (9) does not exhibit any hybridization and reduces to a linear eigenvalue problem: Find s. t.
(L) 
for and where we have introduced the bilinear forms
(16) 
For a simple spherical geometry that is rotationally invariant, the field solution can be expressed as a superposition of the modified Bessel functions of the first and second kind. In the case of a waveguide with a single circular interface , i.e., where is described by a circle with origin and radius , the analytical solution takes the following form.
(17)  
(18) 
where and denote the modified Bessel functions of the first and second kind, respectively, and , and are constants that are determined by the boundary conditions gao14a ; liu16 . Assuming that the conducting film is located on the boundary of the interior circle with radius , we equate the jump conditions (3) of each field component. Then (3) reduces to the following algebraic condition, from which we can retrieve the propagation constant, :
(19) 
where
and
We numerically solve the zeroes of the determinant of the difference of the above matrices via a root finding algorithm for modal orders and . The computed values are then compared against those of the linear problem (L) and of the quartic problem (9) (see Table 1 and Figure 11).
We now validate our numerics with the analytical results. Three validations are carried out: analytical, numerical with linear eigenvalue, and lastly, the quartic eigenvalue problem. For simplicity, we assume the computational domain is isotropic, with material parameters . A conducting interface is coated on the boundary of the interior circle with radius , and choose . The surface conductivity is set to .
The analytic computation of the propagation constant requires finding the complex roots of the determinant of the matrix. We will consider modal orders of and to ease the computation. The computational results displayed in Table 1 deviate by less than 1% from each other. We can thus expect a confidence level of 1% or better in our numerical computations. Figure 11 shows the intensity of the numerically computed electric field,, and a comparison of, both, the analytic and numerical solutions. We thus conclude that our numerical framework is a reliable model that can effectively simulate hybrid plasmonic modes.
5 Numerical results
These eigenvalues are numerically computed using SLEPc slepc05 , a general purpose eigensolver built on top of PETSc petsc04 . The eigensolver provides a number of Krylovspace methods, such as the Arnoldi, Lanczos, KrylovSchur, and conjugategradient methods. For our purposes, we will make use of the KrylovSchur method for its faster convergence.
In this section, we present a number of computational results obtained from solving the quartic eigenvalue problem (9). We examine numerically how the spectrum of a hybridized configuration behaves under modification of spatially dependent material parameters, . We further investigate the relationship between mesh deformation and a quality factor, and study the degree by which the spectrum changes. All numerical computations are carried out with the finite element library deal.II dealII91 . We use a KrylovSchur method to compute solutions of the linearized eigenproblem (LQ) (slepc05, ).
We demonstrate numerically how it is possible to attain an improved quality factor by prescribing the host material with a radiallyvarying refractive index profile. This is methodically carried out in the subsequent subsections. First, a parameter study is conducted to validate our choice of discretization parameters. Second, we solve the quartic eigenvalue problem (LQ) using a number of permittivity functions, and observe how the spectrums differ from those obtained in an isotropic medium. Lastly, we deform our computational domain to demonstrate that our numerical framework is equipped to handle even the most general configuration. The key idea behind this generalization is to show that we can manipulate the spectrum by manipulating the shape. We make note of the evolution of eigenmodes, and how our framework can be used as a basis for shape optimization of gradientindex waveguides.
5.1 Validation of discretization parameters



The computational domain, , is chosen to be the circle with radius 1. A spherical PML is enforced for . The surface conductivity is chosen that is within a realistic parameter range maier17a . Following monk03 , the nonnegative scaling function is chosen to be
(20) 
where we set the free parameter to be in our computations. We carry out a parameter study to test the validity and the sensitivity of discretization parameters. Table 5 summarizes the parameter study quantitatively. As can be seen, the eigenmodes comptued are stable with respect to variations of PML strength , the number of initial refinements , and domain sizes . A spectral transformation is carried out in the form of the Möbius transformation, with the parameters chosen to be
Comments
There are no comments yet.