PDE eigenvalue iterations with applications in two-dimensional photonic crystals

TL;DR

This paper develops iterative methods for solving PDE eigenvalue problems in 2D photonic crystals, including nonlinear cases, demonstrating convergence and efficiency improvements through adaptive refinement and Newton-type iterations.

Contribution

It extends classical eigenvalue algorithms to infinite-dimensional PDE problems and introduces a Newton method for nonlinear eigenproblems in photonic crystal modeling.

Findings

Inverse power method converges on operator level.

Adaptive mesh refinement accelerates computations.

Newton-type iteration achieves local quadratic convergence.

Abstract

We consider PDE eigenvalue problems as they occur in two-dimensional photonic crystal modeling. If the permittivity of the material is frequency-dependent, then the eigenvalue problem becomes nonlinear. In the lossless case, linearization techniques allow an equivalent reformulation as an extended but linear and Hermitian eigenvalue problem, which satisfies a Garding inequality. For this, known iterative schemes for the matrix case such as the inverse power or the Arnoldi method are extended to the infinite-dimensional case. We prove convergence of the inverse power method on operator level and consider its combination with adaptive mesh refinement, leading to substantial computational speed-ups. For more general photonic crystals, which are described by the Drude-Lorentz model, we propose the direct application of a Newton-type iteration. Assuming some a priori knowledge on the eigenpair of interest, we prove local quadratic convergence of the method. Finally, numerical experiments confirm the theoretical findings of the paper.

Figures

Tables

Table 1. Table 6.1 . Runtime comparison (in seconds) until the residual reaches the tolerance tol= 10 − 10 superscript 10 10 10^{-10} in the first numerical experiment.

on fine mesh	$3$ steps per mesh	$5$ steps per mesh	$7$ steps per mesh
129.136	72.396	46.585	42.308

Table 2. Table 6.2 . Runtime comparison (in seconds) until the residual reaches the tolerance tol= 10 − 6 superscript 10 6 10^{-6} in the second numerical experiment.

on fine mesh	$3$ steps per mesh	$5$ steps per mesh	$7$ steps per mesh
100.56	71.96	51.39	34.49

Table 3. Table 6.3 . Material parameters representing 66% porous silicon, used in the third numerical experiment, cf.  [ EKE12 ] .

$\ell$	1	2	3	4	5	6	7
${\xi }_{\ell }^2$	416.6166	352.7054	-339.9124	492.5687	-19.6143	-527.5597	98.0101
${\eta }_{\ell }^2$	92.1086	71.6269	71.4552	227.8301	47.4923	93.5605	121.3762
${\gamma }_{\ell }$	2.7820	0.9597	0.9500	13.1508	9.2697	3.2624	2.2712

Equations

$$a(u,v)=\lambda (u,v).$$

$$a(u,u)\ge \alpha \Vert u{\Vert }_{\mathcal{V}}^2-\beta \Vert u{\Vert }^2$$

$$\Vert u{\Vert }_{\beta }:=a_{\beta }(u,u{)}^{1/2}\ge \sqrt{\alpha }\Vert u{\Vert }_{\mathcal{V}},$$

$$|\mathbf{a}\cdot \overline{\mathbf{b}}+\mathbf{b}\cdot \overline{\mathbf{a}}|\le \frac{1}{\delta }|\mathbf{a}{|}^2+\delta |\mathbf{b}{|}^2,$$

$$0$$

$$a(u,v)=\lambda (u,v).$$

$$a(u,v):={\int }_{\Omega }\nabla u(x)\cdot \overline{\nabla v(x)}+W(x)u(x)\overline{v(x)}\text{d}x$$

$$1\le \epsilon (x,\omega )\le {\epsilon }_{\max }<\infty$$

$$-{\nabla }_{\mathbf{k}}\cdot {\nabla }_{\mathbf{k}}{u}_{\mathbf{k}}(x)={\omega }^2{\mu }_0{\epsilon }_0\epsilon (x)u_{\mathbf{k}}(x)$$

$$a(u,v):={\int }_{\Omega }{\nabla }_{\mathbf{k}}u(x)\cdot \overline{{\nabla }_{\mathbf{k}}v(x)}\text{d}x$$

$$(u,v):={\int }_{\Omega }\epsilon (x)u(x)\overline{v(x)}\text{d}x.$$

$$a(u,v)\leq\int_{\Omega}|\nabla_{\mathbf{k}}u||\nabla_{\mathbf{k}}v|\,\text{d}x\leq\Big{(}\int_{\Omega}|\nabla_{\mathbf{k}}u|^{2}\,\text{d}x\Big{)}^{1/2}\Big{(}\int_{\Omega}|\nabla_{\mathbf{k}}v|^{2}\,\text{d}x\Big{)}^{1/2}.$$

$${\int }_{\Omega }|{\nabla }_{\mathbf{k}}u{|}^2\text{d}x$$

$$\le {\int }_{\Omega }|\nabla u{|}^2+|\mathbf{k}{|}^2|u{|}^2+|\nabla u{|}^2+|\mathbf{k}{|}^2|u{|}^2\text{d}x\le 2\max \{1,|\mathbf{k}{|}^2\}\Vert u{\Vert }_{\mathcal{V}}^2,$$

$$a(u,u)={\int }_{\Omega }|\nabla u{|}^2\text{d}x\ge (1-\beta )\Vert \nabla u{\Vert }_{L^2(\Omega )}^2+\beta \Vert u{\Vert }_{\mathcal{V}}^2-\beta {\int }_{\Omega }\epsilon u\overline{u}\text{d}x\ge \beta \Vert u{\Vert }_{\mathcal{V}}^2-\beta \Vert u{\Vert }^2,$$

$$a(u,u)={\int }_{\Omega }|{\nabla }_{\mathbf{k}}u{|}^2\text{d}x$$

$$\ge {\int }_{\Omega }|\nabla u{|}^2+|\mathbf{k}{|}^2|u{|}^2-\frac{|\nabla u{|}^2}{\delta }-\delta |\mathbf{k}{|}^2|u{|}^2\text{d}x.$$

$$-{\nabla }_{\mathbf{k}}\cdot \frac{1}{\epsilon (x)}{\nabla }_{\mathbf{k}}{u}_{\mathbf{k}}(x)=\lambda u_{\mathbf{k}}(x)$$

$$a(u,v):={\int }_{\Omega }\frac{1}{\epsilon (x)}{\nabla }_{\mathbf{k}}u\cdot \overline{{\nabla }_{\mathbf{k}}v}\text{d}x.$$

$$⟨\mathcal{A}u,v⟩:=a(u,v).$$

$$\mathcal{A}u=\lambda \mathcal{I}u\text{in }{\mathcal{V}}^{\ast }.$$

$${\mathcal{A}}_{\beta }u:=(\mathcal{A}+\beta \mathcal{I})u=(\lambda +\beta )\mathcal{I}u=:\mu \mathcal{I}u\text{in }{\mathcal{V}}^{\ast }.$$

$${\mathcal{A}}_{\beta }{u}^j={\mu }^{j-1}\mathcal{I}{\tilde{u}}^{j-1}\text{in }{\mathcal{V}}^{\ast }.$$

$${\mu }^j:=\frac{⟨{\mathcal{A}}_{\beta }{u}^j,u^j⟩}{⟨\mathcal{I}{u}^j,u^j⟩}=\frac{a_{\beta }(u^j,u^j)}{(u^j,u^j)}=\frac{\Vert u^j{\Vert }_{\beta }^2}{\Vert u^j{\Vert }^2}>0.$$

$$\Vert u^j{\Vert }_{\beta }^2\le {\mu }^{j-1}\Vert u^j\Vert ,{\mu }^{j-1}={\mu }^{j-1}\Vert {\tilde{u}}^{j-1}{\Vert }^2\le \sqrt{{\mu }^{j-1}}\Vert u^j{\Vert }_{\beta }.$$

$$\sqrt{{\mu }^j}=\frac{\Vert u^j{\Vert }_{\beta }}{\Vert u^j\Vert }\le \sqrt{{\mu }^{j-1}}.$$

$$\tilde{u}^{j}\rightharpoonup{\tilde{u}}^{\scalebox{0.5}{$\bigstar$}}\quad\text{in }\mathcal{V},\qquad\tilde{u}^{j}\rightarrow{\tilde{u}}^{\scalebox{0.5}{$\bigstar$}}\quad\text{in }\mathcal{H}.$$

$$\Vert u^j{\Vert }_{\beta }^2\le {\mu }^{j-1}\Vert u^j\Vert \le {\mu }^0\Vert u^j\Vert \le C_{\mathcal{V}\hookrightarrow \mathcal{H}}{\mu }^0\Vert u^j{\Vert }_{\mathcal{V}}\le C_{\mathcal{V}\hookrightarrow \mathcal{H}}{\mu }^0{\alpha }^{-\frac{1}{2}}\Vert u^j{\Vert }_{\beta }.$$

$$u^{j}\rightharpoonup u^{\scalebox{0.5}{$\bigstar$}}\quad\text{in }\mathcal{V},\qquad u^{j}\rightarrow u^{\scalebox{0.5}{$\bigstar$}}\quad\text{in }\mathcal{H}.$$

$${\mu }^j=\frac{a_{\beta }(u^j,u^j)}{\Vert u^j{\Vert }^2}={\mu }^{j-1}\frac{({\tilde{u}}^{j-1},u^j)}{\Vert u^j{\Vert }^2}=\frac{{\mu }^{j-1}}{\Vert u^j\Vert }({\tilde{u}}^{j-1},{\tilde{u}}^j)$$

$$\mu^{\scalebox{0.5}{$\bigstar$}}=\frac{\mu^{\scalebox{0.5}{$\bigstar$}}}{\|u^{\scalebox{0.5}{$\bigstar$}}\|}\|{\tilde{u}}^{\scalebox{0.5}{$\bigstar$}}\|^{2}=\frac{\mu^{\scalebox{0.5}{$\bigstar$}}}{\|u^{\scalebox{0.5}{$\bigstar$}}\|},$$

Full text

∗ Research funded by the Einstein Foundation Berlin in the frame of Einstein Center for Mathematics Berlin ECMath via project OT10 Model Reduction for Nonlinear Parameter-Dependent Eigenvalue Problems in Photonic Crystals and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – The Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID: 390685689).

**Key words. nonlinear eigenvalue problem, photonic crystals, inverse power method, Newton iteration

AMS subject classifications. 65N25, 65J10, 65F15**

1. Introduction

Eigenvalue problems including partial differential equations (PDE) appear in several applications such as structural mechanics [BW73], fluid-solid structures [Vos03], or the simulation of Bose-Einstein condensates [PS03]. In general, such problems are considered in order to optimize certain properties or parameters of the underlying dynamical system [MV04]. In this paper, we focus on applications as they appear in the modeling of photonic crystals [Joh87, Kuc01]. These are special composite materials with a periodic structure that affect the propagation of electromagnetic waves and thus, can be used for trapping and guiding light. As these crystals can be designed and manufactured for industrial applications, the aim is to find so-called photonic band-gaps, which prevent light within a specified frequency range from propagating [JJWM08, Joh12]. Direct applications areas are optical fibers [GH14, Ch. 5], medical technologies with laser guides for cancer surgeries [Tsa12], and thin film solar cells [DJ12].

