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

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)

{\begin{matrix} \begin{matrix} \nabla \times E & = i ω μ H, \nabla \times H & = - i ω ε E, \end{matrix} \end{matrix}

where $E (x)$ and $H (x)$ denote the electric and magnetic field, respectively, and $ω$ is the temporal frequency. $μ (x, y)$ and $ε (x, y)$ 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 \sim e^{i k_{z} z},

and decompose the fields $F = {E, H}$ and the gradient operator, $\nabla$ , into their longitudinal and transverse parts, whence we obtain

F = F_{s} +^z F_{z}, \nabla = \nabla_{s} + i k_{z}^z,

where the subscript s denotes the transverse direction and $^z$ denotes the unit vector in $z$ -direction. In the strong sense, (1) holds true everywhere except on the points comprising the conducting interface, $Σ$ . The surface conductivity $σ^{Σ} (x, y)$ on the conducting interface $Σ$ gives rise to a jump condition on the tangential part of the magnetic field maier17a . In summary, we obtain

(2)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \begin{matrix} {[ν \times (H_{s} +^z H_{z})]}_{Σ} & = σ^{Σ} (ν \times (E_{s} +^z E_{z})) \times ν, {[ν \cdot μ (H_{s} +^z H_{z})]}_{Σ} & = 0, [ν \times (E_{s} +^z E_{z})]_{Σ} & = 0, [ν \cdot ε (E_{s} +^z E_{z})]_{Σ} & = \frac{1}{i ω} \nabla \cdot (σ^{Σ} E), \end{matrix} \end{matrix}

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 \in Σ .

We also fix the notation

F^{+} := lim α ↘ 0 F (x + α ν), F^{-} := lim α ↘ 0 F (x - α ν), for x \in Σ .

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 \to ˘ x = k_{0} x, \nabla_{s} \to {˘ \nabla}_{s} = \frac{1}{k_{0}} \nabla_{s},$
	$μ→μr=1μ0μ,\e→εr=1\e0ε,σΣ→σΣr=√μ0\e0σΣ$
	$E \to ˘ E = \frac{k_{0}^{2}}{ω μ_{0}} E, H \to ˘ H = k_{0} H, k_{z} \to {˘ k}_{z} = \frac{k_{z}}{k_{0}} .$

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)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \begin{matrix} _{Σ} \cdot τ & = & {[\frac{i}{k_{s}^{2}} (k_{z} \partial_{τ} H_{z} + ε_{r} \partial_{ν} E_{z})]}_{Σ} & = σ_{r}^{Σ} E_{z}, [H_{z}]_{Σ} & = - σ_{r}^{Σ} E_{s} \cdot τ, [μ_{r} H_{s}]_{Σ} \cdot ν & = & {[\frac{i μ_{r}}{k_{s}^{2}} (k_{z} \partial_{ν} H_{z} - ε_{r} \partial_{τ} E_{z})]}_{Σ} & = 0, [E_{s}]_{Σ} \cdot τ & = & {[\frac{i}{k_{s}^{2}} (k_{z} \partial_{τ} E_{z} - μ_{r} \partial_{ν} H_{z})]}_{Σ} & = 0, [E_{z}]_{Σ} & = 0, [ε_{r} E_{s}]_{Σ} \cdot ν & = & {[\frac{i ε_{r}}{k_{s}^{2}} (k_{z} \partial_{ν} E_{z} + μ_{r} \partial_{τ} H_{z})]}_{Σ} & = \frac{1}{i} \nabla_{s} \cdot (σ_{r}^{Σ} E_{s}), \end{matrix} \end{matrix}

where $\partial_{τ}$ and $\partial_{ν}$ denote the derivative in the tangential and the normal direction, respectively; $k_{s} (x, y)^{2} = μ_{r} (x, y) ε_{r} (x, y) - k_{z}^{2}$ is a function in the transverse direction; and where we have used the identities (see Appendix A)

(4)

