Evaluation of resonances: adaptivity and AAA rational approximation
of randomly scalarized boundary integral resolvents

We are concerned with the problem of evaluation of resonances supported by open and closed cavities and other scattering structures, which are obtained as solution pairs $(u,k)$ of the eigenvalue problem

with eigenfunction $u$ and eigenvalue $-k^{2}$, posed on an interior or exterior domain $\Omega$, with homogeneous boundary conditions of e.g. Dirichlet, Neumann, or other types. Once approximated by discretized versions of the problem’s boundary integral operators (which is done in this paper on the basis of the open- and closed-curve integral equation algorithms of [9, 23, 14], see also [10]), the resonance-search problem is reduced to the solution of a related Nonlinear Eigenvalue Problem (NEP) for a certain function $F:U\to\mathbb{C}^{d\times d}$, wherein, denoting $F(k)=F_{k}$, a complex value $k$ is sought for which

This contribution focuses on scattering problems, for which the function $F$ in (2) provides a discrete approximation of the associated boundary-integral operator. However the method is general and, indeed, other NEPs unrelated to boundary integral operators are considered in this paper as well.

The proposed approach seeks the desired resonant values $k$, for which (2) holds, as poles of the randomly scalarized resolvent

The desired resolvent poles are obtained as poles of rational approximants of (3); both the rational approximants and their poles are produced numerically by means of the recently introduced AAA rational-approximation algorithm [28]. The proposed eigensolver additionally incorporates an adaptive search strategy and a secant-method termination stage. The adaptive search strategy relies on use of a rectangular subdivision of the resonance search domain which is locally refined to ensure that all resonances in the domain are captured. The adaptivity approach works equally well in the case in which the search domain is a one-dimensional set, such as, e.g., an interval of the real line, in which case the rectangles used degenerate into subintervals of the search domain. The secant-method termination stage, in turn, is an important element in the proposed algorithm, which enables (i) Exponentially fast convergence to near machine precision accuracy starting from AAA-based results of lower accuracy; (ii) Reliable error estimation; and, (iii) A capability to reliably screen the spurious eigenvalues that can (rarely) be produced by the AAA algorithm. In all, the overall proposed approach is simple, easy to implement and rapidly convergent, and it requires limited computation besides the embarrassingly parallelizable evaluation of the scalarized resolvent at various wavenumbers $k$. A variety of numerical results presented in this paper demonstrate the character of the proposed approach: the method yields highly accurate approximations of scattering resonances and solutions of other NEPs, even in cases involving high frequencies.

A significant literature has developed in recent years in connection with the solution of NEPs. As discussed in the survey article [19], solution methods include root finding methods, contour integration methods, and methods based on linearization of $F_{k}$, all of which have been applied to the computation of resonances [2, 3, 16, 27, 31, 34]. In turn, a set of methods for the NEP that, like the present paper, rely on use of AAA rational approximation, have recently been developed [22, 17, 18], and specifically, the contributions [17, 18] apply the AAA algorithm to the scalarized resolvent (3). But in these contributions the AAA algorithm is applied in a manner different to the one we use: in those contexts a rational approximation is employed to produce a linearization of $F_{k}$ whose eigenvalues approximate the desired eigenvalues, whereas the present paper directly uses the poles of the rational approximant of $S$ as approximants of the desired eigenvalues. Closer to our work is the AAA-based algorithm introduced in [7], which considers the transmission of plane waves by a periodic dieletric system. In that approach, numerical solutions of the transmission problem are obtained by means of the finite element method and from such solutions the rational approximant of the coefficient of transmission across the structure is produced—whose poles near the real axis are indicative of the resonant character of the structure, in that they incorporate information concerning some of the structure’s complex eigenvalues.

In view of the strengths of the various algorithmic components it incorporates, the proposed algorithm is flexible, accurate and efficient. The algorithm’s use of rational approximations allows for utilization of scalarized-resolvent data points on arbitrary sets of frequencies, enabling, in particular, the use of arbitrarily distributed data points on the boundaries of rectangles—which greatly facilitates the design of the adaptive strategy proposed in this paper. Further, use of simple intervals of the real line suffices for evaluation of real eigenvalues. The secant-method termination stage, finally, provides significant benefits concerning accuracy and reliability, as mentioned in points (i)-(iii) above.

