# Nonlinear eigenvalue problems for coupled Helmholtz equations modeling gradient-index graphene waveguides

We discuss a quartic eigenvalue problem arising in the context of an optical waveguiding problem involving atomically thick 2D materials. The waveguide configuration we consider consists of a gradient-index (spatially dependent) dielectric equipped with conducting interior interfaces. This leads to a quartic eigenvalue problem with mixed transverse electric and transverse magnetic modes, and strongly coupled electric and magnetic fields. We derive a weak formulation of the quartic eigenvalue problem and introduce a numerical solver based on a quadratification approach in which the quartic eigenvalue problem is transformed to a spectrally equivalent linear eigenvalue problem. We verify our numerical framework against analytical solutions for prototypical geometries. 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 for a family of gradient-index host materials with internal conducting interfaces. We outline how this result lays the groundwork for solving related shape optimization problems.

There are no comments yet.

## Authors

• 1 publication
• 6 publications
• 4 publications
08/30/2021

### Numerical Eigensolver for Solving Eigenmodes of Cavity Resonators Filled With both Electric and Magnetic Lossy, Anisotropic Media

This article presents the numerical eigensolver to find the resonant fre...
05/06/2021

### Vibration Analysis of Piezoelectric Kirchhoff-Love Shells based on Catmull-Clark Subdivision Surfaces

An isogeometric Galerkin approach for analysing the free vibrations of p...
08/19/2019

### Unfitted Nitsche's method for computing wave modes in topological materials

In this paper, we propose an unfitted Nitsche's method for computing wav...
11/05/2021

### Numerical Approximation of Optimal Convex and Rotationally Symmetric Shapes for an Eigenvalue Problem arising in Optimal Insulation

We are interested in the optimization of convex domains under a PDE cons...
03/18/2020

### Computation of Tight Enclosures for Laplacian Eigenvalues

Recently, there has been interest in high-precision approximations of th...
01/05/2020

### Approximation of PDE eigenvalue problems involving parameter dependent matrices

We discuss the solution of eigenvalue problems associated with partial d...
12/01/2020

### Bifurcation Analysis of the Eigenstructure of the Discrete Single-curl Operator in Three-dimensional Maxwell's Equations with Pasteur Media

This paper focuses on studying the bifurcation analysis of the eigenstru...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 two-dimensional 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 complex-valued 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 (gradient-index 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 gradient-index 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 gradient-index structures that couple with graphene xu16 ; moharrami20 . These structures are based on planar and cylindrical graphene-dielectric 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 frequency-domain method is used to analyze eigenmodes of inhomogeneously loaded rectangular waveguides in

dai13a . Another study white02

presents 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 user-defined 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 second-order 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 time-harmonic 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 gradient-index 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 source-free time-harmonic Maxwell’s equations are given by

 (1) {∇×E=iωμH,∇×H=−iωεE,

where and denote the electric and magnetic field, respectively, and is the temporal frequency. and are complex-valued 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

 F∼eikzz,

and decompose the fields and the gradient operator, , into their longitudinal and transverse parts, whence we obtain

 F=Fs+^zFz,∇=∇s+ikz^z,

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) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩[ν×(Hs+^zHz)]Σ=σΣ(ν×(Es+^zEz))×ν,[ν⋅μ(Hs+^zHz)]Σ=0,[ν×(Es+^zEz)]Σ=0,[ν⋅ε(Es+^zEz)]Σ=1iω∇⋅(σΣE),

where is a chosen unit normal vector field on , and denotes the jump over with respect to , viz.,

 [F]Σ(x):=limα↘0(F(x+αν)−F(x−αν))x∈Σ.

We also fix the notation

 F+:=limα↘0F(x+αν),F−:=limα↘0F(x−αν),forx∈Σ.

Next we introduce a convenient rescaling of the system to dimensionless form by setting the characteristic wavenumber of the ambient space to 1 maier17a :

 x→˘x=k0x,∇s→˘∇s=1k0∇s, μ→μr=1μ0μ,\e→εr=1\e0ε,σΣ→σΣr=√μ0\e0σΣ E→˘E=k20ωμ0E,H→˘H=k0H,kz→˘kz=kzk0.

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) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩Σ⋅τ=[ik2s(kz∂τHz+εr∂νEz)]Σ=σΣrEz,[Hz]Σ=−σΣrEs⋅τ,[μrHs]Σ⋅ν=[iμrk2s(kz∂νHz−εr∂τEz)]Σ=0,[Es]Σ⋅τ=[ik2s(kz∂τEz−μr∂νHz)]Σ=0,[Ez]Σ=0,[εrEs]Σ⋅ν=[iεrk2s(kz∂νEz+μr∂τHz)]Σ=1i∇s⋅(σΣrEs),

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) {k2sEs=i(kz∇sEz+μr∇s×Hz),k2sHs=i(kz∇sHz−εr∇s×Ez).

The first-order system (1) can be manipulated in a similar fashion (see Appendix A) to obtain

 (5) {∇s×(μ−1r∇s×^zEz)+ikz∇s×(μ−1r^z×Es)−εrEz=0,∇s×(ε−1r∇s×^zHz)+ikz∇s×(ε−1r^z×Hs)+μrHz=0.

### 2.2 Variational Statement

Let , where , be a simply connected and bounded domain with Lipschitz-continous piecewise smooth boundary, . Assume, in addition, that is a Lipschitz-continuous, 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) ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩−∇s⋅(εrk2s∇sEz)−kz∇s⋅(1k2s∇s×^zHz)−εrEz=0,−∇s⋅(μrk2s∇sHz)+kz∇s⋅(1k2s∇s×^zEz)−μrHz=0.

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) [εr]Σ=0,[μr]Σ=0.

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) (μrε2r∇sEz,∇sφ)+2(εr∇sEz,∇s(¯¯¯εr¯¯¯μr)φ)+kz(μrεr∇s×^zHz,∇sφ)−k2z(εr∇sEz,∇sφ)+2kz(∇s×^zHz,∇s(¯¯¯εr¯¯¯μr)φ)−k3z(∇s×^zHz,∇sφ)−(εrk4sEz,φ)+(εrμ2r∇sHz,∇sψ)+2(μr∇sHz,∇s(¯¯¯εr¯¯¯μr)ψ)−kz(μrεr∇s×^zEz,∇sψ)−k2z(μr∇sHz,∇sψ)−2kz(∇s×^zEz,∇s(¯¯¯εr¯¯¯μr)ψ)+k3z(∇s×^zEz,∇sψ)−(μrk4sHz,ψ)−i⟨σΣrk4sEz,φ⟩Σ=0.