{\begin{matrix} \begin{matrix} k_{s}^{2} E_{s} & = i (k_{z} \nabla_{s} E_{z} + μ_{r} \nabla_{s} \times H_{z}), k_{s}^{2} H_{s} & = i (k_{z} \nabla_{s} H_{z} - ε_{r} \nabla_{s} \times E_{z}) . \end{matrix} \end{matrix}

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 $Ω \subset R^{n}$ , where $n = 2, 3$ , be a simply connected and bounded domain with Lipschitz-continous piecewise smooth boundary, $\partial Ω$ . 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 $\nabla_{s} \times (μ_{r}^{- 1} \nabla_{s} \times^z E_{z}) = - \nabla_{s} \cdot (μ_{r}^{- 1} \nabla_{s}^z E_{z})$ and $\nabla_{s} \times (μ_{r}^{- 1}^z \times E_{s}) =^z \nabla_{s} \cdot (μ_{r}^{- 1} E_{s})$ , which can be used in conjunction with (27) and (5) to obtain:

(6)

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} \begin{matrix} - \nabla_{s} \cdot (\frac{ε_{r}}{k_{s}^{2}} \nabla_{s} E_{z}) - k_{z} \nabla_{s} \cdot (\frac{1}{k_{s}^{2}} \nabla_{s} \times^z H_{z}) - ε_{r} E_{z} & = 0, - \nabla_{s} \cdot (\frac{μ_{r}}{k_{s}^{2}} \nabla_{s} H_{z}) + k_{z} \nabla_{s} \cdot (\frac{1}{k_{s}^{2}} \nabla_{s} \times^z E_{z}) - μ_{r} H_{z} & = 0. \end{matrix} \end{matrix}

We observe that if the domain is homogeneously filled and isotropic, the curl terms vanish, yielding the familiar decoupled Helmholtz equation for $E_{z}$ and $H_{z}$ . 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 $k_{s}^{4}$ and testing the first equation by $φ$ and the second equation by $ψ$ , we obtain

(8)

(μ_{r} ε_{r}^{2} \nabla_{s} E_{z}, \nabla_{s} φ) + 2 (ε_{r} \nabla_{s} E_{z}, \nabla_{s} ({¯ ¯ ¯ ε}_{r} {¯ ¯ ¯ μ}_{r}) φ) + k_{z} (μ_{r} ε_{r} \nabla_{s} \times^z H_{z}, \nabla_{s} φ) - k_{z}^{2} (ε_{r} \nabla_{s} E_{z}, \nabla_{s} φ) + 2 k_{z} (\nabla_{s} \times^z H_{z}, \nabla_{s} ({¯ ¯ ¯ ε}_{r} {¯ ¯ ¯ μ}_{r}) φ) - k_{z}^{3} (\nabla_{s} \times^z H_{z}, \nabla_{s} φ) - (ε_{r} k_{s}^{4} E_{z}, φ) + (ε_{r} μ_{r}^{2} \nabla_{s} H_{z}, \nabla_{s} ψ) + 2 (μ_{r} \nabla_{s} H_{z}, \nabla_{s} ({¯ ¯ ¯ ε}_{r} {¯ ¯ ¯ μ}_{r}) ψ) - k_{z} (μ_{r} ε_{r} \nabla_{s} \times^z E_{z}, \nabla_{s} ψ) - k_{z}^{2} (μ_{r} \nabla_{s} H_{z}, \nabla_{s} ψ) - 2 k_{z} (\nabla_{s} \times^z E_{z}, \nabla_{s} ({¯ ¯ ¯ ε}_{r} {¯ ¯ ¯ μ}_{r}) ψ) + k_{z}^{3} (\nabla_{s} \times^z E_{z}, \nabla_{s} ψ) - (μ_{r} k_{s}^{4} H_{z}, ψ) - i ⟨ σ_{r}^{Σ} k_{s}^{4} E_{z}, φ ⟩_{Σ} = 0.

This shows the following statement.

Proposition 1.

Provided that $k_{s} \neq 0$ and $[ε_{r}]_{Σ} = [μ_{r}]_{Σ} = 0$ , the nonlinear eigenvalue problem

(N)

Find

k_{z} \in C ∖ {0}

and

E_{z}

H_{z}

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 (k_{z}, (E_{z}, H_{z})) [(φ, ψ)] = 4 \sum l = 0 (k_{z})^{l} a_{l} ((E_{z}, H_{z})) [(φ, ψ)],

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 $X_{h} (Ω) \subset X (Ω)$ be a finite element subspace spanned by Lagrange finite elements $Q_{p}$ .

_{h}

)

4 \sum l = 0 k_{z}^{l} a_{l} ((E_{z}, H_{z}), (φ, ψ)) = 0,