The corresponding mathematical model is given by a sequence of nonlinear PDE eigenvalue problems based on the Maxwell equations [SP05, DLP*+*11]. An important role is played by the electric permittivity $\varepsilon$, which is periodic in space and characterizes certain properties of the crystal. If $\varepsilon$ is independent of the frequency, then we obtain a linear eigenvalue problem. In more realistic models, however, the permittivity is approximated by a rational function, which introduces the nonlinearity into the eigenvalue problem.

Numerical methods for computing the spectrum of such materials have been studied intensively. This includes adaptive finite element methods [BKS*+*06, GG12], Newton-type methods [Kre09, HLM16], and linearization techniques [SB11, EKE12]. In the latter case, a spatial discretization is assumed and yields then a linear but extended eigenvalue problem, for which well-known iteration schemes can be applied. The combination of mesh refinement and (inexact) eigenvalue iteration methods has been considered in [MM11, Mie11].

Corresponding iterative methods for the operator case have, so far, not received much attention in the literature. In the first part of this paper, we focus on linear and Hermitian PDE eigenvalue problems in the weak formulation. This corresponds to the case of a frequency-independent permittivity. Convergence of the (inverse) power method for compact operators mapping from a Hilbert space $\mathcal{H}$ to $\mathcal{H}$ was already shown in [ESL95]. In this setting, the proof basically follows the same lines as in the finite-dimensional case. General bounded operators were considered in [EE07] but only together with a power iteration based on the exact (und thus unknown) eigenvalue. Considering the weak formulation, which is more natural in view of spatial discretization methods, we are in a different setting. Nevertheless, the power method converges if an appropriate scaling is included. For the $p$-Laplacian eigenvalue problem, this was shown in [Boz16]. Note that the proven convergence on operator level has the advantage that it is mesh- and basis-independent, cf. [AHP19] for a similar approach applied to nonlinear eigenvector problems. Thus, dependencies on a particular discretization do not play a role. Further, the approach allows to choose the spatial discretization in each iteration step independently.

The second part of the paper focuses on the nonlinear case as it appears in two-dimensional photonic crystal modeling. Here we consider two different paths of either linearizing the problem or applying directly a Newton iteration. In the first case, we adapt the techniques introduced in [SB11, EKE12] in combination with an inverse power iteration applied to the resulting linear problem. Due to the linearization, certain compactness properties get lost, which calls for a novel convergence analysis for the inverse power method. The second strategy translates the local convergence of Newton’s method from [Sch08] to the operator case. An analogous method for infinite-dimensional eigenvalue problems was developed in [AR68], and its local convergence was proven for Fredholm operators with index [math] mapping from a Hilbert space $\mathcal{H}$ to $\mathcal{H}$. We show a similar result when the operators arise from the here considered weak formulation.

The paper is structured as follows. In Section 2 we introduce the problem setting, i.e., the linear PDE eigenvalue problem in its weak and operator formulation. Here we gather all the assumptions on the spaces and included operators. In particular, we assume an underlying Gelfand triple with a compact embedding. Section 3 then considers several iteration schemes including the inverse power method and the Arnoldi method. Two-dimensional photonic crystals with frequency-dependent permittivity are then topic of Sections 4 and 5. First, we consider a special Hermitian case. For this, we apply a linearization and the inverse power method. More realistic models are then discussed in Section 5, for which we prove the local convergence of Newton’s method. Finally, we present three numerical examples in Section 6.

2. Preliminaries

As described in the introduction, we consider the weak formulation of a PDE eigenvalue problem. Given the sesquilinear forms $a\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$ and $(\cdot\,,\cdot)\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}^{+}$, we search for a non-trivial pair $(u,\lambda)\in\mathcal{V}\times\mathbb{C}$ such that for all test functions $v\in\mathcal{V}$ it holds that

[TABLE]

More precisely, considering Hermitian eigenvalue problems, we are interested in the eigenpair corresponding to the smallest eigenvalue. In the following, we gather assumptions on the space $\mathcal{V}$ and the included sesquilinear forms. Afterwards we discuss well-known PDE eigenvalue problems, which fit into the given framework.

2.1. General setting

We start with general assumptions on the involved function spaces.

Assumption 2.1* (Function spaces).*

We assume $\mathcal{V}$ to be a complex, separable, and reflexive Banach space. Furthermore, we assume the existence of a complex and separable Hilbert space $\mathcal{H}$ (the pivot space) such that $\mathcal{V}$, $\mathcal{H}$, $\mathcal{V}^{*}$ form a Gelfand triple, cf. [Zei90, Ch. 23.4] and [Bre10, Ch. 11.4]. This means, in particular, that the embedding $i_{\mathcal{V}\hookrightarrow\mathcal{H}}\colon\mathcal{V}\hookrightarrow\mathcal{H}$ is continuous and dense [Wlo87, Ch. 17.1]. The continuity constant is denoted by $C_{\mathcal{V}\hookrightarrow\mathcal{H}}$.

Assumption 2.2* (Compactness).*

The embedding $\mathcal{V}\hookrightarrow\mathcal{H}$ is compact.

With the pivot space $\mathcal{H}$ in hand, we assume that the sesquilinear form $(\cdot\,,\cdot)$ in (2.1) is also defined for functions in $\mathcal{H}$. We even assume that this defines the inner product in the Hilbert space $\mathcal{H}$ and set $\|\cdot\|:=\|\cdot\|_{\mathcal{H}}=(\cdot\,,\cdot)^{1/2}$. Further, $j_{\mathcal{H}}\colon\mathcal{H}\to\mathcal{H}^{*}$ denotes the Riesz isomorphism. The norm in the space $\mathcal{V}$ is denoted by $\|\cdot\|_{\mathcal{V}}$. For the sesquilinear form $a$ we consider the following assumptions.

Assumption 2.3* (Sesquilinear form).*

The sesquilinear form $a\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$ is assumed to be continuous and Hermitian such that $a(u,u)\in\mathbb{R}$ for all $u\in\mathcal{V}$. Furthermore, $a$ satisfies a Gårding inequality, i.e.,

[TABLE]

for real constants $\alpha>0$, $\beta\geq 0$ and all $u,v\in\mathcal{V}$, see e.g. [QV94].

Remark 2.4*.*

The previous Assumption 2.3 implies that $a_{\beta}(u,v):=a(u,v)+\beta\,(u,v)$ is $\mathcal{V}$-coercive and thus, defines an inner product in $\mathcal{V}$. As a result, $\mathcal{V}$ is actually a Hilbert space and the corresponding norm

[TABLE]

the so-called energy norm, is equivalent to $\|\cdot\|_{\mathcal{V}}$. Further, all eigenvalues $\lambda$ of (2.1) satisfy $a_{\beta}(u,u)=(\lambda+\beta)(u,u)$, which implies $\mathbb{R}\ni(\lambda+\beta)>0$ and thus $\mathbb{R}\ni\lambda>-\beta$.

In order to be well-posed, an eigenvalue problem of the form (2.1) requires boundary conditions. Throughout this paper, we assume that these conditions are included in the space $\mathcal{V}$, cf. the examples in the following subsection. We close this preliminary part with the proof of Young’s inequality in the specific case of complex vectors.

Lemma 2.5** (Young’s inequality).**

Consider $\mathbf{a},\mathbf{b}\in\mathbb{C}^{2}$. Then, for every $\delta>0$ we have an estimate of the form

[TABLE]

where $\cdot$ denotes the real dot product.

Proof.

For any two vectors $\mathbf{c},\mathbf{d}\in\mathbb{C}^{2}$, the following estimates hold:

[TABLE]

As a consequence, we have $|\mathbf{c}\cdot\overline{\mathbf{d}}+\mathbf{d}\cdot\overline{\mathbf{c}}|\leq|\mathbf{c}|^{2}+|\mathbf{d}|^{2}$. The claim then follows by setting $\mathbf{c}=\mathbf{a}/\sqrt{\delta}$, and $\mathbf{d}=\mathbf{b}\sqrt{\delta}$. ∎

2.2. Examples

We present a couple of well-known examples, which fit into the given framework if formulated in the weak setting.

Example 2.6** (Laplace eigenvalue problem).**

Consider the eigenvalue problem $-\Delta u=\lambda u$ in a bounded domain $\Omega$ with homogeneous Dirichlet boundary conditions. For this, we set $\mathcal{V}:=H^{1}_{0}(\Omega)$ with $\|\cdot\|_{\mathcal{V}}:=\|\nabla\cdot\|_{L^{2}(\Omega)}$ and $\mathcal{H}:=L^{2}(\Omega)$ with the standard $L^{2}$ inner product. The corresponding sesquilinear form reads $a(u,v):=\int_{\Omega}\nabla u\cdot\overline{\nabla v}\,\text{d}x$. Note that this implies $a(u,u)=\|u\|_{\mathcal{V}}^{2}$ and thus, $\alpha=1$ and $\beta=0$. The weak form of the Laplace eigenvalue problem then reads: find a non-trivial pair $(u,\lambda)\in\mathcal{V}\times\mathbb{C}$ such that for all $v\in\mathcal{V}$ it holds that

[TABLE]

Example 2.7** (Schrödinger eigenvalue problem).**

The computation of the ground state of the linear Schrödinger operator leads to the sesquilinear form

[TABLE]

with a real-valued potential $W\in L^{\infty}(\Omega)$. For homogeneous Dirichlet boundary conditions this leads to the same spaces $\mathcal{V}$ and $\mathcal{H}$ as in Example 2.6. Further, $a$ satisfies the Gårding inequality with $\alpha=1$ and \beta=\max\big{\{}0,-\inf_{x\in\Omega}W(x)\big{\}}. For periodic boundary conditions one has to replace the space $\mathcal{V}$ accordingly.

In this paper, we focus on applications with photonic crystals. The dynamics of the electromagnetic fields inside such a crystal can be modelled by the Maxwell equations in the whole domain $\mathbb{R}^{d}$, cf. [DLP*+*11, Ch. 1]. These equations combine the magnitudes of the time-harmonic electric and magnetic fields $E$, $H$ and the frequency $\omega$, which takes the role of an eigenvalue.

One crucial parameter within the equations is the relative electric permittivity of the materials inside the crystal. We assume the relative permittivity $\varepsilon$ to be piecewise-constant and periodic in space as well as bounded in the sense that

[TABLE]

for all $x\in\mathbb{R}^{d}$ and $\omega$ in some frequency domain of interest. For the applications in mind, where $\varepsilon$ is given as a rational function, this means that we consider $\omega$ bounded away from the poles. In the two-dimensional case, i.e., when $\varepsilon$ is periodic within a two-dimensional plane and constant along the direction orthogonal to this plane, the Maxwell eigenvalue problem decouples into so-called transverse magnetic (TM) and transverse electric (TE) modes. Thanks to the periodicity of $\varepsilon$, which implies a discrete translational symmetry in the system, a Floquet transformation can be applied to reduce the problem posed in $\mathbb{R}^{2}$ to a family of problems on a bounded domain $\Omega$ called the Wigner–Seitz cell of the crystal lattice, see e.g. [Kuc01, DLP*+*11]. Note that the function $\varepsilon$ is here a unitless quantity expressed relative to $\varepsilon_{0}$, the vacuum permittivity. Within this paper, we assume the magnetic permeability of the crystal to be constant and equal to that of vacuum denoted by $\mu_{0}$.

Example 2.8** (TM mode).**