This paper is organized as follows. Sections 2 and 3 review the integral-equation and associated numerical schemes used in this paper to represent solutions of the Helmholtz equation (1) for open and closed two-dimensional domains. Section 4 provides a brief description of the AAA algorithm. The overall proposed approach for the solution of NEPs is then described in Section 5. A variety of numerical results presented in Section 6 include NEPs unrelated to Laplace eigenvalues, comparisons with well-known methods based on complex contour integration, illustrations concerning low- and high-frequency Laplace eigenvalue problems for open and closed cavities, as well as the asymptotics that result as open cavities approach closed cavities—in all, demonstrating the accuracy provided by the method even for low- and high-frequency states alike.

We consider eigenvalue problems of the form (1), posed in open two dimensional spatial domains $\Omega$ with smooth boundaries $\Gamma$, and with homogeneous boundary conditions on $\Gamma$. Three types of spatial domains are considered in this paper, namely, domains $\Omega$ equal to: (a) The complement $\Gamma^{c}$ of an open arc $\Gamma$ in $\mathbb{R}^{2}$; (b) The region interior to a closed curve $\Gamma$ in $\mathbb{R}^{2}$; and, (c) The region exterior to a closed curve $\Gamma$ in $\mathbb{R}^{2}$. For definiteness this paper mostly concerns eigenvalue problems under homogeneous Dirichlet boundary conditions

for each of these domain-types. Homogeneous Neumann and Zaremba boundary conditions can be handled similarly [23, 2, 1]; in fact, results of Neumann problems produced by the proposed algorithm are briefly mentioned in Section 6.4.

Our treatment of the problem (1)–(4) is based on use of the two dimensional Helmholtz Green function $G_{k}(x,y)\coloneq\frac{i}{4}H^{1}_{0}(k|x-y|)$ (where $H^{1}_{0}$ denotes the Hankel function of the first kind of order $0$) and the associated single-layer potential representation

of the eigenfunction $u$ in terms of a surface density $\psi$. In view of the well known [13, 25] continuity of the single layer potential $u=u(x)$ as a function of $x\in\mathbb{R}^{2}$, up to and including $\Gamma$, we consider the boundary integral operator $F_{k}:H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$ (resp.  $F_{k}:\widetilde{H}^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$) given by the expression

on a closed (resp. open) smooth curve $\Gamma$. See [25] and [23] and references therein for detailed definitions of the Sobolev spaces $H^{\pm 1/2}$ and $\widetilde{H}^{1/2}$ relevant to the closed and open-arc single layer operators, respectively.

As suggested above, the $k$-dependent operator (6) can be used to tackle the interior, exterior and open-arc eigenvalue problems under consideration. Indeed, for any smooth open arc or closed curve $\Gamma$ and for any given density function $\psi$ defined on $\Gamma$, we have [25] (i) The function $u$ given by the representation formula (5) is a solution of equation (1) for all $x\in\Gamma^{c}$; and, (ii) The function $u\neq 0$ satisfies the Dirichlet boundary condition (4) if and only if $\psi\neq 0$ is a solution of the equation $F_{k}[\psi]=0$. It can accordingly be shown that the resolvent operator $(F_{k})^{-1}$ is an analytic function of $k$ for $\Im k>0$, and that that given a complex number $k=\mu$ with $\Im\mu\leq 0$, the number $-\mu^{2}$ is an eigenvalue of the problem (1)–(4) if and only if the value $k=\mu$ is a pole of the resolvent operator $(F_{k})^{-1}$ as a function of $k$. This fact is established in [32, Prop. 7.10] for closed-curve interior problems (for which $\mu$ in the lower half-plane $\Im\mu\leq 0$ must actually be real) and for closed-curve exterior problems (for which $\mu$ must satisfy $\Im\mu<0$). As indicated in what follows, the corresponding result for open-arc problems can be established along lines similar to those in [32, Prop. 7.10].