Our goal is to translate (Q $_{h}$ ) 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 $λ \in C$ , and $x \in C^{n}$ , s. t.

4 \sum k = 0 λ^{k} A_{k} x = 0,

where $A_{k} \in C^{n \times n}$ 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 $λ \in C$ , and $z \in C^{4 n}$ s. t.

(11)

⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} A_{3} & 0 & - I_{n} & 0 A_{1} & 0 & 0 & - I_{n} 0 & - I_{n} & 0 & 0 A_{0} & 0 & 0 & 0 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠ z + λ ⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} A_{4} & 0 & 0 & 0 A_{2} & I_{n} & 0 & 0 0 & 0 & I_{n} & 0 0 & 0 & 0 & I_{n} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠ z = 0.

Here, $I_{n}$ denotes the $n \times n$ identity matrix.

With this result at hand, we rewrite (Q $_{h}$ ) as a linear eigenvalue problem:

(LQ

_{h}

)

\begin{matrix} S z + λ M z = 0, \end{matrix}

where

S = ⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} a_{3} & 0 & - I_{n} & 0 a_{1} & 0 & 0 & - I_{n} 0 & - I_{n} & 0 & 0 a_{0} & 0 & 0 & 0 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠, M = ⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} a_{4} & 0 & 0 & 0 a_{2} & I & 0 & 0 0 & 0 & I & 0 0 & 0 & 0 & I \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠, z = ⎛ ⎜ ⎜ ⎜ ⎝ \begin{matrix} z_{1} z_{2} z_{3} z_{4} \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎠ .

Here, by some abuse of notation $a_{i}$ denotes the corresponding matrix formed by the bilinear form $a_{i} (. .)$ given in (10) and by fixing a basis of $X_{h} (Ω)$

. A quick computation shows that the eigenvectors of the original problem (

_{h}

) and of the final linearized problem (LQ

_{h}

) are related as follows.

Proposition 4.

Let $λ \in C$ and $x \in C^{n}$ be an eigenvalue with corresponding eigenvector of (Q $_{h}$ ). Then, $λ$ and $z \in C^{4 n}$ given by

(12)

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ \begin{matrix} z_{1} & = λ x, z_{2} & = λ^{2} (a_{3} + λ a_{4}) x, z_{3} & = λ (a_{3} + λ a_{4}) x, z_{4} & = - a_{0} x, \end{matrix}

is an eigenvalue with corresponding eigenvector of (LQ $_{h}$ ). Conversely, if $λ \in C$ and $z \in C^{4 n}$ is an eigenvalue and eigenvector pair of (LQ $_{h}$ ), then provided that $det (a_{0}) \neq 0$ and $det (λ a_{4} + a_{3}) \neq 0$ , the vector $z$ is characterized by (12) and $λ$ and $x$ are an eigenvalue and eigenvector pair of (Q $_{h}$ ).

3.2 Perfectly Matched Layer

Figure 4: Schematic of the computational domain of a circular waveguide $Ω$ , with boundary $\partial Ω$ , outer normal vector $ν$ , and tangential vector $τ$ . A conducting ciruclar interface is labeled $Σ$ . A perfectly matched layer (PML) is enforced in the shaded region.

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 $~ ρ = ρ ¯ ¯ ¯ d$ , where

d = 1 + i s (r), ¯ ¯ ¯ d = 1 + i / r \int_{ρ}^{r} s (τ) d τ,

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

{^\nabla}_{s} \cdot (k_{s}^{2} (ε_{r} {^\nabla}_{s} {^E}_{z})) - 2 k_{s}^{2} {^\nabla}_{s} \cdot (ε_{r} {^\nabla}_{s} {^E}_{z}) + k_{z} {^\nabla}_{s} \cdot (k_{s}^{2} {^\nabla}_{s} \times {^H}_{z}) - 2 k_{z} k_{s}^{2} {^\nabla}_{s} \cdot ({^\nabla}_{s} \times {^H}_{z}) + {^\nabla}_{s} \cdot (k_{s}^{2} (μ_{r} {^\nabla}_{s} {^H}_{z})) - 2 k_{s}^{2} {^\nabla}_{s} \cdot (μ_{r} {^\nabla}_{s} {^H}_{z}) + k_{z} {^\nabla}_{s} \cdot (k_{s}^{2} {^\nabla}_{s} \times {^E}_{z}) + 2 k_{z} k_{s}^{2} {^\nabla}_{s} \cdot ({^\nabla}_{s} \times {^E}_{z}) - ε_{r} k_{s}^{4} {^E}_{z} - μ_{r} k_{s}^{4} {^H}_{z} = 0.