In the two-dimensional setting we consider the TM mode with a real-valued frequency-independent function $\varepsilon(x)$. The resulting PDE eigenvalue problem describes the third component of the electric field $E_{3}$, from which one can directly compute the components $H_{1}$ and $H_{2}$. Let ${\mathbf{k}}$ be a fixed wave vector in the so-called irreducible Brillouin zone $\mathcal{K}\subset\mathbb{R}^{2}$, cf. [DLP*+*11, Ch. 1], $u_{\mathbf{k}}$ be the Floquet transform of $E_{3}(x)$ at ${\mathbf{k}}$, and $\nabla_{\mathbf{k}}\coloneqq\nabla+i{\mathbf{k}}$ denote the shifted gradient. Then, $u_{\mathbf{k}}$ satisfies the eigenvalue problem

[TABLE]

for all $x\in\Omega$. For the sake of conciseness, we use in the following by abuse of notation a scaled frequency defined by $\omega\to\omega/\sqrt{\mu_{0}\varepsilon_{0}}$. For the weak formulation of the eigenvalue problem we then define $\lambda\coloneqq\omega^{2}$. Including periodic boundary conditions, we set $\mathcal{V}=H^{1}_{\text{per}}(\Omega)$ with the standard $H^{1}$-norm and $\mathcal{H}=L^{2}(\Omega)$. Note that $\mathcal{V}$ is densely embedded in $\mathcal{H}$ and thus, Assumption 2.1 is satisfied, cf. [Bre10, Ch. 4.4]. The sesquilinear form $a$ and the (weighted) inner product in $\mathcal{H}$ then read

[TABLE]

and

[TABLE]

Due to the boundedness of $\varepsilon$, the sesquilinear form $(\cdot,\cdot)$ defines an inner product in $\mathcal{H}$. The following lemma shows that $a$ satisfies Assumption 2.3.

Lemma 2.9**.**

For a fixed wave vector ${\mathbf{k}}\in\mathbb{R}^{2}$ the sesquilinear form $a$ defined in (2.2) in Example 2.8 is Hermitian, continuous, and satisfies Gårding’s inequality for any $\beta>0$.

Proof.

Clearly, the sesquilinear form $a$ is Hermitian. For the continuity we apply the Cauchy-Schwarz inequality with respect to the complex dot product in $\mathbb{C}^{2}$ as well as the inner product in $\mathcal{H}$ and obtain for all $u,v\in\mathcal{V}$,

[TABLE]

Young’s inequality from Lemma 2.5 with $\delta=1$ then yields

[TABLE]

which proves the continuity of $a$. To show the Gårding inequality, we first consider the case ${\mathbf{k}}=0$. Then, for any $0<\beta\leq 1$, we have

[TABLE]

i.e., the Gårding inequality with $\alpha=\beta$. Otherwise, for ${\mathbf{k}}\neq 0$, we apply once more Lemma 2.5 for some parameter $\delta>0$ and get

[TABLE]

Now assume $\delta>1$ and define $\alpha\coloneqq(1-\delta^{-1})\min\{1,|{\mathbf{k}}|^{2}\}>0$. Then, $a$ satisfies the Gårding inequality with $\beta\coloneqq(\delta-\delta^{-1})\,|{\mathbf{k}}|^{2}>0$. Note that $\beta$ can be chosen arbitrarily small with an appropriate choice of $\delta>1$. ∎

Example 2.10** (TE mode).**

In the eigenvalue problem corresponding to the TE mode, the relative permittivity $\varepsilon$ appears in the differential operator. More precisely, for a fixed wave vector ${\mathbf{k}}$, we search for $u_{\mathbf{k}}$ such that

[TABLE]

for all $x\in\Omega$. Considering again periodic boundary conditions, we set $\mathcal{V}=H^{1}_{\text{per}}(\Omega)$ and $\mathcal{H}=L^{2}(\Omega)$ with the standard inner products. In this case, the sesquilinear form $a$ has the form

[TABLE]

The proof of the Gårding inequality follows similarly as in Lemma 2.9.

More general eigenvalue problems are discussed in Sections 4 and 5. There we consider a TM mode with an electric permittivity, which depends on the frequency $\omega$. This then leads to a nonlinear eigenvalue problem.

2.3. Operator formulation

The weak formulation of the eigenvalue problem (2.1) can be equivalently written as an operator equation in the (conjugate) dual space of $\mathcal{V}$, i.e., the space of conjugate linear and continuous mappings from $\mathcal{V}\to\mathbb{C}$. This then yields a convenient formulation for the introduction of iterative schemes. We introduce the operator $\mathcal{A}\colon\mathcal{V}\to\mathcal{V}^{*}$ by

[TABLE]

Further, we define $\mathcal{I}\colon\mathcal{V}\to\mathcal{V}^{*}$ as the embedding of $\mathcal{V}$ in $\mathcal{V}^{*}$ induced by the Gelfand triple $\mathcal{V}$, $\mathcal{H}$, $\mathcal{V}^{*}$ from Assumption 2.1 with respect to the inner product in $\mathcal{H}$, cf. [Zei90, Ch. 23.4]. To be precise, this means $\langle\mathcal{I}u,v\rangle=(u,v)$ for all $u,v\in\mathcal{V}$. The operator equation corresponding to (2.1) then reads

[TABLE]

Note that this equation is stated in the dual space of $\mathcal{V}$, which means that we consider test functions in $\mathcal{V}$ as in (2.1). Hence, the two formulations (2.1) and (2.3) are equivalent.

Finally, recall the definition of the shifted sesquilinear form $a_{\beta}$ from Remark 2.4. With this, we define the corresponding operator $\mathcal{A}_{\beta}:=\mathcal{A}+\beta\mathcal{I}\colon\mathcal{V}\to\mathcal{V}^{*}$, which is then positive and thus, invertible.

3. Iterative Methods on Operator Level

In this section, we analyze iterative methods to find the smallest eigenvalue as well as the corresponding eigenfunction of the operator eigenvalue problem (2.3). We emphasize that we do not apply any spatial discretization but perform the eigenvalue iteration directly to the operator equation.

We first consider the inverse power method for which we prove the convergence in $\mathcal{V}$ under the assumptions collected in the previous section. For this we consider two variants of the method and discuss the commutativity of spatial discretization and eigenvalue iteration, assuming an appropriate normalization within the scheme. Second, we discuss Arnoldi’s method for the operator case, based on Krylov subspaces as in the matrix case.

3.1. Inverse power method

Power and inverse power methods come in various variants with different kinds of scaling, see e.g. [AK08, Ch. 10.3] or [Saa11, Ch. 4] for the matrix case. One may even consider the scaling with the exact eigenvalue $\lambda$ as done, e.g., in [EE07, AHP18]. Clearly, the latter is only of interest for theoretical observations rather than actual computations.

As for every iteration scheme we assume a starting function $u^{0}\in\mathcal{V}$. In order to permit the iterates to converge to the wanted eigenfunction, one needs an additional assumption on $u^{0}$, e.g., having a non-vanishing component in the direction of this eigenfunction.

3.1.1. Rayleigh quotient iteration

A direct application of the inverse power method for the operator case would apply the inverse of the differential operator over and over again. However, this operator may not be invertible and the Rayleigh quotient $\langle\mathcal{A}u^{j},u^{j}\rangle/\langle\mathcal{I}u^{j},u^{j}\rangle$ is not guaranteed to remain positive. Thus, we consider the shifted eigenvalue problem

[TABLE]

This then leads to the following algorithm: Given an initial function $u^{0}\in\mathcal{V}$, $u^{0}\neq 0$, we solve for $j=1,2,\dots$ the variational problem

[TABLE]

Therein, $\tilde{u}^{j}:=u^{j}/\|u^{j}\|$ includes the normalization in $\mathcal{H}$, which is the natural norm corresponding to the eigenvalue problem (2.3). Further, $\mu^{j}$ denotes the Rayleigh quotient, i.e.,

[TABLE]

We emphasize that the iterates $u^{j}$ of the power iteration (3.2) are not normalized. The presence of the Rayleigh quotient within the iteration directly provides an approximation of the eigenvalue $\lambda$, given by $\lambda^{j}:=\mu^{j}-\beta$. It remains to discuss the convergence of the suggested iteration.

Theorem 3.1**.**

Given Assumptions 2.1-2.3, the power iteration (3.2) converges to an eigenpair (u^{\scalebox{0.5}{\bigstar}},\lambda^{\scalebox{0.5}{\bigstar}}) of (2.3) in the sense that the (sub)sequences $u^{j}$ and $\lambda^{j}:=\mu^{j}-\beta$ satisfy u^{j}\rightarrow u^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}$ and \lambda^{j}\rightarrow\lambda^{\scalebox{0.5}{\bigstar}} with \mathcal{A}u^{\scalebox{0.5}{\bigstar}}=\lambda^{\scalebox{0.5}{\bigstar}}\mathcal{I}u^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}^{*}$.

Proof.

We follow the proof in [Boz16] where the convergence of the inverse power method for the $p$-Laplacian is shown. The first step is to show the monotonic decrease of the sequence $\mu^{j}$. For this, we consider (3.2) with test functions $u^{j}$ and $\tilde{u}^{j-1}$ leading to

[TABLE]

A combination of these two estimates yields $\|u^{j}\|_{\beta}\leq\sqrt{\mu^{j-1}}\|u^{j}\|$ and thus,

[TABLE]

The monotonicity and $\mu^{j}>0$, which follows from the positivity of $\mathcal{A}_{\beta}$, imply that there exists the limit \mu^{\scalebox{0.5}{\bigstar}}:=\lim_{j\to\infty}\mu^{j}\geq 0. It even holds \mu^{\scalebox{0.5}{\bigstar}}>0, since the $\mu^{j}$ are uniformly bounded away from zero.

As a second step, we conclude from estimate (3.3) that $\|\tilde{u}^{j}\|_{\beta}=\sqrt{\mu^{j}}\leq\sqrt{\mu^{0}}$. Since the norms $\|\cdot\|_{\beta}$ and $\|\cdot\|_{\mathcal{V}}$ are equivalent, we obtain that $\tilde{u}^{j}$ is uniformly bounded in $\mathcal{V}$. Thus, there exists a convergent subsequence and an element {\tilde{u}}^{\scalebox{0.5}{\bigstar}}\in\mathcal{V}, which satisfy (without relabeling)

[TABLE]

Note that we have used here the compact embedding $\mathcal{V}\hookrightarrow\mathcal{H}$ from Assumption 2.2. Obviously, we have \|{\tilde{u}}^{\scalebox{0.5}{\bigstar}}\|=\lim_{j\to\infty}\|\tilde{u}^{j}\|=1. Further, we know from previous calculations that

[TABLE]

This means that the sequence $u^{j}$ is uniformly bounded in $\mathcal{V}$. Thus, there exists a limit u^{\scalebox{0.5}{\bigstar}}\in\mathcal{V} such that (again without relabeling of the subsequence)

[TABLE]

In the following we compare the two limits u^{\scalebox{0.5}{\bigstar}} and {\tilde{u}}^{\scalebox{0.5}{\bigstar}}. For this, we consider

[TABLE]

Taking the limit $j\to\infty$ on both sides, we conclude that

[TABLE]