A critical element in the extension of these results to the open-arc case is the injectivity of the mapping $\psi\to u$, which, according to equation (5), maps functions defined on $\Gamma$ to functions defined in $\mathbb{R}^{2}$. The corresponding closed-curve injectivity result for the exterior problem is established in [32] by showing that if $\psi$ satisfies the equation $F_{k}(\psi)=0$, then the associated function $u$ given by (5) is a Laplace eigenfunction in the interior of $\Gamma$ and that, therefore, the corresponding eigenvalue $-\mu^{2}$ must be real, and by subsequently making use of the jump relations for the single- and double-layer potentials across $\Gamma$. The corresponding closed-curve result for the interior problem is established similarly. For the open-arc problem we have no equivalent of the interior region, but the injectivity result can be established nevertheless, simply by using the jump relation for the normal derivative of the single-layer potential. The equivalence between Laplace eigenvalues $-\mu^{2}$ with $\mu$ in the lower half-plane and poles $k=\mu$ of the resolvent $(F_{k})^{-1}$ then follows for the open arc case in a manner analogous to that put forth in [32, Sec. 7] by relying on the second-kind formulation for open problems introduced in [9, 23].

In sum, noting that, without loss of generality, the search for values of $k$ satisfying the eigenvalue problem (1)–(4) may be restricted to the lower half-plane $\Im k\leq 0$, the eigenvalues may be sought as real poles $k=\mu$ of the resolvent operator $(F_{k})^{-1}$ for the eigenvalue problem in the interior of a closed curve $\Gamma$, and as complex poles $k=\mu$ of the same operator, with $\Im\mu<0$, for either the eigenvalue problem in the exterior of a closed curve $\Gamma$ or for the complement of an open arc $\Gamma$. The associated eigenfunctions $u$ then result via equation (5) with density (or, for multiple eigenvalues, densities) $\psi$ in the nullspace of the operator $F_{k}$. In other words, the integral equation setting just described reduces the eigenvalue problem (1)–(4) for interior and exterior closed-curve and open arc problems to an NEP for the single-layer operator (6) with values of $k$ restricted to the lower half-plane.

The numerical implementations utilized in this paper for the Green function-based integral operator (6) are based on the discretization methods presented in [14, Sec. 3.6] (resp. [9, Secs. 3.2, 5.1]) for closed (resp. open) curves in the plane. For simplicity, our computational examples restrict attention to smooth curves $\Gamma$ and boundary conditions of Dirichlet type, although related methods are available [1, 9, 14] that enable corresponding treatments for non-smooth boundaries [1, 14] as well as Neumann, Robin and Zaremba boundary conditions [1, 9]. As indicated in [9], in particular, the numerical implementation of the integral operator (6) for open arcs $\Gamma$ requires consideration of the edge singularities that are incurred by the solutions $\psi$ of problems of the form $F_{k}[\psi]=f$ even for functions $f$ defined on $\Gamma$ which do not contain such singularities.

For both open- and closed-curve problems the methods [14, 9] discretize the single layer operator $F_{k}$ on the basis of Nyström-type methodologies—utilizing a sequence of points along the curve $\Gamma$ for both integration and enforcement of the equation. For closed curves $\Gamma$ the discretization is produced as the image under the curve parametrization of a uniform grid in the parameter interval $[0,2\pi]$; in the case of open curves the discretization results as the image of a Chebyshev mesh in the parameter interval $[-1,1]$. In both cases the unknown functions $\psi$ are expressed in terms of Fourier-based expansions in the parameter intervals, which are then integrated termwise by reducing each integrated term to evaluation of explicit integrals. The edge singularities of the function $\psi$ in the open-arc case are tackled by explicitly factoring out the singular term: using a smooth parametrization $\mathbf{r}=\mathbf{r}(t)$ of the open curve $\Gamma$ ($-1\leq t\leq 1$) and writing $\psi(\mathbf{r}(t))=\phi(\mathbf{r}(t))/\sqrt{1-t^{2}}$, it follows that $\phi$ is a smooth function. Upon introduction of a cosine change of integration variables two desirable effects occur, namely, the function $\phi$ is converted into a periodic and even function which may be expanded in a cosine series with high accuracy, and, at the same time, the square root term in the denominator is cancelled by the Jacobian of the change of variables. The method is completed by exploiting explicit expressions for the integral of products of a logarithmic kernel and the cosine Fourier basis functions. As illustrated in [9] and other contributions mentioned above, these methodologies can produce scattering solutions with accuracies near machine precision on the basis of relatively coarse discretizations, even for configurations involving high spatial frequencies.

