Fast computation of high frequency Dirichlet eigenmodes via the spectral flow of the interior Neumann-to-Dirichlet map

• High-frequency applications may require thousands of modes up to j ∼10^6, motivating efficient methods with controlled errors.Source: §1 · p. 2
• Direct discretization requires at least N growing like k^d unknowns and can suffer pollution, while boundary methods use N = O(k^(d−1)) but require nonlinear frequency searches.Source: §1 · pp. 2–3
• Boundary-based methods can reach high mode numbers accurately, but conventional nonlinear searches cost O(N^3) per eigenfrequency and may omit eigenfrequencies.Source: §1 · p. 3
• The method linearizes the problem through the k-dependent spectrum of an interior NtD operator, whose eigenvalues reach zero at Dirichlet eigenfrequencies.Source: §1 · pp. 3–4
• The weighted formulation yields an error bound |k̂_j − k_j| ≤ C(ε^2/k + ε^3), while higher-order variants achieve empirical O(ε^5) eigenfrequency and O(ε^3) eigenfunction errors.Source: §1 · p. 4
• In two dimensions, the implementation uses potential theory, Nyström discretization, and the Cayley transform, and achieves a practical 10^3 speed-up over standard root-search methods.Source: §1 · p. 4

• The operator Θ is the inverse of the weighted Dirichlet-to-Neumann map and is self-adjoint under the inner product induced by the positive weight (x · n)^−1.Source: §2 · p. 6
• For a Dirichlet eigenmode φ_j, Θ(k) has a zero eigenvalue at k = k_j with boundary eigenfunction f = (x · n)∂_nφ_j.Source: §2 · p. 6
• The Cayley transform makes the spectral dependence analytically tractable away from −1, corresponding for Θ to analytic eigenvalue branches except for poles at Neumann eigenfrequencies.Source: §2 · p. 6
• The extended eigenfunction formulation imposes both weighted Neumann and Dirichlet conditions, with consistency ensuring uniqueness even at Neumann eigenfrequencies.Source: §2 · p. 7
• Eigenvalue branches have positive k-derivatives, including at zero crossings, so nearby Dirichlet frequencies can be inferred from the NtD spectral flow.Source: §2 · pp. 7–8
• The weight (x · n)^−1 equalizes zero-crossing speeds independently of eigenmode structure, unlike the unweighted map whose slopes vary unpredictably.Source: §2 · p. 8

3.1. Reconstructing eigenfrequencies.

• The method extrapolates each small negative NtD eigenvalue along its approximately linear spectral branch to estimate a corresponding Dirichlet eigenfrequency.A frequency window of width O(1) contains O(k^(d−1)) eigenfrequencies, so adjacent windows can recover the spectrum across the frequency axis.Source: §3.1 · p. 8
• In the [90, 100] test window, 492 approximate eigenfrequencies were computed in 13 minutes, with most absolute errors below 10^-4.The discretization used 720 boundary quadrature points, about six per wavelength, and required 1.6 seconds per computed eigenfrequency.Source: §3.1 · p. 9
• The eigenfrequency error has O(ε^2/k) and O(ε^3) contributions with constants independent of k, yielding asymptotically third-order accuracy in the window width.The transition between the two error regimes occurs at ε = O(1/k).Source: §3.1 · pp. 9–10, Figure 2 · p. 9

3.2. Reconstructing eigenfunctions from boundary data.

• Boundary-function errors are first order in the prediction distance, and estimated domain eigenfunction errors cluster near a unit-slope relation with boundary errors.The domain-error study used 327 interior grid points and had statistical error of order ±10%; rigorous control would require bounds on the k-dependence of S(k).Source: §3.2 · p. 11, Figure 3 · p. 11
• 30 times faster than the reference method for extracting boundary eigenfunctions, despite increasing runtime from 1.6 to around 2.3 seconds per mode.Complete diagonalization is required to obtain boundary eigenfunctions, while the reference method also incurs additional mode-extraction time.Source: §3.2 · p. 12

3.3. Numerical computation of spectrum of Θ.

• The planar implementation combines potential-theory layer operators with Nyström quadrature and dense O(N^3) linear algebra to approximate the Cayley spectrum and NtD eigenpairs.For analytic curves, the product quadrature is spectrally accurate, although the paper does not provide a rigorous error analysis for the noncompact operator product and inverse.Source: §3.3 · pp. 12, 14, 15
• The Cayley transform keeps the transformed operator spectrum O(1), avoiding the large dynamic range and divergence affecting direct NtD discretization near Neumann eigenfrequencies.Numerically, the direct method was limited to 8 digits or failed at a Neumann eigenfrequency, whereas the Cayley method achieved 14-digit accuracy uniformly in k.Source: §3.3 · pp. 13–14