i.e., \|u^{\scalebox{0.5}{\bigstar}}\|=1. As a result, the sequence $(\tilde{u}^{j})$ converges to u^{\scalebox{0.5}{\bigstar}}/\|u^{\scalebox{0.5}{\bigstar}}\|=u^{\scalebox{0.5}{\bigstar}}. The uniqueness of the limit then yields u^{\scalebox{0.5}{\bigstar}}={\tilde{u}}^{\scalebox{0.5}{\bigstar}}.

To show that the pair (u^{\scalebox{0.5}{\bigstar}},\lambda^{\scalebox{0.5}{\bigstar}}) with \lambda^{\scalebox{0.5}{\bigstar}}:=\mu^{\scalebox{0.5}{\bigstar}}-\beta is an eigenpair of (2.3), we apply the limit to equation (3.2). We conclude that (u^{\scalebox{0.5}{\bigstar}},\mu^{\scalebox{0.5}{\bigstar}}) indeed solves

[TABLE]

Thus, we have a(u^{\scalebox{0.5}{\bigstar}},v)=\lambda^{\scalebox{0.5}{\bigstar}}(u^{\scalebox{0.5}{\bigstar}},v) for all $v\in\mathcal{V}$ or, in operator form, \mathcal{A}u^{\scalebox{0.5}{\bigstar}}=\lambda^{\scalebox{0.5}{\bigstar}}\mathcal{I}u^{\scalebox{0.5}{\bigstar}}. The weak convergence \tilde{u}^{j}\rightharpoonup{\tilde{u}}^{\scalebox{0.5}{\bigstar}}=u^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}$ additionally implies, using (3.4) with $v=u^{\star}$,

[TABLE]

Note that the inequality is strict, if and only if the convergence is not strong. This implies \tilde{u}^{j},u^{j}\rightarrow u^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}$. If \mu^{\scalebox{0.5}{\bigstar}} is a single eigenvalue, then every convergent subsequence has the same limit. ∎

Remark 3.2*.*

In order to show that (u^{\scalebox{0.5}{\bigstar}},\lambda^{\scalebox{0.5}{\bigstar}}) is the smallest eigenpair, one needs an additional assumption on $u^{0}$. In the real case, one may consider $u^{0}\geq 0$, using the maximum principle, see e.g. [Boz16].

Remark 3.3*.*

A second strategy to prove Theorem 3.1 is to reformulate the eigenvalue problem in terms of the resolvent. Assumption 2.2 shows that the resolvent is a compact operator such that the results in [ESL95] can be applied.

3.1.2. An alternative power method

In Theorem 3.1 we have shown that the inverse power method (3.2) provides in the limit an eigenvalue (the limit of the Rayleigh quotient) and an eigenfunction. However, one may also omit the Rayleigh quotient, which leads to the iteration

[TABLE]

and the following convergence result.

Lemma 3.4**.**

Assume $u^{0}=v^{0}\in\mathcal{V}$ with $\|u^{0}\|=1$. Let $u^{j}$ and $v^{j}$ be the sequences obtained from the iteration procedures (3.2) and (3.5), respectively. Then, we have the relation $u^{j}=\mu^{j-1}v^{j}$ for all $j\geq 1$ with the Rayleigh quotient $\mu^{j}=\|u^{j}\|_{\beta}^{2}/\|u^{j}\|^{2}$.

Proof.

We prove this result by mathematical induction and observe first that

[TABLE]

Now, assuming that $u^{j}=\mu^{j-1}v^{j}$ is true for a fixed but arbitrary index $j$, we obtain

[TABLE]

Note that we have used the fact that $\mu^{j}>0$ and thus $|\mu^{j}|=\mu^{j}$. ∎

Lemma 3.4 directly implies the convergence of the iteration (3.5). In particular, we have v^{j}=u^{j}/\mu^{j-1}\to u^{\scalebox{0.5}{\bigstar}}/\mu^{\scalebox{0.5}{\bigstar}}=:v^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}$. Thus, the limit is not normalized but rather satisfies

[TABLE]

3.1.3. Commutativity of discretization and power iteration

The convergence of the power method in the operator case directly leads to the question whether the application of the power method and the spatial discretization commute. If we discretize the shifted eigenvalue problem (3.1) by finite elements, then we obtain a system of the form

[TABLE]

Therein, $q\in\mathbb{C}^{n}$ encodes a finite-dimensional approximation of the eigenfunction $u\in\mathcal{V}$, e.g., the coefficients w.r.t. a finite element basis. Since we have included the boundary conditions in the space $\mathcal{V}$, we assume that the stiffness matrix $K\in\mathbb{C}^{n,n}$ and the mass matrix $M\in\mathbb{C}^{n,n}$ are Hermitian and positive definite. With $A:=M^{-1}K$ the discrete system is equivalent to the eigenvalue problem $Aq=\mu q$. Seeking for the smallest eigenvalue, we apply the inverse power method with normalization, i.e.,

[TABLE]

where $Z$ is a symmetric and positive definite matrix and $\|q\|_{Z}:=(q^{T}Zq)^{1/2}$ the corresponding norm. As approximation for the eigenvalue we consider $\mu^{j}:=1/\|q^{j}\|_{Z}$. The starting vector is denoted by $q^{0}$. We emphasize that the iteration converges despite of the choice of the norm as long as $q^{0}$ contains a non-zero component in direction of the first eigenvector.

We now consider the spatial discretization of (3.5), i.e., we first apply the inverse power method to the PDE eigenvalue problem (3.1) and then discretize. We emphasize that this allows a different discretization scheme in each iteration step and thus, adaptivity. If we consider, however, the same spatial mesh as before in each iteration step, then the same matrices $K$ and $M$ appear and we get

[TABLE]

Note that the applied normalization is here w.r.t the $M$-norm, since this corresponds to the $L^{2}$-inner product in the infinite-dimensional case. Thus, the resulting iteration equals (3.6) if we choose the normalization matrix $Z=M$.

3.2. A Krylov subspace method

The natural extension of the inverse power method in order to approximate several eigenvalues is a subspace iteration, cf. [Saa11, Ch. 5]. This includes several starting functions, for which a power iteration is applied, and an additional orthogonalization step is performed. The computation of several eigenvalues is also of interest in the calculation of band-gaps of a photonic crystal. Here, we consider the Arnoldi method in the operator setting. For this, we need an extension of the Krylov subspaces used in numerical linear algebra.

3.2.1. Krylov spaces

Krylov subspaces play a crucial role for iterative eigenvalue computations, see, e.g., [Saa11, Ch. 6.1]. In order to generalize these methods to the PDE setting, we need Krylov subspaces for general Hilbert spaces [GHS14].

Let $u^{0}$ be a function in $\mathcal{V}$, e.g., an initial guess for the power method in Section 3.1. With this, we define the Krylov subspace

[TABLE]

Obviously, this defines a closed subspace of $\mathcal{V}$. We emphasize that $\mathfrak{K}_{\beta}^{m}(u^{0})$ is spanned – as in the finite-dimensional case – by the iterates of the power method. To see this, note that $(\mathcal{A}_{\beta}^{-1}\mathcal{I})^{j}u^{0}$ equals the corresponding iterate of the power method up to a multiplicative constant. Thus, if we denote the sequence resulting from the modified inverse power method (3.5) by $u^{j}_{\text{pow}}$, then we have

[TABLE]

3.2.2. Arnoldi’s method

We translate the Arnoldi algorithm from the matrix setting [Saa11, Ch. 6.2] to the present operator case. A similar infinite-dimensional Arnoldi algorithm is considered in [JMM12]. Let $\{v_{1},\dots,v_{m}\}$ denote a basis of $\mathfrak{K}_{\beta}^{m}(u^{0})$, e.g., obtained by a Gram-Schmidt orthogonalization process. Then, the new iterate of the Arnoldi method is given by $u^{m}:=\sum_{j=1}^{m}\zeta_{j}v_{j}\in\mathfrak{K}_{\beta}^{m}(u^{0})$, whose coefficients $\zeta:=[\zeta_{1},\dots,\zeta_{m}]^{T}$ and corresponding $\mu^{m}\in\mathbb{R}$ are derived by the Galerkin projection

[TABLE]

This is equivalent to the $m\times m$ eigenvalue problem $\tilde{K}\zeta=\mu^{m}\tilde{M}\zeta$, for which we search for the smallest eigenpair. Here, $\tilde{K}$ and $\tilde{M}$ are stiffness and mass matrices restricted to the Krylov basis, i.e.,

[TABLE]

Thus, the extra costs going from the power method to the Arnoldi method are identical as in the finite-dimensional setting, namely the solution of a small (but dense) eigenvalue problem. Note that the resulting approximation of the eigenvalue, namely $\mu^{m}$, is again the Rayleigh quotient of the iterate $u^{m}$.

Remark 3.5*.*

The commutativity result for the inverse power method does not carry over to the Arnoldi method, since the discretization of the operator $(\mathcal{A}_{\beta}^{-1}\mathcal{I})^{j}$ leads, in general, not to the matrix $(K^{-1}M)^{j}$.

The obtained pair of the Arnoldi method provides the best-approximation within the Krylov subspace $\mathfrak{K}_{\beta}^{m}(u^{0})$ in the sense that

[TABLE]

for all $v\in\mathfrak{K}_{\beta}^{m}(u^{0})$ and with the residual defined by $\operatorname{Res}(u,\mu)\coloneqq\mathcal{A}_{\beta}u-\mu\mathcal{I}u\in\mathcal{V}^{*}$ for any pair $(u,\mu)\in\mathcal{V}\times\mathbb{C}$. Obviously, this implies that the Arnoldi method is superior to the inverse power method. The norm of the residual may also be used as an error estimator as it can be transformed to the backward error. Thus, small residuals indicate good approximations of the eigenpair [Mie11, Ch. 4].

The gain of the Arnoldi method can also be characterized in terms of the Courant min-max principle in Hilbert spaces, cf. [Cou20] or [WS72, Ch. 1]. This means that the $\ell$th eigenvalue is defined by minimizing over all $\ell$-dimensional subspaces of $\mathcal{V}$, i.e.,

[TABLE]

With this, one shows that computed approximations of the eigenvalues are larger than the exact ones. With the same arguments one can show that the Arnoldi iteration provides better approximations than the inverse power method and thus, converges as well. More precisely, the Arnoldi method computes an approximation $\mu^{m}_{\text{Arnoldi}}$ satisfying

[TABLE]

Note that the inequality holds, since $u^{m-1}_{\text{pow}}$ is an element of the Krylov subspace $\mathfrak{K}_{\beta}^{m}(u^{0})$ and thus, $\operatorname{span}\{u^{m-1}_{\text{pow}}\}$ is a particular one-dimensional subspace.

4. The Inverse Power Method for a Nonlinear Model Problem

In this section, we consider an extension of the TM mode from Example 2.8, in which the relative electric permittivity $\varepsilon$ depends on the frequency and thus, the eigenvalue. This then leads to a nonlinear eigenvalue problem. Assuming that $\varepsilon$ is a rational function in the frequency, we are able to reformulate the eigenvalue problem to a linear one satisfying Assumptions 2.1 and 2.3. This linearization, however, leads to a lack of compactness, which calls for a novel convergence analysis of the inverse power method. This then leads to a slightly weaker convergence result, compared to Section 3.

4.1. A simplified Drude-Lorentz model