This can be rewritten as

\frac{1}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (k_{s}^{2} ε_{r} A \nabla_{s} E_{z}) - \frac{2 k_{s}^{2}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (ε_{r} A \nabla_{s} E_{z}) + \frac{k_{z}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (k_{s}^{2} \nabla_{s} \times H_{z}) - \frac{2 k_{z} k_{s}^{2}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (\nabla_{s} \times H_{z}) + \frac{1}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (k_{s}^{2} ε_{r} A \nabla_{s} E_{z}) - \frac{2 k_{s}^{2}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (ε_{r} A \nabla_{s} E_{z}) + \frac{k_{z}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (k_{s}^{2} \nabla_{s} \times H_{z}) + \frac{2 k_{z} k_{s}^{2}}{d ¯ ¯ ¯ d} \nabla_{s} \cdot (\nabla_{s} \times E_{z}) + ε_{r} k_{s}^{4} E_{z} + μ_{r} k_{s}^{4} H_{z} = 0,

where $A = T_{e_{x}, e_{r}}^{- 1} diag (¯ ¯ ¯ d / d, d / ¯ ¯ ¯ d) T_{e_{x}, e_{r}}$ , and $T_{e_{x}, e_{r}}$ is the rotation matrix that rotates $e_{r}$ onto $e_{x}$ . We enforce the condition that the material parameters are constant outside the PML, i.e., $ε_{r}$ and $μ_{r}$ 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, ${~ a}_{i}$ , are

(13)

\begin{matrix} ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ \begin{matrix} {~ a}_{0} ((E_{z}, H_{z}), (φ, ψ)) = (ε_{r}^{2} A \nabla_{s} E_{z}, \nabla_{s} φ) + 2 (ε_{r} A \nabla_{s} E_{z}, (\nabla_{s} {¯ ¯ ¯ ε}_{r}) φ) - d ¯ ¯ ¯ d (ε_{r}^{3} E_{z}, φ) + (ε_{r} A \nabla_{s} H_{z}, \nabla_{s} ψ) + 2 (A \nabla_{s} H_{z}, (\nabla_{s} {¯ ¯ ¯ ε}_{r}) ψ) - d ¯ ¯ ¯ d (ε_{r}^{2} H_{z}, ψ) - i ⟨ σ_{r}^{Σ} ε_{r}^{2} E_{z}, φ ⟩_{Σ}, {~ a}_{1} ((E_{z}, H_{z}), (φ, ψ)) = (ε_{r} \nabla_{s} \times^z H_{z}, \nabla_{s} φ) + 2 (\nabla_{s} \times^z H_{z}, (\nabla_{s} {¯ ¯ ¯ ε}_{r}) φ) - (ε_{r} \nabla_{s} \times^z E_{z}, \nabla_{s} ψ) - 2 (\nabla_{s} \times^z E_{z}, (\nabla_{s} {¯ ¯ ¯ ε}_{r}) ψ) {~ a}_{2} ((E_{z}, H_{z}), (φ, ψ)) = - (ε_{r} A \nabla_{s} E_{z}, \nabla_{s} φ) + 2 d ¯ ¯ ¯ d (ε_{r}^{2} E_{z}, φ) - (A \nabla_{s} H_{z}, \nabla_{s} ψ) + 2 d ¯ ¯ ¯ d (ε_{r} H_{z}, ψ) + 2 i ⟨ σ_{r}^{Σ} ε_{r} E_{z}, φ ⟩_{Σ}, {~ a}_{3} ((E_{z}, H_{z}), (φ, ψ)) = - (\nabla_{s} \times^z H_{z}, \nabla_{s} φ) + (\nabla_{s} \times^z E_{z}, \nabla_{s} ψ), {~ a}_{4} ((E_{z}, H_{z}), (φ, ψ)) = - d ¯ ¯ ¯ d (ε_{r} E_{z}, φ) - d ¯ ¯ ¯ d (H_{z}, ψ) - i ⟨ σ_{r}^{Σ} E_{z}, φ ⟩_{Σ} . \end{matrix} \end{matrix}

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 $0$ (or maximal), so that they can be selectively computed with conventional Krylov-space iteration techniques. It is a conformal mapping, defined as follows.

λ \mapsto \frac{a λ + b}{c λ + d},

where $λ, a, b, c, d \in C$

. 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)

S z = k_{z} M z,

for a complex-valued vector $z$ and appropriate complex-valued matrices $S$ and $M$ . The implementation of the PML discussed in 3.2 necessitates changes to the definition of $a_{1}, a_{2}, a_{3}$ , and $a_{4}$ . The idea is to use the Möbius transform to map points near the origin to target values

(15)

(a S + b M) z = {~ k}_{z} (c S + d M) z,

where $a, b, c, d$ are the Möbius transform parameters. The original eigenvalue can be retrieved via the inverse Möbius transform $k_{z} = (- b + d {~ k}_{z}) / (a - c {~ k}_{z})$ .

4 Validation of weak formulation

	Eigenvalues from (L)	Eigenvalues from (9)	Eigenvalues from (19)
Mode	$k_{z}$	$k_{z}$	$k_{z}$
1	$9.447 \pm 0.090468 i$	$9.447 \pm 0.090467 i$	$9.447 \pm 0.090467 i$
2	$13.00 \pm 0.090641 i$	$13.00 \pm 0.090640 i$	$13.00 \pm 0.090641 i$
3	$16.17 \pm 0.099812 i$	$16.17 \pm 0.099812 i$	$16.17 \pm 0.099813 i$

Table 1: Validation of computed eigenvalues form the linear problem (16), the quartic problem (9), and the analytical approach (19) for a single-layer waveguide. Note that the values obtained are in agreement with confidence level of less than 1%. Parameters used are

R = 2.0, σ_{r}^{Σ} = 0.002 + 0.20 i, ε_{r} = μ_{r} \equiv 1,

and

ρ_{i} = 0.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 $u \in H^{1} (Ω, C)$ s. t.

(L)

L (u) [φ] := a (u, φ) + k_{z}^{2} m (u, φ),

for $φ \in H^{1} (Ω, C)$ and where we have introduced the bilinear forms

(16)

\begin{matrix} {\begin{matrix} a (u, φ) = \int_{Ω} \nabla_{s} u \cdot \nabla_{s} ¯ ¯¯ ¯ φ d x - \int_{Ω} μ_{r} ε_{r} u ¯ ¯¯ ¯ φ d x - i \int_{Σ} μ_{r} σ_{r}^{Σ} u ¯ ¯¯ ¯ φ d o_{x}, m (u, φ) = \int_{Ω} u ¯ ¯¯ ¯ φ d x + i \int_{Σ} σ_{r}^{Σ} ε_{r}^{- 1} u ¯ ¯¯ ¯ φ d o_{x} . \end{matrix} \end{matrix}

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 $0$ and radius $R$ , the analytical solution takes the following form.

(17)		$E_{z}$	$= {\begin{matrix} \begin{matrix} A_{m} I_{m} (i k_{s} ρ) e^{i m θ} e^{i k_{z} z} & for ρ < R, B_{m} K_{m} (i k_{s} ρ) e^{i m θ} e^{i k_{z} z} & for ρ > R, \end{matrix} \end{matrix}$
(18)		$H_{z}$	$= {\begin{matrix} \begin{matrix} C_{m} I_{m} (i k_{s} ρ) e^{i m θ} e^{i k_{z} z} & for ρ < R, D_{m} K_{m} (i k_{s} ρ) e^{i m θ} e^{i k_{z} z} & for ρ > R, \end{matrix} \end{matrix}$

where $I_{m}$ and $K_{m}$ denote the modified Bessel functions of the first and second kind, respectively, and $A_{m}, B_{m}, C_{m}$ , and $D_{m}$ 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 $R$ , 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, $k_{z}$ :

(19)

det (M - σ_{r}^{Σ} N)

= 0,

where

M

:= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} I_{m} (i k_{s} R) & 0 & - K_{m} (i k_{s} R) & 0 \frac{k_{z} m I_{m} (i k_{s} R)}{R k_{s}^{2}} & - \frac{μ I_{m}^{'} (i k_{s} R)}{k_{s}} & - \frac{k_{z} m K_{m} (i k_{s} R)}{R k_{s}^{2}} & \frac{μ K_{m}^{'} (i k_{s} R)}{k_{s}} \frac{ϵ_{1} I_{m}^{'} (i k_{s} R)}{k_{s}} & \frac{k_{z} m I_{m} (i k_{s} R)}{R k_{s}^{2}} & - \frac{ϵ_{2} K_{m}^{'} (i k_{s} R)}{k_{s}} & - \frac{k_{z} m K_{m} (i k_{s} R)}{R k_{s}^{2}} 0 & I_{m} (i k_{s} R) & 0 & - K_{m} (i k_{s} R) \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

and

N

:= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 0 & 0 & 0 & 0 0 & 0 & 0 & 0 I_{m} (i k_{s} R) & 0 & 0 & 0 - \frac{k_{z} m I_{m} (i k_{s} R)}{R k_{s}^{2}} & \frac{μ I_{m}^{'} (i k_{s} R)}{k_{s}} & 0 & 0 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ .

We numerically solve the zeroes of the determinant of the difference of the above matrices via a root finding algorithm for modal orders $m = 1, 2$ and $3$ . 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 $ε_{r} = μ_{r} \equiv 1$ . A conducting interface is coated on the boundary of the interior circle with radius $ρ_{i} = 0.3$ , and choose $R = 2.0$ . The surface conductivity is set to $σ_{r}^{Σ} = 0.002 + 0.20 i$ .

The analytic computation of the propagation constant $k_{z}$ requires finding the complex roots of the determinant of the $4 \times 4$ matrix. We will consider modal orders of $m = 1, 2$ and $3$ 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, $∥ E_{z} ∥_{2}$ , 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, $ε_{r} (ρ, θ)$ . 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 $_{h}$ ) (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 $_{h}$ ) 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

	$s_{o} = 1.5, R = 1.0, k = 4, DoF = 19138$			$s_{o} = 1.5, R = 1.0, k = 6, DoF = 303874$
Mode	Re $k_{z}$	Im $k_{z}$	Re $k_{z} / Im k_{z}$	Re $k_{z}$	Im $k_{z}$	Re $k_{z} / Im k_{z}$
1	30.5400	0.17646	173.070	30.4160	0.17540	173.409
2	40.5520	0.21898	185.186	40.5043	0.21881	185.112

(a) Vary the number of refinements,

k

	$s_{o} = 1.5, R = 1.0, k = 5, DoF = 76162$
Mode	Re $k_{z}$	Im $k_{z}$	Re $k_{z} / Im k_{z}$
1	30.4160	0.17582	172.995
2	40.5043	0.21881	185.112

(b) The control case for the parameter studies.

	$s_{o} = 1.5, R = 0.5, k = 5, DoF = 76162$			$s_{o} = 1.5, R = 2.0, k = 5, DoF = 76162$
Mode	Re $k_{z}$	Im $k_{z}$	Re $k_{z} / Im k_{z}$	Re $k_{z}$	Im $k_{z}$	Re $k_{z} / Im k_{z}$
1	30.4242	0.17528	173.575	30.8567	0.17331	178.043
2	40.4905	0.21876	185.091	40.5043	0.21881	185.112

R

Table 5: Parameter study with permittivity function

ε_{r} (x) = 3 \cdot χ_{| x | < ρ_{i}} + 1 \cdot χ_{ρ_{i} < | x | < R}

μ_{r} (x) \equiv 1

and

σ_{r}^{Σ} = 0.002 + 0.20 i

. (b) is the control discretization parameters. (a) differs from (b) in the number of refinements, and (c) differs from (b) in the size of the domain.

The computational domain, $Ω$ , is chosen to be the circle with radius 1. A spherical PML is enforced for $ρ > 0.8$ . The surface conductivity $σ_{r}^{Σ} = 0.002 + 0.20 i$ is chosen that is within a realistic parameter range maier17a . Following monk03 , the nonnegative scaling function $s (ρ)$ is chosen to be

(20)

s (ρ) = s_{0} \frac{(ρ - 0.8 R)^{2}}{(R - 0.8 R)^{2}},

where we set the free parameter $s_{0}$ to be $s_{0} = 1.5$ 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 $s_{0}$ , the number of initial refinements $k$ , and domain sizes $R$ . A spectral transformation is carried out in the form of the Möbius transformation, with the parameters chosen to be

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

Authors

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

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

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

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

Computation of Tight Enclosures for Laplacian Eigenvalues

Approximation of PDE eigenvalue problems involving parameter dependent matrices

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