• The accuracy result is unconditional, unlike the higher-order method whose sharper bounds require absence of spectral concentration.Source: §4 · p. 15
• The weighted NtD eigenvalue flow deviates from linearity by a controlled amount over intervals [k_p - epsilon, k_p].The deviation bound is proved using estimates on the eigenvalue branch and integration in k.Source: §4 · pp. 15, 17, 18
• The basic estimator has error |k_hat - k_p| <= C(epsilon^2/k + epsilon^3), with constants independent of k.The estimate applies to sufficiently large frequencies within a fixed window around the reference frequency.Source: §4 · p. 15
• For a fixed window width epsilon, O(k^(d-1)) eigenfrequencies are recovered, and the basic method is effectively O(epsilon^3) accurate for almost all at high k.The epsilon^2/k term becomes negligible as k grows, while the constant remains independent of k.Source: §4 · p. 15
• The proof derives the eigenvalue-flow estimate from boundary identities and tangential-gradient bounds for high-frequency Helmholtz solutions.The argument uses the tangential field W, estimates its associated terms, and integrates the resulting differential control in k.Source: §4 · pp. 16–17

5.1. Higher-order eigenfrequency approximation.

• The eigenfrequency Riccati estimator freezes boundary-function norms and solves a constant-coefficient Riccati equation after changing variables to log k.The choice of freezing point based on the mean of k* and the linear estimator gives slightly smaller errors than kz = k*.Source: §5.1 · pp. 19–20
• 103 to 105 times better eigenfrequency errors are observed for the Riccati estimator than for the linear estimator, with error O(ε^3/k^2 + ε^5) and about 10^-7 absolute error at ε = 0.1.The ε^5 term dominates for windows large enough to collect many eigenfrequencies, and its constant is independent of k.Source: Figure 5 · p. 19, §5.1 · p. 20

5.2. Higher-order reconstruction of eigenfunctions.

• The higher-order eigenfunction reconstruction differentiates eigenfunction branches of Θ(k), assumes a simple eigenspace, and uses first- and second-derivative formulas to build linear and quadratic estimators.The quadratic construction retains leading terms because explicitly applying Θ(kp)^-1 at every eigenfrequency would destroy O(N^2) complexity per mode.Source: §5.2 · pp. 21, 23, 24
• A quadratic eigenfunction estimator gives typically one extra digit over the linearized estimator and up to four extra digits over the trivial estimator.Its empirical error is O(ε^2k^-1/2 + ε^3k^1/2), reducing to O(ε^3k^1/2) in the fast regime.Source: Figure 6 · p. 21, §5.2 · p. 24
• The higher-order eigenfunction method requires ǫ < k^-1/6 to keep its fast-regime error bounded as k increases, reducing the speed-up from O(N) to O(N^5/6) in d = 2.This restriction is contrasted with ǫ = O(1), which is allowed for bounded eigenfrequency error.Source: §5.2 · p. 24

5.3. Performance in a domain with abundant degeneracies.

• For the pentafoil domain, degeneracies have no deleterious effect on NtD error levels, and the method pays no penalty for close or exact degeneracies.The reference method slows by roughly a factor of three when resolving eigenfrequency pairs, whereas NtD speed-up factors become around three times better.Source: §5.3 · pp. 25–26
• At k ≈ 1000, a matrix of size N = 9000 took 3.6 hours to diagonalize and returned 51 modes, with an estimated speed-up of 10^3 over the reference method.The reference computation was not run at this scale; the speed-up estimate uses O(N^3) scaling from previous tests.Source: Figure 8 · p. 26, §5.3 · p. 27

6.1. A spectral nonconcentration assumption.

• In d=2, absence of spectral concentration isolates a negative NtD eigenvalue branch from neighboring branches across the estimation interval.The assumption requires the target eigenvalue to be alone in a prescribed interval and therefore simple.Source: §6.1 · pp. 27–28
• Average Dirichlet eigenfrequency spacing is ∼1/k in two dimensions, so the corresponding NtD eigenvalues are typically spaced ∼1/k^2.This motivates expecting the nonconcentration assumption to hold for most sufficiently small-scale eigenvalues.Source: §6.1 · pp. 27–28

6.2. Error estimate for second-order eigenfunction reconstruction.

• Under the same assumption, the second-order eigenfunction reconstruction estimator satisfies the error bound stated in Proposition 6.4.The estimate controls reconstruction of the weighted normal derivative at the true eigenfrequency.Source: §6.2 · pp. 28–29

6.3. Error estimate for the Riccati eigenfrequency estimator.