We consider a photonic crystal made up of two different materials. For this, we decompose the computational domain $\Omega$ into two subdomains $\Omega=\Omega_{1}\cup\Omega_{2}$, each representing one material. We define the corresponding indicator functions on $\Omega_{j}$ by $\chi_{j}\colon\Omega\to\{1,0\}$, $j=1,2$. On both subdomains we assume the relative electric permittivity to be constant in space and thus,

[TABLE]

For the sake of brevity, the material contained in $\Omega_{1}$ is assumed to be linear, i.e., we set the relative permittivity in this subdomain to a constant $\varepsilon_{1}(\omega)\equiv\alpha_{1}>0$. For the frequency dependence in the second material we consider a simplified version of the Drude-Lorentz model, see e.g. [LL10] or [Jac99, Ch. 7.5]. More precisely, we assume the electric permittivity to be of the form

[TABLE]

with a positive constant $\alpha_{2}>0$ and real parameters $\eta_{\ell},\xi_{\ell}$ such that $\eta_{\ell}^{2},\xi_{\ell}^{2}\geq 0$. This corresponds to the ’lossless’ case considered in [Eng10]. In order to stay bounded, we only consider $\omega$ away from the poles given by $\eta_{\ell}$. Since $\omega$ appears in (4.2) only squared, we introduce $\lambda\coloneqq\omega^{2}$ and obtain

[TABLE]

with

[TABLE]

Due to the inclusion of $\Xi$, the fractional terms are in a strictly proper form, i.e., the degree of the polynomial in terms of $\lambda$ in the numerator is strictly smaller than the degree in the denominator. All in all, this leads to the nonlinear eigenvalue problem

[TABLE]

with $\nabla_{\mathbf{k}}$ denoting again the shifted gradient introduced in Example 2.8. Our aim is to turn this into a linear eigenvalue problem in order to apply the iterative methods from the previous section. For this, we follow the ideas presented in [SB11, EKE12, Eff13], which consider the corresponding finite-dimensional case. The main clue is to rewrite the sum, which may be regarded as a transfer function, by means of a realization, i.e.,

[TABLE]

Due to the simple structure of the permittivity, we can directly read off the vector $\mathbf{b}=[\xi_{1}\eta_{1},\ldots,\xi_{L}\eta_{L}]^{T}\in\mathbb{\mathbb{R}}^{L}$ and $\mathbf{A}\in\mathbb{R}^{L\times L}$ as the diagonal (and thus Hermitian) matrix with $\mathbf{A}_{j,j}=\eta_{j}^{2}$ for $j=1,\ldots,L$. By $\mathbf{I}\in\mathbb{R}^{L\times L}$ we denote the identity matrix.

Remark 4.1*.*

For other models of the permittivity, the choice (and even the dimension) of $\mathbf{A}$ and $\mathbf{b}$ may not be as straightforward. Further, the identity matrix $\mathbf{I}$ may be replaced by another positive Hermitian matrix. Proper realizations for such cases may be found using the techniques discussed in [SUBG18].

4.2. Spaces and embeddings

For the weak formulation and linearization of the eigenvalue problem (4.1) we introduce several function spaces. First, we introduce

[TABLE]

These spaces form Hilbert spaces equipped with the inner products

[TABLE]

for $u$ and $v$ in the respective spaces $H$, $V$, or $X$. Second, we define the product spaces

[TABLE]

with $X^{L}:=X\otimes\dots\otimes X$. Also these product spaces are Hilbert spaces. In $\mathcal{H}$ we consider the inner product

[TABLE]

for $\bm{z}_{1}=[u;\bm{x}]$, $\bm{z}_{2}=[v;\bm{y}]\in\mathcal{H}$ consisting of $u,v\in H$ and $\bm{x},\bm{y}\in X^{L}$ with $\bm{x}=[x_{1};\dots;x_{L}]$, $\bm{y}=[y_{1};\dots;y_{L}]$. Note that we use here the notation $[u;\bm{x}]:=[u,\bm{x}^{T}]^{T}$. Analogously, we define an inner product in $\mathcal{V}$ by replacing $(u,v)$ by $(u,v)_{V}$. We also define the corresponding norms $\|\bm{z}\|_{\mathcal{H}}^{2}:=(\bm{z},\bm{z})_{\mathcal{H}}$ and $\|\bm{z}\|_{\mathcal{V}}^{2}:=(\bm{z},\bm{z})_{\mathcal{V}}$.

Remark 4.2*.*

Although the embedding $V\hookrightarrow H$ is compact, the embedding $\mathcal{V}\hookrightarrow\mathcal{H}$ is not. This is due to the fact that the identity operator is not compact in infinite dimensions.

For the weak formulation we further need several embeddings. First, $\mathcal{I}\colon V\to V^{*}$ denotes the continuous inclusion map defined by the Gelfand triple $V$, $H$, $V^{\ast}$, cf. Section 2.3. Second, we define the extension of the mapping $u\mapsto(u,\cdot\,)_{X}$ as $\mathcal{I}_{2}\colon V\to V^{*}$, i.e.,

[TABLE]

In the same manner, but based on the indicator function $\chi_{1}$, we define $\mathcal{I}_{1}\colon V\to V^{*}$. The weighted combination of these two embeddings yields $\mathcal{I}_{\alpha}\colon V\to V^{*}$, given by

[TABLE]

Finally, we introduce the embedding $\overline{\mathcal{I}}_{2}\colon X\to V^{\ast}$. For $u\in X$ this is defined by

[TABLE]

The corresponding dual operator $\overline{\mathcal{I}}_{2}^{*}\colon V\to X^{*}$ satisfies for $v\in V$ and $u\in X$ that

[TABLE]

Lemma 4.3**.**

With the Riesz isomorphism $j_{X}\colon X\to X^{\ast}$, the introduced embeddings satisfy the identity $\overline{\mathcal{I}}_{2}\,j_{X}^{-1}\,\overline{\mathcal{I}}_{2}^{\ast}=\mathcal{I}_{2}\colon V\to V^{*}$.

Proof.

Consider $u,v\in V$. Then, the claimed identity can be seen by

[TABLE]

Remark 4.4*.*

In the remainder of this paper, we often omit to write the Riesz isomorphisms $j_{H}\colon H\to H^{\ast}$ or $j_{X}\colon X\to X^{\ast}$ if their presence is clear from the context. Thus, we may write $\overline{\mathcal{I}}_{2}\overline{\mathcal{I}}_{2}^{\ast}=\mathcal{I}_{2}$.

4.3. Weak formulation and linearization

This subsection is devoted to the transformation of the nonlinear eigenvalue problem into a linear one by the introduction of new variables. Thus, we aim to write (4.1) in the form

[TABLE]

Based on the proper form of the permittivity given in (4.4), we obtain the following weak form of the eigenvalue problem. For a given (and fixed) wave vector ${\mathbf{k}}\in\mathcal{K}$, find a non-trivial pair $(u_{\mathbf{k}},\lambda)\in V\times\mathbb{R}$ such that

[TABLE]

We emphasize that this equation is stated in $V^{*}$. The operator $\mathcal{A}_{{\mathbf{k}}}\colon V\to V^{*}$ denotes the weak form of the shifted Laplacian, cf. Example 2.8. For the first order formulation of (4.5) we introduce a new variable

[TABLE]

Note that this includes a hidden application of $j_{X}^{-1}$. With Lemma 4.3 this leads to a linear eigenvalue problem where we search for a pair $(\bm{z},\lambda)$ with $\bm{z}:=[u_{\mathbf{k}};\bm{x}]\in\mathcal{V}$ such that

[TABLE]

This formulation consists of two equations stated in the dual spaces of $V$ and $X^{L}$, respectively. The Riesz isomorphism $j_{X}$ should be understood as a componentwise application.

Remark 4.5*.*

The operator on the left-hand side of the first order eigenvalue problem (4.7) has a generalized saddle point structure.

Lemma 4.6**.**

The eigenvalue problems (4.5) and (4.7) are equivalent. This means that an eigenpair $(u_{\mathbf{k}},\lambda)$ of (4.5) defines a solution of (4.7) by $([u_{\mathbf{k}};\bm{x}],\lambda)$ with $\bm{x}$ defined as in (4.6) and vice versa.

Proof.

The second block row of (4.7) is given by $(\mathbf{A}-\lambda\mathbf{I})j_{X}\bm{x}=\mathbf{b}\,\overline{\mathcal{I}}_{2}^{\ast}u_{\mathbf{k}}$, which implies (4.6). Substituting this formula for $\bm{x}$ into the first block row yields together with $\overline{\mathcal{I}}_{2}\overline{\mathcal{I}}_{2}^{\ast}=\mathcal{I}_{2}$ from Lemma 4.3 that (4.7) is indeed equivalent to the nonlinear eigenvalue problem (4.5). ∎

Defining $\mathbb{A},\mathbb{I}\colon\mathcal{V}\to\mathcal{V}^{*}$ in an obvious manner, we can write (4.7) in the form $\mathbb{A}\bm{z}=\lambda\,\mathbb{I}\bm{z}$. The sesquilinear form corresponding to $\mathbb{I}$ reads $\mathfrak{i}\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$,

[TABLE]

for $\bm{z}_{1}=[u;\bm{x}_{1}]$, $\bm{z}_{2}=[v;\bm{x}_{2}]\in\mathcal{V}$. Note that we may also consider $\mathfrak{i}$ as a sesquilinear form mapping from $\mathcal{H}\times\mathcal{H}$ to $\mathbb{C}$. We now show that $\mathfrak{i}$ defines an inner product in $\mathcal{H}$.

Lemma 4.7**.**

The sesquilinear form $\mathfrak{i}\colon\mathcal{H}\times\mathcal{H}\to\mathbb{C}$, defines an inner product in $\mathcal{H}$.

Proof.

Obviously, $\mathfrak{i}$ is Hermitian and sesquilinear. Further, for any $\bm{z}=[u;\bm{x}]\in\mathcal{H}$ and $\underline{\alpha}:=\min\{\alpha_{1},\alpha_{2}\}>0$, it holds that

[TABLE]

which proves its positivity. ∎

4.4. Shifted eigenvalue problem

In the previous subsection, the nonlinear eigenvalue problem was brought into the form $\mathbb{A}\bm{z}=\lambda\,\mathbb{I}\bm{z}$. In order to use the framework presented in Section 2, we need to apply a shift to gain positivity of the differential operator. Recall from the discussion of the linear case that the operator $\mathcal{A}_{{\mathbf{k}}}$ is not elliptic. From Lemma 2.9 we know, however, that $\mathcal{A}_{{\mathbf{k}},\tilde{\beta}}:=\mathcal{A}_{{\mathbf{k}}}+\tilde{\beta}\mathcal{I}$ is elliptic for every $\tilde{\beta}>0$.

For fixed $\tilde{\beta}>0$ and $\underline{\alpha}:=\min\{\alpha_{1},\alpha_{2}\}$ we introduce $\beta:=\tilde{\beta}/\underline{\alpha}$ and shift the linearized eigenvalue problem by $\beta\mathbb{I}\bm{z}$. This then provides the appearance of $\mathcal{A}_{{\mathbf{k}},\tilde{\beta}}$. More precisely, we obtain the shifted operator

[TABLE]

Due to $\mathcal{I}=\mathcal{I}_{1}+\mathcal{I}_{2}$ we have

[TABLE]