In particular, these procedures produce highly accurate numerical approximations of the integral operator $F_{k}$ in (5)—which we exploit in the context of this paper to produce accurate numerical evaluations of eigenvalues and eigenfunctions. As indicated in Section 5, the poles of (a randomly scalarized version of) this integral operator, which, per the discussion in Section 2, correspond to Laplace eigenvalues in the various cases considered, are then obtained as poles of associated AAA rational approximants. The corresponding eigenfunctions are obtained via consideration of a Gaussian elimination-based de-singularization procedure described in Section 5. For added accuracy, the methods in that section propose two alternatives, namely, the use of either iterated AAA rational approximants on one hand, or application of the secant method after an initial evaluation of poles via the AAA approach, on the other hand. In practice we have observed that, without exceptions, accuracies near machine precision are obtained for both eigenvalues and eigenfunctions on the basis of the overall proposed methodology; a few related illustrations are presented in Section 6.

The AAA algorithm is a greedy procedure for the construction of a rational approximant to a given complex-valued function $f$ on the basis of its values on an $N$-point set $Z_{N}\subset\mathbb{C}$ (full details can be found in [28]). Given the set $\{(z,f(z))\ |\ z\in Z_{N}\}$, the algorithm proceeds by selecting a sequence of points $z_{j}\in Z_{N}$, starting with some point $z_{1}$, which in principle can selected arbitrarily, but which the Matlab implementation [15] takes as a $z\in Z_{N}$ for which the function value $f(z)$ is farthest from the mean of the set $\{f(z)\ |\ z\in Z_{N}\}$. The remaining points are then selected inductively. Once points $z_{j}\in Z_{N}$ ($1\leq j\leq m$) have been chosen, with corresponding function values $f_{j}\equiv f(z_{j})$, for a suitably chosen vector $w^{m}=(w^{m}_{1},\dots,w^{m}_{m})$ of complex weights $w^{m}_{j}$ satisfying $\sum|w_{j}|^{2}=1$, the procedure to obtain $z_{m+1}$ starts by constructing the barycentric-form rational function

where, as suggested by the notation used, $n^{m}(z)$ and $d^{m}(z)$ denote the numerator and denominator in the right-hand expression in (7). The vector $w^{m}$ is selected as follows: calling $A_{m}=\{(\widetilde{w}_{1},\dots,\widetilde{w}_{m})\in\mathbb{C}^{m}\ |\sum|\widetilde{w}_{j}|^{2}=1\}$, $w^{m}$ is defined as a minimizer of the least-squares problem

where $Z^{m}_{N}=Z_{N}\setminus\{z_{j}\ |\ j=1,\dots,m\}$, $n_{\widetilde{w}}(z)=\sum_{j=1}^{m}\frac{\widetilde{w}_{j}f_{j}}{z-z_{j}}$ and $d_{\widetilde{w}}(z)=\sum_{j=1}^{m}\frac{\widetilde{w}_{j}}{z-z_{j}}$. Once the minimizer $w^{m}$ has been computed, $z_{m+1}$ is defined as a point $z\in Z^{m}_{N}$ for which $|f(z)-n^{m}(z)/d^{m}(z)|$ is maximum relative to $\max\{f(z)\ |\ z\in Z_{N}\}$. The algorithm terminates when this maximum is less than or equal to a specified tolerance, and the last rational function (7) obtained as part of the $z_{j}$ selection process provides the desired rational approximant.

All of the numerical illustrations presented in this paper utilize the AAA implementation included with Chebfun [15].

The discussion in Sections 2 and 3 reduces the eigenvalue problem (1)–(4) to NEPs for the single-layer operator (6) (and the corresponding discrete approximate operator introduced in Section 3), with values of $k$ restricted to the lower half-plane. This section presents numerical algorithms for the solution of this NEP and, indeed, of general NEPs for which the resolvent operator $(F_{k})^{-1}$ is a meromorphic function of $k$.