This shows the following statement.

###### Proposition 1.

Provided that and , the nonlinear eigenvalue problem

 (N) Find kz∈C∖{0} and Ez, Hz s. t. (6) and (3) are satisfied

can be restated as a quartic eigenvalue problem

 (9) Find kz∈C∖{0} and Ez, Hz∈X(Ω)={(u,v):u,v∈H1(Ω,C)} s.\,t.Q(kz,(Ez,Hz))[(φ,ψ)]=0for all φ, ψ∈H1(Ω,C),

where

 Q(kz,(Ez,Hz))[(φ,ψ)]=4∑l=0(kz)lal((Ez,Hz))[(φ,ψ)],

and

 (10)

## 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 non-uniform quadrilateral mesh.

###### Proposition 2.

Let be a finite element subspace spanned by Lagrange finite elements .

 (Qh) 4∑l=0klzal((Ez,Hz),(φ,ψ))=0,

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.

 4∑k=0λkAkx=0,

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) ⎛⎜ ⎜ ⎜⎝A30−In0A100−In0−In00A0000⎞⎟ ⎟ ⎟⎠z+λ⎛⎜ ⎜ ⎜⎝A4000A2In0000In0000In⎞⎟ ⎟ ⎟⎠z=0.

Here, denotes the identity matrix.

With this result at hand, we rewrite (Q) as a linear eigenvalue problem:

 (LQh) Sz+λMz=0,

where

 S=⎛⎜ ⎜ ⎜⎝a30−In0a100−In0−In00a0000⎞⎟ ⎟ ⎟⎠,M=⎛⎜ ⎜ ⎜⎝a4000a2I0000I0000I⎞⎟ ⎟ ⎟⎠,z=⎛⎜ ⎜ ⎜⎝z1z2z3z4⎞⎟ ⎟ ⎟⎠.

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) ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩z1=λx,z2=λ2(a3+λa4)x,z3=λ(a3+λa4)x,z4=−a0x,

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 so-called 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 real-valued coordinates to a domain with complex-valued coordinates. We refer the reader to collino98 for details. For a spherical absorption layer, we define the transformation , where

 d=1+is(r),¯¯¯d=1+i/r∫rρs(τ)dτ,