• The Riccati eigenfrequency estimator satisfies the rigorous error estimate stated in Proposition 6.5 under the d=2 nonconcentration assumption.The proof bounds accumulated errors in the estimator terms and converts them to frequency error using dβ/dk comparable to 1/k.Source: §6.3 · p. 30
• For ε=O(1), the rigorous high-frequency bound is dominated by O(kε^6/η), although empirical errors are dominated by O(ε^5).Using a modified frozen frequency also removes the empirical O(ε^4/k) term at intermediate ε values.Source: §6.3 · pp. 30–31

7.1. Sketch of the scaling method.

• The scaling method linearizes a self-adjoint dilation-based operator using global Helmholtz basis functions and a generalized eigenvalue problem.Its basis functions may be plane waves, fundamental solutions, or Fourier-Bessel wedge solutions, depending on the domain.Source: §7 · p. 31, §7.1 · p. 31
• Ill-conditioning of global bases makes the scaling method’s mass and derivative matrices share a numerical nullspace, requiring regularization.The regularization parameter is not fully understood, and reconstructed eigenfrequency accuracy is reported empirically as O(ε^3).Source: §7.1 · p. 32

7.2. Connecting scaling and NtD methods via dilation.

• The scaling and NtD methods are connected because the dilation operator is asymptotically related to the inverse weighted NtD operator near a Dirichlet eigenfrequency.This relation predicts comparable eigenfrequency accuracy, while the second-order flows also agree to highest order for eigenfunction reconstruction.Source: §7.2 · p. 32

7.3. Advantages of NtD method over the scaling method.

• Both methods share an O(N) acceleration and similar eigenfunction error, but the NtD method permits higher-order frequency accuracy and gives 3 to 5 extra correct digits in practice.The NtD formulation also reconstructs modes without continuing basis functions outside the domain.Source: §7.2 · p. 33, §7.3 · p. 33
• The NtD method offers rigorous error estimates, robustness across choices of k*, and spectrally accurate boundary-integral discretizations, whereas scaling accuracy can degrade near Dirichlet eigenfrequencies.The comparison also notes that global-basis efficiency in scaling can be higher in some cases, so the advantage is not universal across every implementation detail.Source: §7.3 · p. 33

8. Conclusions

• The method linearizes weighted NtD spectral flow over frequency windows, enabling an O(k^(d−1)) = O(N) speed-up without root-searching.The weight (x·n)^−1 equalizes flow gradients and prevents avoided crossings; independent windows are embarrassingly parallel.Source: §8 · p. 33
• The method proves neither spurious nor missing modes, third-order absolute eigenfrequency accuracy, and conditional third-order L2 boundary-mode error.The higher-order estimator is empirically fifth-order for eigenfrequencies and third-order for modes.Source: §8 · pp. 33–34
• The algorithm’s operator-level acceleration is independent of additional solver acceleration, but the benefit of fast solvers remains uncertain for oscillatory high-frequency kernels.The authors identify testing block iterative, fast multipole, and fast direct approaches as future work.Source: §8 · p. 34
• Future work includes removing the star-shaped restriction, treating other boundary conditions and nonsmooth domains, and improving or analyzing higher-order estimators.The stated extensions include Neumann or Robin conditions, corners, boundary error bounds, and high-order extrapolation.Source: §8 · p. 34

Appendix A. Smoothness of eigenvalues and eigenprojections in k

• The boundary-value problem has a unique holomorphic solution near the positive real axis, so the Cayley transform R(k) is analytic there.The proof uses compactness and analytic Fredholm theory after establishing triviality of the relevant null space.Source: Appendix A · pp. 34–35
• Analyticity of R(k) implies analytic eigenvalues and eigenspaces for Θ(k) away from poles and zero crossings.The exceptional points include finitely many branches with poles at Neumann eigenfrequencies and crossings where NtD eigenvalues reach zero.Source: Appendix A · p. 36
• A zero eigenvalue of Θ(k*) generates analytic branches approaching zero from below, with multiplicity matching the Dirichlet eigenspace multiplicity.Conversely, any branch tending to zero identifies k* as a Dirichlet eigenfrequency.Source: Appendix A · pp. 36–37
• The limiting eigenspaces are spanned by weighted normal derivatives of Dirichlet eigenfunctions, and simple branches have eigenfunctions that are C2 up to k=k*.This regularity supports tracking eigenprojections and eigenfunction flows through the target frequency.Source: Appendix A · pp. 36–38
• The proof obtains nonzero Dirichlet limits by compactness of extended eigenfunctions and strong convergence in the interior and on the boundary.For multiple branches, orthogonal sequences yield a Dirichlet eigenspace whose dimension matches the branch multiplicity.Source: Appendix A · p. 37