with \big{(}{\underline{\alpha}}^{-1}\mathcal{I}_{\alpha}-\mathcal{I}\big{)}\geq 0. The shifted eigenvalue problem has the form $\mathbb{A}_{\beta}\bm{z}=(\lambda+\beta)\,\mathbb{I}\,\bm{z}$. Corresponding to the operator $\mathbb{A}_{\beta}$, we define $\mathfrak{a}_{\beta}\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$ by

[TABLE]

for $\bm{z}_{1}=[u;\bm{x}_{1}]$, $\bm{z}_{2}=[v;\bm{x}_{2}]\in\mathcal{V}$.

Lemma 4.8**.**

The sesquilinear form $\mathfrak{a}\colon\mathcal{V}\times\mathcal{V}\to\mathbb{C}$ defined through $\mathfrak{a}(\cdot\,,\cdot):=\langle\mathbb{A}\,\cdot\,,\cdot\rangle$ satisfies Assumption 2.3.

Proof.

Due to the given structure of $\mathbb{A}$ and the fact that $\mathbf{A}$ is Hermitian, it is easy to see that $\mathfrak{a}$ is continuous and Hermitian. It remains to show that $\mathfrak{a}_{\beta}$ is positive for some shift $\beta>0$. For this, we consider $\bm{z}=[u;\bm{x}]\in\mathcal{V}$ and note that

[TABLE]

The definition of $\mathbf{b}=[\xi_{1}\eta_{1},\ldots,\xi_{L}\eta_{L}]^{T}$ yields the estimate

[TABLE]

and thus, with $c_{\text{ell}}$ denoting the ellipticity constant of $\mathcal{A}_{{\mathbf{k}},\tilde{\beta}}$,

[TABLE]

Assuming $\beta>\max\eta_{\ell}^{2}-\min\eta_{\ell}^{2}$, we conclude the positivity of $\mathfrak{a}_{\beta}$. ∎

4.5. Convergence of the inverse power method

We apply the inverse power method from Section 3.1 to the shifted eigenvalue problem $\mathbb{A}_{\beta}\bm{z}=(\lambda+\beta)\,\mathbb{I}\,\bm{z}$. For an initial function $\bm{z}^{0}=[u^{0};\bm{x}^{0}]\in\mathcal{V}$ we thus consider the iteration

[TABLE]

with the Rayleigh quotient

[TABLE]

Note that we use here the norms $\|\bm{z}\|_{\beta}^{2}=\mathfrak{a}_{\beta}(\bm{z},\bm{z})$ and $\|\bm{z}\|^{2}=\mathfrak{i}(\bm{z},\bm{z})$ in line with Section 2.1. Further, $\tilde{}\bm{z}^{j-1}$ denotes the normalization in the $\mathfrak{i}$-norm, i.e., we define $\tilde{}\bm{z}^{j-1}:=\bm{z}^{j-1}/\|\bm{z}^{j-1}\|$. Note, however, that in the present setting the embedding $\mathcal{V}\hookrightarrow\mathcal{H}$ is not compact and thus, Assumption 2.2 is not satisfied. Hence, Theorem 3.1 is not applicable but we are able to show the following result.

Theorem 4.9**.**

Consider the nonlinear eigenvalue problem (4.5) for a fixed wave vector ${\mathbf{k}}$ and a starting function $u^{0}\in V$. Set $\bm{z}^{0}:=[u^{0};\bm{x}^{0}]\in\mathcal{V}$ with any $\bm{x}^{0}\in X^{L}$. Then, the power iteration (4.9) converges in the sense that there exists a subsequence of $u^{j}$, which satisfies u^{j}\rightarrow u^{\scalebox{0.5}{\bigstar}} in $H$ with u^{\scalebox{0.5}{\bigstar}} being an eigenfunction of (4.5).

Proof.

We proceed as in the proof of Theorem 3.1 and test (4.9) with $\bm{z}^{j}$ and $\tilde{}\bm{z}^{j-1}$, respectively. This then yields the estimates

[TABLE]

Due to $\mu^{j}\geq 0$, we conclude the existence of a limit \mu^{\scalebox{0.5}{\bigstar}}:=\lim_{j\to\infty}\mu^{j}\geq 0. This also implies the uniform bounds

[TABLE]

For the last estimate we have used the continuity of the embedding $\mathcal{V}\hookrightarrow\mathcal{H}$ and the ellipticity of $\mathfrak{a}_{\beta}$ shown in Lemma 4.8. Thus, there exist convergent subsequences and limits {\bm{z}}^{\scalebox{0.5}{\bigstar}},{\tilde{}\bm{z}}^{\scalebox{0.5}{\bigstar}}\in\mathcal{V}, which satisfy (without relabeling) \bm{z}^{j}\rightharpoonup{\bm{z}}^{\scalebox{0.5}{\bigstar}},\tilde{}\bm{z}^{j}\rightharpoonup{\tilde{}\bm{z}}^{\scalebox{0.5}{\bigstar}} in $\mathcal{V}$. We emphasize that the two limits can only differ by a multiplicative constant, i.e., {\tilde{}\bm{z}}^{\scalebox{0.5}{\bigstar}}=c\,{\bm{z}}^{\scalebox{0.5}{\bigstar}}. A componentwise consideration with $\bm{z}^{j}=[u^{j};\bm{x}^{j}]$, $\tilde{}\bm{z}^{j}=[\tilde{u}^{j};\tilde{\bm{x}}^{j}]$, and {\bm{z}}^{\scalebox{0.5}{\bigstar}}=[u^{\scalebox{0.5}{\bigstar}};{\bm{x}}^{\scalebox{0.5}{\bigstar}}] then yields

[TABLE]

Using the compact embedding $V\hookrightarrow H$, we conclude that the first component converges strongly in $H$, i.e., u^{j}\rightarrow u^{\scalebox{0.5}{\bigstar}} and \tilde{u}^{j}\rightarrow c\,u^{\scalebox{0.5}{\bigstar}} in $H$. We now show that the limit pair (u^{\scalebox{0.5}{\bigstar}},\lambda^{\scalebox{0.5}{\bigstar}}) with \lambda^{\scalebox{0.5}{\bigstar}}:=c\,\mu^{\scalebox{0.5}{\bigstar}}-\beta solves the nonlinear eigenvalue problem (4.5). For this, we apply to (4.9) a test function $[v;0]\in\mathcal{V}$ with arbitrary $v\in V$ and consider the limit $j\to\infty$,

[TABLE]

By the definitions of $\mathfrak{a}_{\beta}$ and $\mathfrak{i}$ we conclude that

[TABLE]

On the other hand, taking the limit in (4.9) with a test function $[0;\bm{y}]\in\mathcal{V}$, we obtain

[TABLE]

These two equations together then give

[TABLE]

Thus, the pair (u^{\scalebox{0.5}{\bigstar}},\lambda^{\scalebox{0.5}{\bigstar}}) is an eigenpair of the nonlinear eigenvalue problem (4.5). ∎

This result shows that the convergence in the eigenfunction is maintained for this particular case coming from the linearization of a nonlinear eigenvalue problem. In contrast to Section 3, however, we only showed the convergence in $H$ and not in $V$ due to the lack of compactness.

5. A Newton Method for Nonlinear Eigenvalue Problems Arising in Photonic Crystals

Similarly as in Section 4, we consider an extension of Example 2.8 where the electric permittivity $\varepsilon$ is frequency-dependent and has the form (4.1). Here, however, we allow for a more general class of models for the electric permittivity $\varepsilon_{2}(\omega)$. Indeed, the simplified Drude-Lorentz model considered in Section 4 has the particular property that it can be linearized in a way that ensures Assumption 2.3. But other – more realistic – models may not have this nice property, if complex poles exist in $\varepsilon_{2}(\omega)$. Indeed, this rules out the possibility to obtain a realization analogous to (4.4) where $\mathbf{A}$ and $\mathbf{I}$ are Hermitian matrices and $\mathbf{I}$ is positive definite. Because of this, we follow here an alternative strategy and tackle the nonlinear problem directly by a Newton method.

We are especially interested in the general causality-preserving model described in [GVDZ17] and the full Drude-Lorentz model, see e.g. [LL10] or [Jac99, Ch. 7.5], i.e.,

[TABLE]

for some empirically determined parameters $A_{\ell},B_{\ell}$ and, respectively, $\alpha_{\text{DL}}$, $\xi_{\ell}$, $\eta_{\ell}$, $\gamma_{\ell}$ for $\ell=1,\dots,L$. Here we restrict ourselves to the real part of these models so that the resulting eigenvalue problem (5.3) is self-adjoint, i.e., we consider functions of the form

[TABLE]

and

[TABLE]

The information about the imaginary part can be included later through perturbation theory as a post-processing step, see [RF11]. After discussing the properties of the resulting nonlinear eigenvalue problem in details, we show how to solve it with a Newton iteration, provided some a priori knowledge on the eigenpair of interest is given. First, however, we discuss the application of Newton’s method in the linear case.

5.1. Newton’s method for linear eigenvalue problems

Even for the linear case $\mathcal{A}u=\lambda\,\mathcal{I}u$ one may apply Newton’s method, cf. [MV04, Sch08]. For this, we rewrite the eigenvalue problem as

[TABLE]

with a fixed function $y\in\mathcal{V}$, $(y,u)\neq 0$, serving as normalization. Further, the second equation balances the number of equations and variables. The resulting Newton iteration reads

[TABLE]

This means that all iterates are normalized to $(y,u^{j})=1$ and $u^{j}$ satisfies

[TABLE]

The similarities to the power method are obvious. More precisely, Newton’s method leads to a shifted inverse iteration with the shift given by the previous eigenvalue approximation. Considering the constant shift zero, we obtain – up to scaling – the power iteration given in (3.2). The precise algorithm is given in [MV04, Alg. I]. In the matrix case, this then leads to a local third-order convergence of the eigenvalue, cf. [Osb64].

5.2. Definition of the nonlinear eigenvalue problem

In the remainder of this section, we assume that $\varepsilon_{1}(\omega)\equiv\alpha_{1}\in\mathbb{R}^{+}$ and that $\varepsilon_{2}$ is a real-valued function which is analytic on an open connected set $S\subset\mathbb{R}$ where it satisfies

[TABLE]

for all $\omega\in S$.

Given $\mathcal{V}=H^{1}_{\text{per}}(\Omega)$, $\mathcal{H}=L^{2}(\Omega)$ with the standard $H^{1}$ and $L^{2}$-norm, respectively, we consider the nonlinear eigenvalue problem

[TABLE]

where $\mathcal{A}\colon\mathcal{V}\to\mathcal{V}^{*}$ is defined by the sesquilinear form (2.2) via $\langle\mathcal{A}u,v\rangle_{\mathcal{V}^{*},\mathcal{V}}=a(u,v)$ for all $u,v\in\mathcal{V}$ and $\mathcal{B}(\omega)\colon S\to\mathcal{L}(\mathcal{V},\mathcal{V}^{*})$ is defined by

[TABLE]

Here, $\mathcal{L}(X,Y)$ denotes the set of bounded linear operators mapping from $X$ to $Y$.

An eigenvalue problem similar to (5.3), but formulated in terms of operators mapping from $\mathcal{V}$ into itself, was already analyzed in [Eng10]. Some properties of the operator-valued function $\mathcal{T}(\omega)$ can be directly derived from there.

Lemma 5.1**.**