As indicated in Section 1, the proposed NEP algorithm obtains the eigenvalues $k$ as the poles of a AAA rational approximant associated with the randomly scalarized resolvent $S(k)$ in equation (3). While, in principle, values of the scalarized resolvent $S(k)$ on any given subset of the complex plane may be used to obtain eigenvalues, throughout this paper we restrict attention to algorithms based on use of values $S(k)$ for $k$ on a given curve $\mathcal{C}$ in the complex plane. Both open and closed curves $\mathcal{C}$ may be used, such as e.g. the closed curves equal to the boundaries of either a rectangle or a circle, or the open curve equal to an interval $[a,b]$ contained in the real axis. Proceeding on the basis of such data, Algorithm 1 below evaluates either eigenvalues contained on the union $\widetilde{\mathcal{C}}$ of $\mathcal{C}$ and its interior, if $\mathcal{C}$ is a closed curve, or eigenvalues along the curve $\mathcal{C}$, if $\mathcal{C}$ is an open curve. Algorithm 2, in turn, seeks to find all such eigenvalues, namely, all eigenvalues contained within $\widetilde{\mathcal{C}}$ if $\mathcal{C}$ is a closed curve, and all eigenvalues along $\mathcal{C}$ if $\mathcal{C}$ is an open curve.

In detail, the pseudo-code Algorithm 1 proceeds by first evaluating the scalarized resolvent $S=S(k)$ at a suitably selected set $Z_{N}=\{k_{1},\dots,k_{N}\}$ of points along $\mathcal{C}$ to obtain the set $G_{N}=\{(k_{j},s_{j})\ |\ s_{j}=S(k_{j}),\ j=1,\dots,N\}$. Using this set of pairs Algorithm 1 then obtains a rational approximant $r=r(k)$ by means of the AAA algorithm, and, for a closed curve $\mathcal{C}$, the poles of $r(k)$ contained in $\widetilde{\mathcal{C}}$ (resp. for an open curve $\mathcal{C}$, the poles on or sufficiently close to $\mathcal{C}$) are returned as approximations of eigenvalues of $F_{k}$.

Typically Algorithm 1 finds all of the eigenvalues within $\mathcal{C}$ or $\widetilde{\mathcal{C}}$, as relevant, provided the number of such eigenvalues is not too large. This observation suggests the development of an adaptive version of Algorithm 1, which, as detailed in the pseudo-code titled Algorithm 2, is specialized to searches for eigenvalues within either (i) A rectangular region $\widetilde{\mathcal{C}}$, or, (ii) A real interval $\mathcal{C}=[a,b]$. The algorithm proceeds by finding all eigenvalues within $\mathcal{C}$ or $\widetilde{\mathcal{C}}$, and then reapplying Algorithm 1 upon recursive dyadic subdivision along each dimension. The recursion ends when the number of poles obtained does not change upon subdivision. Clearly, the fact that Algorithm 1 is applicable to arbitrary curves is highly beneficial in this context, as the use of rectangular regions lends itself for the adaptive subdivision strategy used in Algorithm 2.

It is important to note that the AAA algorithm may fail to accurately produce the eigenvalues within $\mathcal{C}$ or $\widetilde{\mathcal{C}}$ if the set $G_{N}$ does not adequately represent $S(k)$ along the curve $\mathcal{C}$. An inadequate representation generally manifests itself through the appearance of false poles (see [28]) near $\mathcal{C}$, which in practice may be detected via a convergence analysis as points are added to the set $G_{N}$ and $N$ grows accordingly. Without exception, in practice we have found that, provided the curve $\mathcal{C}$ encloses a sufficiently small number of eigenvalues, once convergence to a fixed set of poles has occurred, the poles found within $\mathcal{C}$ or $\widetilde{\mathcal{C}}$ correspond in a one-to-one fashion to all the eigenvalues in that region. Starting from such poles, eigenvalues with near machine precision accuracy can be obtained by means of one of two possible methods, namely (i) Use of a subsequent “localized” AAA approximation applied to a few points around each eigenvalue obtained in the initial application of AAA; or, (ii) Use of the secant method applied to $1/S(k)=1/(u^{*}F_{k}^{-1}v)$ starting at each one of the eigenvalues obtained in the initial application of AAA. In regard to method (i), we have found that use of localized AAA approximations over a circle of radius $\rho_{\mathrm{conv}}=10^{-5}$ works well in many cases, but generally tuning of the parameter $\rho_{\mathrm{conv}}$ is necessary for convergence and, indeed, it is difficult in practice to determine whether convergence to a given tolerance has occurred, except in test cases for which the eigenvalues are known beforehand. Fortunately, method (ii) does not present this difficulty and, as it happens, it provides an additional significant advantage. Indeed, in practice we have found that, without exception, use of method (ii) results in convergence to an eigenvalue, and, clearly, the convergence history provides an error estimate for the eigenvalue found—and we thus recommend this approach as a completion procedure for each eigenvalue found. The numerical examples in Section 6 illustrate the performance of Algorithm 2 augmented by means of each one of these localized accuracy-improvement procedures.

