---
abstract: 'We propose an algorithm for general nonlinear eigenvalue problems to compute eigenvalues within a chosen contour and to compute the corresponding eigenvectors. Eigenvalue information is explored by contour integration incorporating different weight functions. The gathered information is processed by solving a nonlinear system of equations of small dimension. No auxiliary functions have to be introduced for linearization. The numerical implementation of the approach is straightforward and the algorithm allows for parallelization. We apply the method to two examples from physics. Resonant states of a one-dimensional quantum mechanical system and resonant states of a three-dimensional photonic nanoantenna are computed.'
author:
- 'Felix Binkowski[^1]'
- 'Lin Zschiedrich[^2]'
- Sven Burger
title: |
    A Riesz-projection-based method for\
    nonlinear eigenvalue problems
---

nonlinear eigenvalue problems, contour integration, Riesz projection, photonic nanoantenna, resonant states

Introduction {#sec:Intro}
============

The numerical treatment of nonlinear eigenvalue problems (NLEVPs) is a highly relevant research field in applied mathematics [@Guettel_NLEVP_2017; @Mackey_2015; @Mehrmann_GAMM_2004]. Fundamental solution techniques from numerical linear algebra which are used for solving linear eigenproblems are not available in the case of NLEVPs. This leads to challenges for the development of suitable algorithms. We address the most general problem class of NLEVPs $$\begin{aligned}
    T(\lambda)v = 0, \label{eq:NLEVP}\end{aligned}$$ where $T(\lambda) \in \mathbb{C}^{n \times n}$ is the eigenvector residual function, is an eigenvalue and is an eigenvector corresponding to $\lambda$.

In physics, NLEVPs occur in several fields, from the dynamic analysis of macroscopic structures [@Lancaster_BookVibrSys_2002; @Przemieniecki_BookMatrixStruct_2012] to the investigation of photonic resonators on the nanoscale . NLEVPs often have the form $ A(\lambda)v = \lambda B(\lambda)v,  $ which can be brought into the above form with $T(\lambda) = A(\lambda) - \lambda B(\lambda)$. In many applications, $T(\lambda)$ has a sparse matrix structure and is very large, while just a few eigenpairs $(\lambda,v)$ are responsible for the physical behavior of the described problem. For rational residual functions $T(\lambda)$, the NLEVP can be cast into a linear form. Such a linearization introduces auxiliary functions and increases the dimension of the problem [@Mehrmann_GAMM_2004; @Tisseur_SIAMrev_2001]. Material dispersion is often significant in physical systems and modeled by measured material data. To apply the approach of linearization, material data have to be fitted by rational functions and numerical costs grow with the order and number of poles of the fit.

Standard solvers for NLEVPs contain iterative projection methods [@Saad_Book_NumMeth_Eig_2011], such as the Arnoldi or the Jacobi-Davidson method, which take advantage of the sparse matrix structures. These methods are applied either to the linearized systems or to the NLEVPs directly [@Voss_nonlinArnoldi_2004; @Voss_nonlinJacobiDavid_2007]. For the computation of one or more specific eigenpairs, the choice of an appropriate initial guess is essential [@Mehrmann_GAMM_2004]. A further problematic task is to avoid convergence to the same eigenpair repeatedly. In the case of linear eigenproblems, this can be ensured by a Schur decomposition, which is not available in the nonlinear case [@Voss_nonlinJacobiDavid_2007]. Contour integral methods for solving NLEVPs [@Asakura_JSIAM_2009; @Beyn_LAAppl_2012; @Gavin_JCompPhy_2018; @Yokota_JSIAM_2013] address these challenges by the construction of an approximate subspace corresponding to specific eigenvectors whose eigenvalues are located inside a chosen region in the complex plane. The subspace is computed by contour integration where typically random vectors are projected onto the eigenspace of interest. In [@Gavin_JCompPhy_2018; @Yokota_JSIAM_2013], the Rayleigh-Ritz method is then used to achieve approximations of eigenpairs. The nonlinear structure is still inherited to the lower dimensional projected system and needs to be given explicitly in a rational form. The methods of [@Asakura_JSIAM_2009; @Beyn_LAAppl_2012] extract eigenvalue and eigenvector information by applying a singular-value decomposition to the subspace generated by the contour integration and then solving a linear eigenproblem. An alternative way of extracting this information is based on canonical polyadic tensor decomposition [@vanBarel_ContInt_JoCompApplMath_2016].

In this work, we present a numerical method for solving NLEVPs which circumvents linearization in any stage of the procedure. Data for the eigenvalues within a chosen contour are generated by contour integration incorporating different weight functions and fitted to a nonlinear model based on Cauchy’s residue theorem. Essentially, the algorithm only relies on solving linear systems of equations $T(\lambda)^{-1}y$ along the contour, where $y$ is a random vector. The method is algorithmically simple as these systems can be regarded as a blackbox. Any material dispersion can be included.

Method for solving NLEVPs {#sec:RP}
=========================

This section derives an approach to compute eigenpairs $(\lambda,v)$ fulfilling where the eigenvalues are located inside a chosen contour. To start with, notation and theoretical background on elements of complex analysis is introduced [@Nolting_BookElectrodyn3_2016 Section 4.4].

We consider with a regular matrix function , where $\Omega \subset \mathbb{C}$. Let $x,y \in \mathbb{C}^n$ be random vectors. Let $\Gamma_k \subset \Omega$ be a contour which encloses one eigenvalue $\lambda_k$ of the residual function $T(\lambda)$ and on which the function $x^HT(\lambda)^{-1}y$ is holomorphic. The eigenvalue $\lambda_k$ is a pole of $x^HT(\lambda)^{-1}y$ and the pole is assumed to be of order $p$. Then, the Laurent series for $x^HT(\lambda)^{-1}y$ about $\lambda_k$ is given by $$\begin{aligned}
    x^H T(\lambda)^{-1}y &= \sum_{n=-p}^{\infty} a_n(\lambda-\lambda_k)^n,\hspace{0.5cm}
    a_{n}(\lambda_k) := \frac{1}{2 \pi i} \oint \limits_{\Gamma_k}
    \frac{x^H T(\xi)^{-1}y}{(\xi-\lambda_k)^{n+1}}\,
    d\xi  \in \mathbb{C}. \label{eq:Laurent}\end{aligned}$$ The coefficient $a_{-1}(\lambda_k)$ is the so-called residue $$\begin{aligned}
    {\text{Res}}_{\lambda_k}\left(x^H T(\lambda)^{-1}y\right) =
    \frac{1}{2 \pi i} \oint \limits_{\Gamma_k} x^H T(\lambda)^{-1}y
    \, d\lambda \label{eq:residuum}\end{aligned}$$ of $x^H T(\lambda)^{-1}y$ at $\lambda_k$.

Sketch of the approach {#sec:RP_sketch}
----------------------

To show the idea of this work, a simple eigenvalue $\lambda_k$ of $T(\lambda)$ is assumed, i.e., $\lambda_k$ is a pole of $x^HT(\lambda)^{-1}y$ and has the order $p=1$. With the aim of extracting eigenvalue information from , the scalar function $f(\lambda) = \lambda$ is introduced. Then, Cauchy’s integral formula leads to $$\begin{aligned}
    {\text{Res}}_{\lambda_k} \,\left(\lambda x^H T(\lambda)^{-1}y\right) &= 
    \frac{1}{2 \pi i} \oint \limits_{\Gamma_k} \lambda x^H T(\lambda)^{-1}y
    \, d\lambda   \nonumber \\  &=
    \frac{1}{2 \pi i} \oint \limits_{\Gamma_k} \frac{ \lambda}{\lambda-\lambda_k}
    a_{-1}(\lambda_k)
    \, d\lambda \nonumber \\  &= 
    \lambda_k {\text{Res}}_{\lambda_k}\,\left(x^H T(\lambda)^{-1}y\right), \nonumber\end{aligned}$$ where the regular part of the Laurent series vanishes and only the principal part remains. Rearranging yields the eigenvalue $$\begin{aligned}
    \lambda_k = 
 \frac{     {\text{Res}}_{\lambda_k}\left(\lambda x^H T(\lambda)^{-1}y\right)}
    {   {\text{Res}}_{\lambda_k}\left(x^H T(\lambda)^{-1}y\right)}\nonumber\end{aligned}$$ inside of the contour $\Gamma_k$. The eigenvector $v_k$ corresponding to $\lambda_k$ can be obtained by the Riesz projection $$\begin{aligned}
    P \coloneqq \frac{1}{2 \pi i} \oint \limits_{\Gamma_k} T(\lambda)^{-1}
    \, d\lambda \nonumber\end{aligned}$$ for $T(\lambda)$ and $\Gamma_k$, which projects vectors onto the eigenspace associated with the eigenvalue inside of $\Gamma_k$ [@Hislop_BookSpecTheoSchroed_1996]. The eigenvector is given by $v_k =Py$.

Generalized approach {#sec:numeric}
--------------------

The idea of the previous subsection can be generalized. Firstly, we assume that the pole $\lambda_k$ has an order $p \geq 1$ and consider a function $f: \Omega \rightarrow \mathbb{C}$ which is holomorphic on $\Gamma_k$ and inside of $\Gamma_k$ yielding $$\begin{aligned}
\frac{1}{2 \pi i} \oint \limits_{\Gamma_k} f(\lambda) x^H T(\lambda)^{-1}y
\, d\lambda   &=
\sum_{n=-p}^{-1} a_{n}(\lambda_k)\frac{1}{2 \pi i} \oint \limits_{\Gamma_k}
f(\lambda)(\lambda-\lambda_k)^{n}   \, d\lambda  \label{eq:one_pole} \\ &= 
\sum_{n=-p}^{-1} a_{n}(\lambda_k) \frac{f(\lambda_k)^{(-n-1)}}{(-n-1)!}. \nonumber\end{aligned}$$ The coefficients $a_{-p}(\lambda_k),\dots,a_{-1}(\lambda_k)$, are from the Laurent series for $x^HT(\lambda)^{-1}y$ in . Secondly, we choose a contour on which $x^HT(\lambda)^{-1}y$ is holomorphic and which encloses finitely many poles $\lambda_1, \dots, \lambda_m$ of $x^HT(\lambda)^{-1}y$. Cauchy’s residue theorem is used [@Nolting_BookElectrodyn3_2016] so that can be extended to $$\begin{aligned}
    &\frac{1}{2 \pi i} \oint \limits_{\Gamma} f(\lambda) x^H T(\lambda)^{-1}y
    \, d\lambda   = \sum_{j=1}^{m}
    \sum_{n=-p}^{-1} a_{n}(\lambda_j) \frac{f(\lambda_j)^{(-n-1)}}{(-n-1)!},
     \label{eq:residue_theorem_many_poles}\end{aligned}$$ where $f(\lambda)$ has to be holomorphic on $\Gamma$ and inside of $\Gamma$. contains the eigenvalues $\lambda_1, \dots, \lambda_m$ of the residual function $T(\lambda)$. To explore the information for these $m$ unknowns, we introduce several *weight functions* $f_1(\lambda), \dots, f_M(\lambda)$ and construct the nonlinear system of equations (NLSE) $$\begin{aligned}
    \mu_k = F_k(\lambda_1,\dots,\lambda_m),
    \hspace{0.5cm} k=1,\dots,M, \label{eq:nonlin_sys}\end{aligned}$$ where $$\begin{aligned}
    \mu_k &\coloneqq \frac{1}{2 \pi i} \oint \limits_{\Gamma}
    f_k(\lambda)x^H T(\lambda)^{-1}y
    \, d\lambda,\nonumber \\
    F_k(\lambda_1,\dots,\lambda_m) &=  \sum_{j=1}^{m}
    \sum_{n=-p}^{-1} a_{n}(\lambda_j) \frac{f_k(\lambda_j)^{(-n-1)}}{(-n-1)!}.
    \nonumber\end{aligned}$$ Solving this NLSE yields the eigenvalues inside of the contour $\Gamma$. The corresponding eigenvectors $v_1,\dots,v_m$ are given by $$\begin{aligned}
    v_k = \frac{1}{2 \pi i} \oint \limits_{\Gamma}
    h_k(\lambda)T(\lambda)^{-1}y \, d\lambda, \hspace{0.5cm}
    k=1,\dots,m,
    \label{eq:filter_func}\end{aligned}$$ where $$\begin{aligned}
    h_k(\lambda) = \prod_{\substack{j = 1 \\ j \neq k}}^m
    \left(\lambda - \lambda_j\right)^{p_j}
    \nonumber\end{aligned}$$ are *filter functions* and $p_j$ is the order of the corresponding pole. The Riesz projection is applied for the function $h_k(\lambda)T(\lambda)^{-1}$ where certain eigenvalues are filtered out.

Algorithm
---------

To solve NLEVPs based on the presented theoretical considerations, we propose . The algorithm can be sketched as follows. A contour $\Gamma$ which encloses the eigenvalues of interest has to be chosen (Step 1). The contour integrals in are computed with a suitable quadrature rule (Step 2-4). The evaluation of the integrand at the integration points essentially requires to solve linear systems of equations $T(\lambda)^{-1}y$. The results can be stored and reused to compute the eigenvectors.

The calculated integrals serve as the input data for solving the NLSE given by (Step 5). The nonlinear model $F_k(\lambda_1,\dots,\lambda_{{m}})$ is chosen based on the expected order of the poles and based on the expected number of poles inside of the selected contour $\Gamma$. The number $m$ of unknowns for the nonlinear model has to be greater or equal than the number of eigenvalues inside of $\Gamma$. To solve the NLSE, an initial guess is required. If no a priori information about the eigenvalues is available, then randomly chosen numbers inside of the contour $\Gamma$ are a possible choice.

The result of solving the NLSE are approximations to the eigenvalues inside of $\Gamma$. Approximations of the corresponding eigenvectors can be achieved by computing the contour integrals in using the calculated eigenvalues (Step 6-7). The algorithm returns approximate eigenvalues $\lambda_1,\dots,\lambda_m$ and eigenvectors $v_1,\dots,v_m$ (Step 8).

is parallelizable on two levels. Firstly, the complete algorithm can be performed for different contours $\Gamma$ simultaneously. This can be useful if eigenvalues in different regions are of interest. Secondly, solving the linear systems for the numerical integration can be done in parallel.

eigenvalues $\lambda_1,\dots,\lambda_m$, eigenvectors $v_1,\dots,v_m$

Numerical realization {#sec:numer_real}
---------------------

We realize the numerical integration in with an $N$-point trapezoidal rule. The integration path $\Gamma$ is a circular contour, which leads to exponential convergence [@Trefethen_SIAM_Trapz_2014]. The equidistant integration points are given by $\hat{\lambda}_k = \lambda_0 + r e^{2 \pi i k/N}$, , where $\lambda_0$ and $r$ are the center and the radius of $\Gamma$, respectively. Note that recently, methods based on rational filters for computing eigenvalues of real symmetric matrices have been developed [@Austin_RatFilt_RealSymm_2015] and rational filter functions for contour integral discretization have been designed using optimization techniques [@vanBarel_ratFilt_LinAlgAppl_2016]. However, for the sake of simplicity, we use a trapezoidal rule. To solve the linear systems of equations $T(\hat{\lambda}_k)^{-1}y$, , an LU decomposition can be used. The LU decomposition needs not to be updated at each integration point. Instead, for sufficiently small changes in $\lambda$, the LU decomposition of a previous evaluation can be used as a preconditioner for iterative solving. This leads to a more efficient numerical implementation.

In the following section, we consider physical examples with simple eigenvalues, i.e., the nonlinear model in simplifies to $$\begin{aligned}
F_k(\lambda_1,\dots,\lambda_m) &=  \sum_{j=1}^{m}
a_{-1}(\lambda_j) f_k(\lambda_j).
\nonumber\end{aligned}$$ is solved with the nonlinear solver *fsolve* from MATLAB. We regard the residues $a_{-1}(\lambda_1),\dots,a_{-1}(\lambda_m)$ as unknown variables themselves and set $M=2m$ to construct non-underdetermined NLSEs. This handling of the residues allows for a simpler numerical realization. If $m$ is greater than the number of eigenvalues inside of $\Gamma$, then *fsolve* returns results where some residues are very small. These results can be discarded as they do not correspond to eigenvalues. For the weight functions, we choose the scaled polynomials $$\begin{aligned}
f_k(\lambda) = \left(\frac{\lambda - \lambda_0}{r}\right)^{k-1},
\hspace{0.5cm} k = 1,\dots,M. \nonumber\end{aligned}$$ Due to the fact that the weight functions are known and that we treat the residues as unknowns, also the Jacobians can be provided for the nonlinear solver.

Application of the method {#sec:Appl}
=========================

We apply \[alg:RP\_method\] to NLEVPs resulting from the Schrödinger equation and Maxwell’s equations. These quantum mechanical and quantum optical examples are open systems, which are described by non-Hermitian operators. In physics, the eigensolutions of such problems are usually called resonant states or quasinormal modes [@Lalanne_QNMReview_2018; @Zworski_Resonances_1999]. Material dispersion is omnipresent in such systems and the physical understanding of the resonance phenomena through numerical simulations is an active research topic [@Franke_QNMS_QED_2018; @Yan_PRB_2018].

Resonant states in an open quantum system {#sec:Appl1}
-----------------------------------------

Propagation of a quantum particle of effective mass $m^*$ through a one-dimensional potential $V(x)$ can be described by the time-independent Schrödinger equation $$\begin{aligned}
-\frac{\hbar}{2} \frac{\partial}{\partial x}\left( \frac{1}{m^*}
\frac{\partial \Psi(x)}{\partial x}  \right) + V(x) \Psi(x)
= E \Psi(x), \label{eq:Schroedinger}\end{aligned}$$ where $\hbar$ is the Planck constant, $E$ is the energy and $\Psi(x)$ is the unknown wave function. For a detailed description and motivation of this example, see [@Shao_JapplPhys_1995]. To compute the resonant states of this problem in the domain $[-L,L]$, we use the approach of [@Gavin_JCompPhy_2018] to scale and discretize yielding the quadratic NLEVP $$\begin{aligned}
    T(\lambda)v = \left(\lambda^2A_2 + i\lambda A_1 - A_0\right)
    v = 0, \nonumber\end{aligned}$$ where $$\begin{aligned}
A_2 =
\frac{h}{6}\begin{bmatrix}
2       & 1     & 0       & 0        & \dots & 0        \\
1       & 4     & 1       & 0        & \dots & 0        \\
0       & 1     & 4       & 1        & \dots & 0        \\
\vdots  &       &         &\ddots    &       & \vdots   \\
0       &\dots  & 0       & 1        & 4     & 1        \\
0       &\dots  & 0       & 0        & 1     & 2        
  \end{bmatrix}, \hspace{0.25cm}
A_1 =
\begin{bmatrix}
1       & 0       & \dots   & 0          & 0        \\
0       & 0           & \dots   & 0      & 0        \\
\vdots  & \vdots& \vdots  & \vdots       & \vdots       \\
0       & 0           & \dots   & 0      & 0        \\
0       & 0           & \dots   & 0      & 1        \\  
\end{bmatrix}, \nonumber\end{aligned}$$ $$\begin{aligned}
A_0 =
\frac{1}{h}\begin{bmatrix}
1       & -1    & 0     & 0     & \dots & 0         \\
-1      & 2     & -1    & 0     & \dots & 0         \\
0       & -1    & 2     & -1    & \dots & 0         \\
\vdots  &       &       &\ddots &       & \vdots    \\
0       &\dots  & 0     & -1    & 2     & -1        \\
0       &\dots  & 0     & 0     & -1    & 1     
\end{bmatrix} - V_0 A_2 \in \mathbb{R}^{n+2\times n+2}. \nonumber\end{aligned}$$ The same parameters as in [@Gavin_JCompPhy_2018] are chosen, where $L=\pi/ \sqrt{2}$, and the spatial step size is $h = 2L/(n+1)$ with $n=302$.

We perform a set of computations in order to demonstrate a numerical realization of \[alg:RP\_method\]. For the calculation of certain eigenpairs, a circular contour $\Gamma$ with the center $\lambda_0=2$ and the radius $r=3$ is considered. Solutions for an increasing number $N$ of integration points and for different numbers of unknown eigenvalues, $m=7,\dots,11$, are computed. The maximum residual $\max_k||T(\lambda_k)v_k||_2$, where $||v_k||_2 = 1$, is shown in . Exponential convergence up to a certain accuracy is obtained. It can be further observed that for small $N$ and an increasing $m$, the residuals become smaller. Here, the nonlinear solver can exploit more degrees of freedom and thus compute more accurate solutions. The results of \[alg:RP\_method\] can be compared with solutions of the linearized problem [@Gavin_JCompPhy_2018] computed with MATLAB’s function *eigs*. These reference solutions and the eigenvalues computed by \[alg:RP\_method\] with $N=100$ and $m=10$ are shown in . Seven eigenvalues inside of $\Gamma$ are obtained by solving the linearized system. The eigenvalues resulting from \[alg:RP\_method\] coincide with these reference solutions.

Photonic nanoantenna {#sec:Appl2}
--------------------

In the second numerical experiment, we consider a nanooptical structure. Nanoantennas allow, e.g., for realizing single-photon emitters for quantum technology devices [@Koenderink_SinglePhoton_2017]. We apply \[alg:RP\_method\] to an example from [@Zschiedrich_PRA_2018], where a defect in a diamond nanodisk is considered as solid-state single-photon emitter. In the steady-state regime, the light-matter interaction of such a structure can be described by the time-harmonic Maxwell’s equations in the second-order form $$\begin{aligned}
    \nabla \times \mu(\mathbf{r},\omega)^{-1} \nabla \times \mathbf{E}(\mathbf{r},\omega) -
    \omega^2\epsilon(\mathbf{r},\omega) \mathbf{E}(\mathbf{r},\omega) =
    i\omega\mathbf{J}(\mathbf{r}) \label{eq:Maxwells_eq}\end{aligned}$$ with the electric field $\mathbf{E}(\mathbf{r},\omega) \in \mathbb{C}^3$ and the source term $\mathbf{J}(\mathbf{r}) \in \mathbb{C}^3$ as impressed current, where $\mathbf{r} \in \mathbb{R}^3$ is the position. The permittivity tensor $\epsilon(\mathbf{r},\omega)$ characterizes the spatial distribution of materials and, through its dependence on the complex angular frequency $\omega \in \mathbb{C}$, the material dispersion. In the regime of optical frequencies, the permeability tensor $\mu(\mathbf{r},\omega)$ can typically be set to the vacuum permeability $\mu_0$. Here, we are interested in solutions to in absence of the source term $\mathbf{J}(\mathbf{r})$.

To compute eigenvalues of the diamond nanostructure, is slightly adapted. At each integration point, is solved for $i\omega \mathbf{J}(\mathbf{r})$ depending on $\omega$. This solution of scattering problems replaces Step 3 in . We integrate the scattering solutions with weight functions along a contour in the and evaluate the resulting fields at one chosen point. This emphasizes the blackbox character of \[alg:RP\_method\]. Step 3 is essentially replaced by applying $\mathcal{L}\left( T(\lambda)^{-1}y\right)$ along the contour, where $\mathcal{L}$ is a linear operator evaluating the solution fields at the chosen point. The resulting eigenvalue is a frequency , where $c$ is the speed of light and $\lambda_k$ is a wavelength. The computations are performed using the FEM solver JCMsuite. The outgoing radiation conditions for the diamond nanoresonator are realized with perfectly matched layers [@Berenger_1994; @Hohage_PoleCondPartII_2003]. We refer to [@Pomplum_NanoopticFEM_2007; @Zschiedrich_PRA_2018] for details on the FEM implementation and for details on the diamond nanoresonator.

For this numerical experiment, a circular contour $\Gamma$ is considered where the center is and the radius is . We apply \[alg:RP\_method\] while setting $m=5$ for solving the NLSE, which yields three eigenvalues inside of $\Gamma$. The maximum relative error of their real and imaginary parts with respect to the number $N$ of integration points is shown in . The reference solutions are computed with the eigensolver from , which applies an Arnoldi algorithm using auxiliary fields. Here, we choose as initial guess. We observe convergence to the reference solutions up to a relative error of about $10^{-8}$. The accuracy limitation can be attributed to the accuracy of the scattering problem solver from JCMsuite. The calculated eigenvalues of and the reference solutions are shown in .

Conclusions and outlook {#sec:Concl}
=======================

We have presented an algorithm based on contour integration for solving general NLEVPs. A numerical realization has been applied to non-Hermitian problems from the fields of quantum mechanics and photonics. We have considered the numerical solution of $T(\lambda)^{-1}y$ as a blackbox and extracted eigenvalue information by introducing the meromorphic function $x^HT(\lambda)^{-1}y$. This can be generalized by $\mathcal{L}\left( T(\lambda)^{-1}y\right)$, where the linear operator $\mathcal{L}$ can be a physical observable, e.g., a point evaluation as in . Instead of global eigenfunctions, modal contributions of the linear observable, e.g., modal Purcell factors [@Zschiedrich_PRA_2018], can be computed. In this way, eigenvalues can still be extracted without the need of a global approximation of the solution field $T(\lambda)^{-1}y$. Therefore, we expect that the algorithm will prove especially useful for numerical methods without a vector representation of the solution field, such as semi-analytical approaches.

Recently, the here presented method has been compared with standard eigensolvers for computing resonant states [@Lalanne_QNM_Benchmark_2018]. Future research will address the numerical treatment of the NLSE in where the impact of the choice of the weight functions is an interesting question. Further, the effect of alternative quadrature formulas will be investigated.

Acknowledgments {#acknowledgments .unnumbered}
===============

We acknowledge support from Einstein Foundation Berlin within the framework of MATHEON (ECMath project OT9).

[10]{}

, [*[A numerical method for nonlinear eigenvalue problems using contour integrals]{}*]{}, JSIAM Lett., 1 (2009), pp. 52–55.

, [*[Computing Eigenvalues of Real Symmetric Matrices with Rational Filters in Real Arithmetic]{}*]{}, SIAM J. Sci. Comput., 37 (2015), pp. A1365–A1387.

, [*[Designing rational filter functions for solving eigenvalue problems by contour integration]{}*]{}, [Linear Algebra Its Appl.]{}, 502 (2016), pp. 346–365.

, [*[Nonlinear eigenvalue problems and contour integrals]{}*]{}, J. Comput. Appl. Math., 292 (2016), pp. 526–540.

, [*[A Perfectly Matched Layer for the Absorption of Electromagnetic Waves]{}*]{}, J. Comput. Phys., 114 (1994), pp. 185–200.

, [*[An integral method for solving nonlinear eigenvalue problems]{}*]{}, Linear Algebra Its Appl., 436 (2012), pp. 3839 – 3863.

, [*[Optical Antennas]{}*]{}, Adv. Opt. Photon., 1 (2009), pp. 438–483.

, [*[Quantization of quasinormal modes for open cavities and plasmonic cavity-QED]{}*]{}, preprint, arXiv:1808.06392 (2018).

, [*[FEAST eigensolver for nonlinear eigenvalue problems]{}*]{}, J. Comput. Sci, 27 (2018), pp. 107–117.

, [*[The Nonlinear Eigenvalue Problem]{}*]{}, Acta Numer., 26 (2017), pp. 1–94.

, [*[Introduction to spectral theory. With applications to Schrödinger operators]{}*]{}, Springer, New York, 1996.

, [*[Solving Time-Harmonic Scattering Problems Based on the Pole Condition II: Convergence of the PML Method]{}*]{}, SIAM J. Math. Anal., 35 (2003), pp. 547–560.

, [*[Single-Photon Nanoantennas]{}*]{}, ACS Photonics, 4 (2017), pp. 710–722.

, [*Quasinormal mode solvers for resonators with dispersive materials*]{}, preprint, arXiv:1811.11751 (2018).

, [*[Light Interaction with Photonic and Plasmonic Resonances]{}*]{}, Laser Photonics Rev., 12 (2018), 1700113.

, [*[Lambda-matrices and Vibrating Systems]{}*]{}, Dover Publications, Mineola (New York), 2002.

, [*Polynomial eigenvalue problems: Theory, computation, and structure*]{}, in Numerical Algebra, Matrix Theory, Differential-Algebraic Equations and Control Theory, Springer, Cham, 2015, pp. 319–348.

, [*[Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods]{}*]{}, GAMM Mitt., 27 (2005), pp. 121–152.

, [*[Theoretical physics 3. Electrodynamics]{}*]{}, Springer, Cham, 2016.

, [*[Adaptive finite element method for simulation of optical nano structures]{}*]{}, Phys. Status Solidi B, 244 (2007), pp. 3419–3434.

, [*[Theory of Matrix Structural Analysis]{}*]{}, Dover Publications, Mineola (New York), 2012.

, [*[Numerical Methods for Large Eigenvalue Problems]{}*]{}, SIAM, Philadelphia, 2011.

, [*[An eigenvalue method for open-boundary quantum transmission problems]{}*]{}, J. Appl. Phys, 78 (1995), pp. 2177–2186.

, [*[The Quadratic Eigenvalue Problem]{}*]{}, SIAM Rev., 43 (2001), pp. 235–286.

, [*[The Exponentially Convergent Trapezoidal Rule]{}*]{}, SIAM Rev., 56 (2014), pp. 385–458.

, [*[[Optical microcavities]{}]{}*]{}, Nature, 424 (2003), pp. 839–846.

, [*[An Arnoldi Method for Nonlinear Eigenvalue Problems]{}*]{}, [BIT Numer. Math.]{}, 44 (2004), pp. 387–401.

, [*[A Jacobi-Davidson method for nonlinear and nonsymmetric eigenproblems]{}*]{}, Comput. Struct., 85 (2007), pp. 1284–1292.

, [*[Rigorous modal analysis of plasmonic nanoresonators]{}*]{}, Phys. Rev. B, 97 (2018), 205422.

, [*[A projection method for nonlinear eigenvalue problems using contour integrals]{}*]{}, JSIAM Lett., 5 (2013), pp. 41–44.

, [*[Riesz-projection-based theory of light-matter interaction in dispersive nanoresonators]{}*]{}, Phys. Rev. A, 98 (2018), 043806.

, [*[Resonances in physics and geometry]{}*]{}, Notices Amer. Math. Soc., 46 (1999), pp. 319–328.

[^1]: [Zuse Institute Berlin, Takustra[ß]{}e 7, 14195 Berlin, Germany]{}.

[^2]: [JCMwave GmbH, Bolivarallee 22, 14050 Berlin, Germany]{}.