The operator-valued function $\mathcal{T}(\omega)$ can be extended to a self-adjoint holomorphic operator-valued function in some neighborhood $D{\subset\mathbb{C}}$ of $S$, where $D$ is symmetric with respect to the real axis. Moreover, the spectrum of $\mathcal{T}(\omega)$ consists of isolated eigenvalues of finite multiplicity.

Proof.

See the proof of Lemma 4.4 in [Eng10] and the following remark. ∎

Lemma 5.2**.**

For all $\omega\in D$ the operator $\mathcal{T}(\omega)\colon\mathcal{V}\to\mathcal{V}^{*}$ is Fredholm with index 0.

Proof.

In [Eng10] it is shown that the operator $T(\omega)\colon\mathcal{V}\to\mathcal{V}$ defined by $(T(\omega)u,v)_{\mathcal{V}}:=\langle\mathcal{T}(\omega)u,v\rangle_{\mathcal{V}^{*},\mathcal{V}}$ for all $u,v\in\mathcal{V}$ is Fredholm with index 0. Moreover, the Riesz isomorphism $j_{\mathcal{V}}$ between $\mathcal{V}$ and $\mathcal{V}^{*}$ is a bounded and bijective linear operator, thus it is Fredholm. Further, its index is 0 since $\dim(\ker(j_{\mathcal{V}}))=\text{codim}(\text{range}(j_{\mathcal{V}}))=0$. It follows that $\mathcal{T}(\omega)=j_{\mathcal{V}}\,T(\omega)$ is a composition of two Fredholm operators with index [math] and thus, $\mathcal{T}(\omega)$ is Fredholm with index [math] itself. ∎

We denote by (u^{\scalebox{0.5}{\bigstar}},\omega^{\scalebox{0.5}{\bigstar}})\in\mathcal{V}\times S an eigenpair of (5.3) and assume the multiplicity of the eigenvalue \omega^{\scalebox{0.5}{\bigstar}} to be $1$.

5.3. A Newton iteration for the nonlinear eigenvalue problems

In this subsection, we define a Newton iteration. We will prove convergence to the eigenpair (u^{\scalebox{0.5}{\bigstar}},\omega^{\scalebox{0.5}{\bigstar}}) under certain conditions discussed in the next subsection. For this, we extend the results of [Sch08, Ch. 4] to the infinite-dimensional case and therefore define the iteration following a similar reasoning. We start with rewriting the eigenvalue problem (5.3) as

[TABLE]

where the functional $\mathcal{P}_{[y]}\colon\mathcal{V}\to\mathbb{C}$ is defined by $\mathcal{P}_{[y]}u\coloneqq(y,u)$ for all $u\in\mathcal{V}$ and the normalizing vector $y\in\mathcal{V}$ has to be chosen such that \mathcal{P}_{[y]}u^{\scalebox{0.5}{\bigstar}}\neq 0. Note that the self-adjointness of $\mathcal{T}(\omega)$ for all $\omega\in S$ implies that the eigenvalues of $\mathcal{T}$ are real. However, since $\mathcal{V}$ is defined over the field $\mathbb{C}$, the well-posedness of the Newton iteration requires an extension of the domain of definition of $\mathcal{T}$ to an open and convex neighborhood $D\subset\mathbb{C}$ of $S$, cf. Lemma 5.1.

As $\mathcal{F}_{[y]}$ depends linearly on the first argument $u$ and depends on the second argument $\omega$ only through the holomorphic function $\mathcal{T}$, it follows that $\mathcal{F}_{[y]}$ is thrice continuously Fréchet differentiable on $E\coloneqq G\times D$ where $G$ is an open and convex subspace of $\mathcal{V}$ containing u^{\scalebox{0.5}{\bigstar}}. A Taylor expansion of $\mathcal{F}_{[y]}$ around $(u,\omega)\in E$ yields

[TABLE]

where $\|\left(u,\ \omega\right)\|_{\mathcal{W}}^{2}\coloneqq\|u\|_{\mathcal{V}}^{2}+|\omega|^{2}$ for all $(u,\omega)\in\mathcal{W}\coloneqq\mathcal{V}\times\mathbb{C}$. The corresponding Jacobian is given by

[TABLE]

where $\mathcal{T}^{\prime}$ denotes the derivative of $\mathcal{T}$ with respect to $\omega$. In order to shorten the notation, we introduce the variables \Delta u\coloneqq u^{\scalebox{0.5}{\bigstar}}-u and \Delta\omega\coloneqq\omega^{\scalebox{0.5}{\bigstar}}-\omega. The second Fréchet derivative of $\mathcal{F}_{[y]}$ along the direction $(\Delta u,\,\Delta\omega)$ reads

[TABLE]

Following the strategy of the standard formulation of Newton’s method, we define the following iterative process: Given an initial vector $(u^{0},\,\omega^{0})\in E$ satisfying $\mathcal{P}_{[y]}u^{0}=1$, the iterates are defined by $u^{j+1}:=u^{j}+s^{j}$, $\omega^{j+1}:=\omega^{j}+\nu^{j}$ with

[TABLE]

for $j=0,1,2,\ldots$.

We show in the following subsection that (5.5) is uniquely solvable for $(s^{j},\,\nu^{j})$ when $(u^{j},\,\omega^{j})$ lies in a neighbourhood of the eigenpair (u^{\scalebox{0.5}{\bigstar}},\,\omega^{\scalebox{0.5}{\bigstar}}).

5.4. Convergence of the Newton iteration

This subsection is devoted to the proof of local convergence of the iteration introduced in (5.5). We adopt the ideas from [AR68], where operators mapping from $\mathcal{V}$ to $\mathcal{V}$ are considered. For $v\in\mathcal{V}$ and $\omega\in\mathbb{C}$ we denote by $B_{\tau}(v,\omega)\subset\mathcal{W}$ the open ball of radius $\tau$ in the $\mathcal{W}$-norm centred at $(v,\omega)$ and by $B_{\tau}(v)\subset\mathcal{V}$ the open ball of radius $\tau$ in the $\mathcal{H}$-norm centred at $v$. We base the proof on the Kantorovich theorem as stated and proved in [Deu04, Th. 2.1].

Theorem 5.3**.**

Let $\mathcal{F}_{[y]}\colon E\to\mathcal{W}^{*}$ be the continuously Fréchet differentiable mapping defined in (5.4). For a starting point $z^{0}=(u^{0},\omega^{0})\in E$ let $\partial\mathcal{F}_{[y]}(z^{0})$ be invertible. Further, we assume that $\overline{B_{\rho}(z^{0})}\subset E$ with $\rho\coloneqq(1-\sqrt{1-2h_{0}})/\kappa_{0}$ and $h_{0}\coloneqq\alpha\kappa_{0}\leq 0.5$, where $\alpha,\kappa_{0}<\infty$ are defined through

[TABLE]

for all $z,\tilde{z}\in E$. Then, the sequence $(u^{j},\omega^{j})$ obtained from iteration (5.5) is well-defined, remains in $\overline{B_{\rho}(z^{0})}$, and converges to some (u^{\scalebox{0.5}{\bigstar}},\omega^{\scalebox{0.5}{\bigstar}}) with \mathcal{F}_{[y]}(u^{\scalebox{0.5}{\bigstar}},\omega^{\scalebox{0.5}{\bigstar}})=0. For $h_{0}<0.5$, the convergence is quadratic.

In order to use Theorem 5.3, we need to determine under which conditions the Jacobian $\partial\mathcal{F}_{y}(z^{0})$ is invertible.

Lemma 5.4**.**

There exists $\tau>0$ such that

[TABLE]

and

[TABLE]

hold for any u\in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}), where \mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]}\colon\mathcal{V}^{*}\to\mathbb{C} is defined by \mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]}:=\langle\,\cdot\,,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}.

Proof.

By assumption the (geometric) multiplicity of the eigenvalue \omega^{\scalebox{0.5}{\bigstar}} is 1, which means that \mbox{dim}\big{(}\ker(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}))\big{)}=1. From Lemma 5.2 we know that \mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}) is Fredholm with index 0, which implies

[TABLE]

Further, the algebraic multiplicity of the eigenvalue \omega^{\scalebox{0.5}{\bigstar}} is also 1 by assumption, thus the Jordan chains of $\mathcal{T}$ at \omega^{\scalebox{0.5}{\bigstar}} are all of length 1 and in particular \mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u^{\scalebox{0.5}{\bigstar}}\not\in\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})), see e.g. [LGMC07, Chap. 7]. The last statement stays true in a neighbourhood of u^{\scalebox{0.5}{\bigstar}} since \mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}}) is a linear operator, i.e., there exists $\tau>0$ such that \mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u\not\in\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})) for all u\in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}). It follows that we can decompose \mathcal{V}^{*}=\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}))\oplus\operatorname{span}\{\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u\}.

In the remainder of this proof, let $u$ be any vector in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}). We note that the hermiticity of \mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}) implies that \langle\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}=\langle\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})u^{\scalebox{0.5}{\bigstar}},v\rangle_{\mathcal{V}^{*},\mathcal{V}}=0 for all $v\in\mathcal{V}$. Now, the eigenfunction u^{\scalebox{0.5}{\bigstar}}\not\equiv 0 per definition, therefore there exist some $h\in\mathcal{V}^{*}$, $\alpha\in\mathbb{C}$, and $v\in\mathcal{V}$ with h=\alpha\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u+\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v such that

[TABLE]

Therefore, \langle\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}\neq 0.

It remains to show that \ker(\mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]})=\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})). We start by showing the inclusion ”$\subseteq$”. Given any h\in\ker(\mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]})\subset\mathcal{V}^{*}, we can write h=\alpha\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u+\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v for some $\alpha\in\mathbb{C}$ and $v\in\mathcal{V}$. Thus, 0=\mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]}(h)=\langle\alpha\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}+\langle\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}=\alpha\langle\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}} implies $\alpha=0$ and therefore h\in\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})). The reverse inclusion follows from the fact that for any h\in\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})) there exists some $v\in\mathcal{V}$ such that h=\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v and thus, \mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]}(h)=\langle h,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}=\langle\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})v,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}=0. ∎

Lemma 5.5**.**

Consider u\in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}) with $\tau$ as in Lemma 5.4. Then \partial\mathcal{F}_{[y]}(u,\omega^{\scalebox{0.5}{\bigstar}}) has a bounded inverse.

Proof.

We prove the bijectivity of \partial\mathcal{F}_{[y]}(u,\omega^{\scalebox{0.5}{\bigstar}}) by showing that the equation

[TABLE]

with \omega=\omega^{\scalebox{0.5}{\bigstar}} and u\in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}) has a unique solution for all $f\in\mathcal{V}^{*}$, $g\in\mathbb{C}$. Since \mathcal{P}_{[y]}u^{\scalebox{0.5}{\bigstar}}\neq 0 by assumption, we can decompose

[TABLE]

Indeed, any $v\in\mathcal{V}$ can be written as the sum v=\alpha_{v}u^{\scalebox{0.5}{\bigstar}}+v_{\ker} with \alpha_{v}\coloneqq\left(\mathcal{P}_{[y]}u^{\scalebox{0.5}{\bigstar}}\right)^{-1}\mathcal{P}_{[y]}v\in\mathbb{C} and v_{\ker}=v-\alpha_{v}u^{\scalebox{0.5}{\bigstar}}\in\ker(\mathcal{P}_{[y]}). Therefore, the lower row of (5.6) implies g=\mathcal{P}_{[y]}s=\alpha_{s}\mathcal{P}_{[y]}u^{\scalebox{0.5}{\bigstar}} and thus \alpha_{s}=(\mathcal{P}_{[y]}u^{\scalebox{0.5}{\bigstar}})^{-1}g. Using the self-adjointness of \mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}), we find that the upper row of (5.6) tested with u^{\scalebox{0.5}{\bigstar}} yields