Eigenvectors corresponding to each eigenvalue can finally be obtained via evaluation, based on Gaussian elimination, of the nullspace of the matrix $F_{k}$ at each eigenvalues $k$ obtained. In detail, using Gaussian elimination with pivoting leads to an LU decomposition of the form $F_{k}=LU$. Using this decomposition the nullspace of $F_{k}$ can be computed by a simple two-step procedure consisting of (i) Selection of a set of canonical-basis vectors that are mapped to zero by the rows in the matrix $U$ that are associated with the nonzero pivots; and, (ii) For each zero pivot in the matrix $U$, construction and solution of the reduced systems that result from the elimination of the rows and columns in $U$ containing the zero pivots, and with right-hand sides equal to the negatives of each of the eliminated columns but excluding the column element in the eliminated row. Null vectors could also be computed by means of the singular value decomposition for sufficiently small problems such as the two dimensional examples considered in this paper, but the singular value decomposition method would be problematic for larger problems such as those arising from scattering in three dimensions, whereas, in view of [30], GMRES-based iterative methods related to the LU-based approach mentioned above can be envisioned. An alternative considered in [7, Sec. 2.2] produces eigenvectors on the basis of the rational approximation itself, albeit, according to our experiments, at the expense of some loss of accuracy. In this contribution the aforementioned Gaussian elimination procedure is used, as it requires no additional approximation after application of the secant based refinement method, and as it results in accuracies comparable to those enjoyed by the eigenvalues themselves.

A few numerical illustrations of the proposed methodology are presented in what follows, including a demonstration of the performance of Algorithm 2 on NEPs unrelated to Laplace eigenvalues (Section 6.1); a set of examples concerning low-frequency Laplace eigenvalues (Section 6.2) and a comparison to complex contour integration methods (Section 6.3); an exploration of the rate at which scattering poles associated with open arcs converge to the corresponding interior eigenvalues as the opening closes (Section 6.4); and finally, a set of examples concerning high-frequency Laplace eigenvalues and eigenfunctions (Section 6.5).

This paper introduced a novel numerical algorithm for the evaluation of real and complex eigenvalues and eigenfunctions of general NEPs, including NEPs associated with Laplace eigenvalue and scattering-resonance problems for open and closed domains. Based on use of adaptive eigenvalue search methods based on AAA rational approximation combined with secant-method refinement, the algorithm produces highly accurate eigenpairs for challenging eigenproblems at both low and high frequencies. Comparisons with well-known contour-integration methods demonstrated a number of advantages of the proposed approach, including fast convergence and geometric flexibility. The latter characteristic is significant, in that it greatly facilitates use of the algorithm in a rectangular refinement-based adaptive strategy—resulting in automatic evaluation of all eigenvalues contained within a given region with near machine precision accuracy.

OB gratefully acknowledges support from NSF and AFOSR under contracts DMS-2109831 and FA9550-21-1-0373. MS acknowledges support from the National Science Foundation Graduate Research Fellowship under Grant No. 2139433. The authors gladly acknowledge helpful discussions with the following individuals: Stefan Güttel, Jonathan Keating, Yuji Nakatsukasa, Ory Schnitzer, Euan Spence, and Maciej Zworski.