Reference eigenfrequency computation

• The reference criterion is equivalent to a nontrivial nullspace of (1/2 − D*(k)), whose vectors represent boundary normal derivatives of Dirichlet eigenfunctions.This characterization follows from jump relations and uniqueness of the exterior Helmholtz Neumann problem.Source: Appendix B · p. 38
• O(N^3) effort per mode is required because each minimum-singular-value evaluation and eigenmode extraction uses dense linear algebra.The method searches near-zeros of (1/2 − D*(k)); each evaluation and singular-vector computation costs O(N^3).Source: Appendix B · pp. 38–39
• Eigenfrequencies are identified as minima of a series of V-shaped smallest-singular-value curves, refined by grid searches or parabolic fitting.Close eigenfrequencies can create W-shaped features, so the implementation adaptively refines the grid using higher singular values.Source: Appendix B · p. 38
• The implementation finds all eigenfrequencies in an interval with about 15 evaluations per simple eigenfrequency and typical errors of 10^-13 or less.For degenerate eigenfrequencies, the cost rises to roughly 50 evaluations per multiple eigenfrequency at tolerance 10^-6, with typical errors of 10^-7.Source: Appendix B · p. 39

Reference eigenmode computation

• After locating each frequency, the last right singular vector approximates the boundary normal derivative, which is normalized and used to reconstruct the eigenfunction.For a p-fold degeneracy, the last p right singular vectors are used.Source: Appendix B · p. 39

Spectral cutoff and elliptic localization

• The proof introduces a spectral cutoff to separate high-energy boundary frequencies before estimating tangential derivatives.The cutoff is supported in the elliptic region and is treated as a semiclassical pseudodifferential operator of order (0, 0).Source: Appendix C · p. 39

Conclusion of the proof

• Combining the cutoff, pseudodifferential, and boundary estimates proves the stated high-energy bound for k ≥ K.The final step uses the estimate 8(3k^2 + 4C^2)/7 ≤ 4k^2.Source: Appendix C · pp. 39, 41

Boundary-frequency estimates

• Low-frequency components are controlled because the complementary cutoff removes tangential frequencies of order at least 2k, while commutator bounds remain uniform in k.Summing over an orthonormal tangential vector-field basis converts the estimate into a weighted boundary Laplacian bound.Source: Appendix C · pp. 39–40

Semiclassical operator analysis

• In the elliptic region, S(k) and D(k) have order (−1, −1), while Q(k) has order two higher with principal symbol −1 − |η|^2.The apparent order discrepancy for D(k) is explained by the boundary-supported half-identity term in the jump formula.Source: Appendix C · pp. 40–41
• The weighted operator identities reduce the difficult high-energy term to a negative square plus lower-order contributions, allowing the dominant term to be discarded in the estimate.The crucial sign yields a decomposition involving −B*B + A2, with A2 of lower order.Source: Appendix C · pp. 40–41

Assumptions and scope

• The analysis is developed for near-Dirichlet Robin boundary conditions and relies on smoothness, spectral separation, and weighted boundary estimates.The proof extends earlier arguments because the relevant functions do not vanish on the boundary.Source: Appendix D · pp. 42–43

Supporting operator bound

• Boundary estimates for approximate eigenfunctions control the operator mapping boundary data to normal derivatives, completing the proof of the theorem through TT* norm identities.The argument uses the spectral window [(k − 1)^2, (k + 1)^2] and the equality ||TT*|| = ||T||^2.Source: Appendix D · p. 43

Resolvent estimate

• The main estimate bounds (Θ(k) − β)^−1 using the distance d(Eβ, σ*) from the relevant energy to the nearest spectrum.Under the stated assumptions, the spectral distance is at least Cη, yielding a remainder bound proportional to k^2/η.Source: Appendix D · pp. 44–45
• The proof projects away from the Θ(k) eigenfunction associated with β, solves a near-Dirichlet boundary problem, and decomposes the solution into v − v′ and v′.Negativity of β guarantees uniqueness for the auxiliary problem.Source: Appendix D · p. 44

Worst-case behavior

• The worst-case resolvent growth occurs when the input aligns with a nearby boundary normal derivative associated with a distinct but close eigenvalue.This scenario produces the factor k^2/d(Eβ, σ*).Source: Appendix D · p. 45

Typical-case behavior

• Typically, the resolvent is expected to scale like k^3/2/d(Eβ, σ*), making the Θ(k)^−1(mf) term smaller by about k^−1/2 in the quadratic estimator.This is presented as a heuristic justification for dropping that term, not as a worst-case bound.Source: Appendix D · pp. 45–46