[TABLE]

from which it follows that

[TABLE]

since the condition \langle\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}\neq 0 holds for all u\in B_{\tau}(u^{\scalebox{0.5}{\bigstar}}). Moreover, we note that the restriction of \mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}}) to $\ker(\mathcal{P}_{[y]})$ is bijective as mapping from $\ker(\mathcal{P}_{[y]})$ to its range. Indeed,

[TABLE]

implies the injectivity and the surjectivity is obvious from the definition. Therefore we obtain from the upper row of (5.6) that

[TABLE]

which implies

[TABLE]

if f-\nu\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u\in\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})|_{\ker(\mathcal{P}_{[y]})}). In order to show that this last condition is satisfied, we use the fact that \ker(\mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]})=\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})) and the decomposition (5.7). This finally implies f-\left(\langle\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}\right)^{-1}\langle f,u^{\scalebox{0.5}{\bigstar}}\rangle_{\mathcal{V}^{*},\mathcal{V}}\mathcal{T}^{\prime}(\omega^{\scalebox{0.5}{\bigstar}})u\in\ker(\mathcal{Q}_{[u^{\scalebox{0.5}{\bigstar}}]})=\operatorname{im}(\mathcal{T}(\omega^{\scalebox{0.5}{\bigstar}})|_{\ker(\mathcal{P}_{[y]})}). ∎

Lemma 5.6**.**

Let \omega\neq\omega^{\scalebox{0.5}{\bigstar}}\in D be close enough to \omega^{\scalebox{0.5}{\bigstar}} so that $\mathcal{T}(\omega)$ is invertible, and let $u\in\mathcal{V}$ satisfy $\mathcal{P}_{[y]}\mathcal{T}(\omega)^{-1}\mathcal{T}^{\prime}(\omega)u\neq 0$. Then $\partial\mathcal{F}_{[y]}(u,\omega)$ has a bounded inverse.

Proof.

Note that there exists a neighborhood around \omega^{\scalebox{0.5}{\bigstar}} in which $\mathcal{T}(\omega)$ is invertible, since the spectrum of $\mathcal{T}$ consists of isolated eigenvalues. We prove the bijectivity of $\partial\mathcal{F}_{[y]}(u,\omega)$ by showing that (5.6) has a unique solution for all $f\in\mathcal{V}^{*},g\in\mathbb{C}$.

The upper row of (5.6) implies

[TABLE]

whereas the lower row implies

[TABLE]

Thus, $\nu\in\mathbb{C}$ is uniquely defined by

[TABLE]

as long as the condition $\mathcal{P}_{[y]}\mathcal{T}^{-1}(\omega)\mathcal{T}^{\prime}(\omega)u\neq 0$ is satisfied. As a result, $\partial\mathcal{F}_{[y]}(u,\omega)$ is linear, bijective, and bounded, which implies the existence of a bounded inverse. ∎

The previous lemmata show that choosing $(u^{0},\,\omega^{0})$ close enough to (u^{\scalebox{0.5}{\bigstar}},\omega^{\scalebox{0.5}{\bigstar}}) in $E$ with

[TABLE]

insures that $\partial\mathcal{F}(u^{0},\omega^{0})$ has a bounded inverse. In this case, the first iteration step of (5.5) is well-defined. An estimate for the bound $\kappa_{0}$ appearing in Theorem 5.3 can be derived from the fact that $\mathcal{F}_{[y]}$ is twice continuously Fréchet differentiable and that $\partial^{2}\mathcal{F}_{[y]}(u,\omega)$ is bounded for all $(u,\omega)\in E$, since

[TABLE]

To sum up, the results of this section show the local convergence of the Newton iteration (5.5) under the assumption that the starting point $(u^{0},\,\omega^{0})\in E$ is close enough to the eigenpair and satisfies the condition (5.8).

6. Numerical Experiments

In this final section we consider the proposed PDE eigenvalue iterations numerically. The first two examples in Sections 6.1 and 6.2 consider the application of the inverse power method. Here, the essential idea is to combine the iterative solver and mesh refinement. Since we have proven the convergence on operator level, we can choose the spatial discretization in each step independently. Thus, we can start the computations on coarse grids and consider refinements as the eigenvalue iteration goes on, leading to substantial runtime improvements. Finally, Section 6.3 shows results for the application of the Newton-type iteration, which was analyzed in Section 5.

In all three experiments the computational domain equals the unit square $\Omega=(0,1)^{2}$ with a disk of radius $0.3$ in the middle defining $\Omega_{2}$. The outer material is air, i.e., $\varepsilon_{1}(\omega)\equiv\alpha_{1}=1$, whereas the relative electrical permittivity of the inner material varies. The spatial discretization has been performed with the software library Concepts [FL02, Con19] using piecewise polynomial finite element basis functions of degree $2$ on a quadrilateral mesh with curved cells. We refer to [SK09, SK10] for more details on the application of finite elements to eigenvalue problems arising in photonic crystals simulations using Concepts.

6.1. Inverse power method for linear model

We consider the problem described in Example 2.8 for a fixed wave vector ${\mathbf{k}}=[\pi/2,\pi]^{T}$ and a constant model $\varepsilon_{2}(\omega)\equiv 8$ for the relative permittivity in $\Omega_{2}$. In order to find the smallest eigenvalue, we take as shift $\beta=1$ and apply the Rayleigh quotient iteration (3.2).

The obtained convergence history for the eigenvalue is shown in the top plot of Figure 6.1. More precisely, we show the relative error $|\mu^{j}-{\mu^{\scalebox{0.5}{\text{ref}}}}|/|{\mu^{\scalebox{0.5}{\text{ref}}}}|$, where $\mu^{\scalebox{0.5}{\text{ref}}}$ denotes the reference solution (as approximation of \mu^{\scalebox{0.5}{\bigstar}}) obtained on a finer spatial mesh with a large number of inverse iteration steps. Here, the reference mesh corresponds to the sixth level of uniform refinement. Note that \mu^{\scalebox{0.5}{\text{ref}}}>\mu^{\scalebox{0.5}{\bigstar}} due to the spatial discretization error. We compare the results for $3$, $5$, and $7$ iteration steps per mesh. This means that we perform a fixed number of iteration steps before we consider a uniform refinement of the mesh. When the fine mesh (corresponding here to the fifth level of refinement) is reached, the inverse power iterations continue without further mesh refinement. The corresponding convergence on this fine mesh is depicted by the dashed line in Figure 6.1.

One observation is that the results for $3$ iteration steps ( ) per mesh are comparable to the results obtained on the fine mesh. We emphasize that the first $15$ iteration steps are roughly as expensive as a single step on the fine mesh but leading to a much better accuracy. Further note that the saturation of the convergence of the eigenvalue is due to the discretization error, since $\mu^{\scalebox{0.5}{\text{ref}}}$ is computed on a finer spatial mesh than the iterates $\mu^{j}$. The computational gain is even more distinct for the convergence of the residuum (and thus the eigenfunction), which is shown in the bottom plot of Figure 6.1. Here, the residuum $\operatorname{Res}\colon\mathcal{V}\times\mathbb{C}\to\mathcal{V}^{*}$ is defined by

[TABLE]

and measured in the dual norm, i.e.,

[TABLE]

In the discrete setting, the dual norm is given by $\|\cdot\|_{M(K+M)^{-1}M}$. The results for $5$ ( ) and $7$ ( ) iteration steps per mesh exhibit a similar behavior, i.e., a given tolerance is reached faster, since the first iteration steps are performed on coarse meshes, cf. Table 6.1.

6.2. Inverse power method for nonlinear model

Motivated by the numerical example in [EKE12] we consider a material within the disk $\Omega_{2}$ having a frequency-dependent permittivity described by a two-term Lorentz model, cf. (4.2), with positive and real parameters

[TABLE]

Recall that we consider here relative permittivities as discussed in Section 2.2. As in the previous example, we are interested in the lower-most eigenvalue of the corresponding nonlinear eigenvalue problem (4.1) for the fixed wave vector ${\mathbf{k}}=[\pi/2,\pi]^{T}$.

For the numerical solution of the eigenvalue problem, we apply the inverse power method (4.9) to the linearized system with shift $\beta=\eta_{2}^{2}-\eta_{1}^{2}+1$. This shift guarantees the positivity of the involved operator, cf. Lemma 4.8. In Figure 6.2 we compare again the results for a fixed number of iteration steps per mesh. Here, the computational results are even more convincing as the plots for $3$ and $7$ iteration steps per mesh are very similar to the results obtained on the fine mesh, although the computational costs are significantly smaller. Note that the iterations are based on the linearization (4.7) but the residuals are computed in terms of the nonlinear eigenvalue problem (4.1).

A comparison of the runtimes, showing the speed-up due to the combination of the iterative solver and mesh refinement, is shown in Table 6.2.

Finally, the eigenmode corresponding to the first eigenvalue obtained by the inverse power method (4.9) applied to the linearized system is depicted in Figure 6.3.

6.3. Newton method for full Drude-Lorentz model

In this final experiment we consider 66% porous silicon within the disk $\Omega_{2}$ and model its relative permittivity with a function of the form (5.1) using $\alpha_{\text{DL}}=1.143$ and the material parameters given in Table 6.3.

We apply the Newton iteration (5.5) with a randomly chosen normalizing vector $y$ to the nonlinear eigenvalue problem (5.3) for the fixed wave vector ${\mathbf{k}}=[\pi/2,\pi]^{T}$. Since this eigenvalue problem depends only on the square of $\omega$, we set $\lambda\coloneqq\omega^{2}$ and actually consider $\lambda$ to be the eigenvalue instead of $\omega$.

In order to produce starting values $u^{0}$ and $\lambda^{0}$ for the Newton iteration, we approximate $\varepsilon_{\text{DL}}^{\Re}(\omega)$ by the constant value $2$ in the disk $\Omega_{2}$ and apply 8 Rayleigh quotient iteration steps, similarly as in Section 6.1.

The resulting convergence history of the eigenvalue and the residual are shown in Figure 6.4. One can observe that already 6 iterations on the fine mesh (corresponding to the fourth level of refinement) are sufficient to converge to an eigenpair for which the residual is of the order of magnitude of the machine precision. With such a fast converging method, a mesh refinement should be applied after no more than 2 iterations per mesh.

7. Conclusion

In this paper, we have considered iterative methods for linear Hermitian as well as specific nonlinear eigenvalue problems arising in photonic crystal modeling. In case the electric permittivity is given by a Drude-Lorentz model with real coefficients and no dissipation, we are able to linearize the problem to obtain a linear and Hermitian eigenvalue problem. For this, we show the convergence of the inverse power method.

For more realistic models taking dissipation into account, the same procedure would lead to a linear but non-Hermitian eigenvalue problem. Thus, instead of a linearization we directly apply Newton’s method, for which we prove local convergence.

Acknowledgements

The authors thank Christoph Zimmer (TU Berlin) for the helpful discussions on the paper.