and is an appropriately chosen, nonnegative scaling function. Applying the above transformation, the quartic eigenvalue problem takes the following rescaled form within the PML:

 ^∇s⋅(k2s(εr^∇s^Ez))−2k2s^∇s⋅(εr^∇s^Ez)+kz^∇s⋅(k2s^∇s×^Hz)−2kzk2s^∇s⋅(^∇s×^Hz)+^∇s⋅(k2s(μr^∇s^Hz))−2k2s^∇s⋅(μr^∇s^Hz)+kz^∇s⋅(k2s^∇s×^Ez)+2kzk2s^∇s⋅(^∇s×^Ez)−εrk4s^Ez−μrk4s^Hz=0.

This can be rewritten as

 1d¯¯¯d∇s⋅(k2sεrA∇sEz)−2k2sd¯¯¯d∇s⋅(εrA∇sEz)+kzd¯¯¯d∇s⋅(k2s∇s×Hz)−2kzk2sd¯¯¯d∇s⋅(∇s×Hz)+1d¯¯¯d∇s⋅(k2sεrA∇sEz)−2k2sd¯¯¯d∇s⋅(εrA∇sEz)+kzd¯¯¯d∇s⋅(k2s∇s×Hz)+2kzk2sd¯¯¯d∇s⋅(∇s×Ez)+εrk4sEz+μrk4sHz=0,

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) ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩~a0((Ez,Hz),(φ,ψ))=(ε2rA∇sEz,∇sφ)+2(εrA∇sEz,(∇s¯¯¯εr)φ)−d¯¯¯d(ε3rEz,φ)+(εrA∇sHz,∇sψ)+2(A∇sHz,(∇s¯¯¯εr)ψ)−d¯¯¯d(ε2rHz,ψ)−i⟨σΣrε2rEz,φ⟩Σ,~a1((Ez,Hz),(φ,ψ))=(εr∇s×^zHz,∇sφ)+2(∇s×^zHz,(∇s¯¯¯εr)φ)−(εr∇s×^zEz,∇sψ)−2(∇s×^zEz,(∇s¯¯¯εr)ψ)~a2((Ez,Hz),(φ,ψ))=−(εrA∇sEz,∇sφ)+2d¯¯¯d(ε2rEz,φ)−(A∇sHz,∇sψ)+2d¯¯¯d(εrHz,ψ)+2i⟨σΣrεrEz,φ⟩Σ,~a3((Ez,Hz),(φ,ψ))=−(∇s×^zHz,∇sφ)+(∇s×^zEz,∇sψ),~a4((Ez,Hz),(φ,ψ))=−d¯¯¯d(εrEz,φ)−d¯¯¯d(Hz,ψ)−i⟨σΣrEz,φ⟩Σ.

### 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 Krylov-space iteration techniques. It is a conformal mapping, defined as follows.

 λ↦aλ+bcλ+d,

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 skew-symmetry

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) Sz=kzMz,

for a complex-valued vector and appropriate complex-valued 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) (aS+bM)z=~kz(cS+dM)z,

where are the Möbius transform parameters. The original eigenvalue can be retrieved via the inverse Möbius transform .

## 4 Validation of weak formulation

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) L(u)[φ]:=a(u,φ)+k2zm(u,φ),

for and where we have introduced the bilinear forms

 (16) {a(u,φ)=∫Ω∇su⋅∇s¯¯¯¯φdx−∫Ωμrεru¯¯¯¯φdx−i∫ΣμrσΣru¯¯¯¯φdox,m(u,φ)=∫Ωu¯¯¯¯φdx+i∫ΣσΣrε−1ru¯¯¯¯φdox.

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) Ez ={AmIm(iksρ)eimθeikzzfor ρR, (18) Hz ={CmIm(iksρ)eimθeikzzfor ρR,

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) det(M−σΣrN) =0,

where

 M :=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝Im(iksR)0−Km(iksR)0kzmIm(iksR)Rk2s−μI′m(iksR)ks−kzmKm(iksR)Rk2sμK′m(iksR)ksϵ1I′m(iksR)kskzmIm(iksR)Rk2s−ϵ2K′m(iksR)ks−kzmKm(iksR)Rk2s0Im(iksR)0−Km(iksR)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠

and

 N :=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝00000000Im(iksR)000−kzmIm(iksR)Rk2sμI′m(iksR)ks00⎞⎟ ⎟ ⎟ ⎟ ⎟⎠.

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 Krylov-space methods, such as the Arnoldi, Lanczos, Krylov-Schur, and conjugate-gradient methods. For our purposes, we will make use of the Krylov-Schur 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 Krylov-Schur 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 radially-varying 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 gradient-index 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) s(ρ)=s0(ρ−0.8R)2(R−0.8R)2,

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