Computation of eigenfrequency sensitivities using Riesz projections for efficient optimization of nanophotonic resonators

Resonances are omnipresent in physics and essential for the description of wave phenomena. We present an approach for computing eigenfrequency sensitivities of resonances. The theory is based on Riesz projections and the approach can be applied to compute partial derivatives of the complex eigenfrequencies of any resonance problem. Here, the method is derived for Maxwell’s equations. Its numerical realization essentially relies on direct differentiation of scattering problems. We use a numerical implementation to demonstrate the performance of the approach compared to differentiation using finite differences. The method is applied for the efficient optimization of the quality factor of a nanophotonic resonator. Resonances are ubiquitous in physics and hold important functionalities in engineering wave propagation and interference effects. This work proposes an approach for computing sensitivities, i.e., partial derivatives, of complex eigenfrequencies in any resonance problem, which here is applied to efficiently optimize nanophotonic resonators and to obtain an improved quality factor.

R esonance phenomena are ubiquitous in nanophotonics and play an important role for tailoring light-matter interactions 1,2 . They are exploited in, e.g., single-photon sources for quantum technology 3 , biosensors 4 , nanolasers 5 , or solar energy devices 6,7 . All these applications rely on the highly localized electromagnetic field energies in the vicinity of the underlying nanoresonators 8 . A central figure of merit for the description of resonance effects is the quality (Q) factor, which quantifies, in the case of low-loss systems, the relation between stored and radiated field energies of the resonances 9 . Nanoresonators with low energy dissipation, i.e., with high Q-factors, have been proposed to improve the efficiencies of nanophotonic devices 2,10 . For example, high-Q resonators can boost the brightness of quantum emitters, the sensitivity of sensors, or the emission processes in plasmonic lasers 11 . Designing devices with numerical optimization is a time and cost-effective approach. The resonances are numerically computed by solving the source-free Maxwell's equations equipped with open boundary conditions 12 . This yields non-Hermitian eigenproblems and the solutions are eigenmodes with complex-valued eigenfrequencies. In this context, the Q-factor is defined as the scaled ratio of the real and imaginary parts of the eigenfrequency.
Nanoresonators with high Q-factors have been theoretically presented, but fabrication of these resonators is a limiting task 11 . The sensitivity analysis of eigenfrequencies can show a way to reduce the sensitivities of the Q-factors. This can support the nanofabrication processes. Furthermore, the sensitivity analysis of eigenfrequencies is essential for numerical simulation. For example, the numerical accuracies of the calculated eigenfrequencies are strongly influenced by the sensitivities of the eigenfrequencies when the systems are subject to small perturbations 13,14 . In particular, for high-Q resonators, the accuracy requirements are demanding since the real and imaginary parts of the eigenfrequencies differ by several orders of magnitude. Sensitivities are also directly exploited in numerical optimization algorithms using gradients 15 , for gradient-enhanced surrogate modeling 16 , and for local sensitivity analyses 17 . The computation of eigenfrequency sensitivities is usually based on perturbation theory 18,19 , where the sensitivity of the underlying operator, the left and the right eigenmodes, and a proper normalization of the eigenmodes are required. The solution of the perturbed systems, on the other hand, is not necessary. For resonance problems, left and right eigenmodes are in general not identical, which increases the computational effort, and normalization requires additional attention. There are specialized approaches that, e.g., exploit magnetic fields for extracting the left eigenmodes 20 , introduce an adjoint system for computing sensitivities 21 , or that rely on internal and external electric fields at the boundaries of the nanoresonators 22 . It is also possible to completely omit the use of eigenmodes for sensitivity analysis 23 . A further approach is the straightforward application of finite differences. However, this also includes the solution of the perturbed resonance problems, which increases the computational effort.
In this work, we present an approach for computing eigenfrequency sensitivities that completely avoids solving resonance problems. The approach is based on Riesz projections given by contour integrals in the complex frequency plane. The contour integrals are numerically accessed by solving Maxwell's equations with a source term enabling an efficient numerical realization using direct differentiation. The numerical experiments show a significant reduction in computational effort compared to applying finite differences. A Bayesian optimization algorithm with the incorporation of eigenfrequency sensitivities is used to optimize a resonator hosting a resonance with a high Q-factor.

Results
Theoretical background and numerical realization. We start with an introduction of the theoretical background on resonance phenomena occurring in nanophotonics. Based on this, Riesz projections for computing eigenfrequency sensitivities and an efficient approach for its numerical realization are presented.
Resonances in nanophotonics. In nanophotonics, in the steadystate regime, light-matter interactions can be described by the time-harmonic Maxwell's equations in second-order form, where Eðr; ω 0 Þ 2 C 3 is the electric field, r 2 R 3 is the position, ω 0 2 R is the angular frequency, and JðrÞ 2 C 3 is the electric current density corresponding to a light source. In the optical regime, the permeability tensor μ(r, ω 0 ) typically equals the vacuum permeability μ 0 . The permittivity tensor ϵ(r, ω 0 ) = ϵ r (r, ω 0 )ϵ 0 , where ϵ r (r, ω 0 ) is the relative permittivity and ϵ 0 the vacuum permittivity, describes the spatial distribution of material and the material dispersion. Solutions to Eq. (1) are called scattering solutions as light from a source is scattered by a material system. Resonances are solutions to Eq. (1) without a source term, i.e., J(r) = 0, and with transparent boundary conditions. The boundary conditions lead to non-Hermitian eigenproblems, and, if material dispersion is also present, the eigenproblems become nonlinear. The electric field distribution of an eigenmode is denoted byẼðrÞ 2 C 3 and the corresponding complex-valued eigenfrequency byω 2 C. The Q-factor of a resonance is defined by Q ¼ ReðωÞ À2ImðωÞ and describes its spectral confinement. In the limiting case of vanishing losses, this definition agrees with the energy definition, according to which the Q-factor quantifies the relation between stored and dissipated electromagnetic field energy of a resonance 9 .
In the following, a nanophotonic resonator supporting a resonance with a high Q-factor is investigated. We compute the eigenfrequency sensitivities with respect to various parameters to optimize the Q-factor of the nanoresonator. Figure 1 sketches the applied framework for an exemplary problem, a one-dimensional resonator defined by layers with different permittivities. Changes δp of the parameter p lead to changes in the eigenmodeẼ and in the corresponding eigenfrequencyω, which describes the sensitivity ofẼ andω with respect to the parameter p. To compute the eigenfrequency sensitivity, we introduce a contourintegral-based approach using Riesz projections, where physical observables are extracted from scattering problems. Solving the scattering problems, which are linear systems, can be regarded as a blackbox 24,25 .
Riesz projections for eigenfrequency sensitivities. To derive a Rieszprojection-based approach for computing eigenfrequency sensitivities, which are the partial derivatives of the eigenfrequency, we consider the electric field Eðr; ω 0 2 RÞ as a solution of Eq. (1) and Eðr; ω 2 CÞ as an analytical continuation of E(r, ω 0 ) into the complex frequency plane. The field E(r, ω) is a meromorphic function with resonance poles at the eigenfrequencies. To simplify the notation, we omit the spatial and frequency dependency of the electric field and write E when we mean E(r, ω).
Let LðEÞ be a physical observable, where L : C 3 ! C is a linear functional, andC be a contour enclosing the poleω of the order m and no other poles. Then, the Laurent expansion of LðEÞ aboutω is given by The coefficient a À1 ðωÞ is the so-called residue of LðEÞ atω. Using Eq. (2) with the assumption thatω has the order m = 1 and applying Cauchy's integral formula yield Ĩ where, due to the closed integral in the complex plane, the regular terms in the expansion vanish. With this, the eigenfrequencyω is given byω The contour integrals in this equation are essentially Riesz projections for LðEÞ andC 24 . Partial differentiation with respect to a parameter p directly gives the desired expression for the eigenfrequency sensitivity, For the interchangeability of integral and derivative, E and ∂E/ ∂p are assumed to be continuously differentiable with respect to the frequency ω and the parameter p. The eigenmodeẼ and its sensitivity ∂Ẽ=∂p can be represented by the contour integrals respectively, which are Riesz projections applied to Maxwell's equations given by Eq. (1). This approach can be generalized for multiple eigenfrequencies inside a contour as well as for higher order poles; cf. Binkowski et al. 24 . Note that Riesz projections can also be used to compute modal expansions of physical observables, where scattering solutions are expanded into weighted sums of eigenmodes 26 .
Numerical realization and direct differentiation. For the numerical realization of the presented approach, the finite element method (FEM) is applied. Scattering problems are solved by applying the solver JCMSUITE 27 . The FEM discretization of Eq.
(1) leads to the linear system of equations AE = f, where A 2 C n n is the system matrix, E 2 C n is the scattered electric field in a finite-dimensional FEM basis, and f 2 C n contains the source term. The solver employs adaptive meshing and higher order polynomial ansatz functions. In all subsequent simulations, it is ensured that sufficient accuracies are achieved with respect to the FEM discretization parameters. Note that also other methods can be used for numerical discretization. In the field of nanophotonics, common approaches are, e.g., the finite-difference timedomain method, the Fourier modal method, or the boundary element method 12,28 .
In order to calculate eigenfrequenciesω and their sensitivities ∂ω=∂p i with respect to parameters p i , the electric fields E and their sensitivities ∂E/∂p i are computed for complex frequencies ω 2 C on the contours given in Eqs. (3) and (4). For the calculation of ∂E/∂p i , we apply an approach based on directly using the FEM system matrix 29,30 . With this direct differentiation method, the sensitivities of scattering solutions can be computed by In a first step, instead of directly computing A −1 , an LUdecomposition of A, which can be seen as the matrix variant of Gaussian elimination, is computed to efficiently solve the linear system AE = f. In the FEM context, this step is usually a computationally expensive step in solving scattering problems, so reusing an LU-decomposition can significantly reduce computational costs. In a second step, the partial derivatives of the system matrix, ∂A/∂p i , and of the source term, ∂f/∂p i , are obtained quasi analytically, i.e., with negligible computational effort. Then, A = LU, E, ∂A/∂p i , and ∂f/∂p i are used to compute ∂E/∂p i in Eq. (5). The LUdecomposition can be used to obtain both E and ∂E/∂p i .
For the calculation of the contour integrals, a numerical integration with a circular integration contour and a trapezoidal rule is used, which leads to an exponential convergence behavior with respect to the integration points 31 . At each integration point, we calculate E and ∂E/∂p i by solving Eq. (1) with oblique incident plane waves as source terms. The linear functional LðEÞ corresponds to a spatial point evaluation of one component of the electric field, which can be understood as physical observable. Note that, with Eqs. (3) and (4), an eigenfrequencyω and its Fig. 1 Schematic representation of computing eigenfrequency sensitivities of a resonator using contour integration. The system is defined by layers with different permittivities ϵ 1 and ϵ 2 and is described by the one-dimensional Helmholtz equation ÀΔẼ Àω 2 ϵẼ ¼ 0. A solution to the resonance problem is given by the eigenmodeẼ and the corresponding complex-valued eigenfrequencyω 2 C. The real part of the electric field of the eigenmode is sketched with the solid black curve. A perturbation δp of the middle layer width p leads to a perturbed electric field, represented by the dashed red curve, and to a perturbation δω of the eigenfrequency. Computing contour integrals by solving linear systems AE = f and ∂=∂p AE ¼ f ½ in the complex frequency plane yields the eigenfrequency sensitivity ∂ω=∂p. Solving the linear systems is considered as a blackbox.
sensitivity ∂ω=∂p i can be calculated without solving resonance problems ∇ μ À1 ∇ Ẽ Àω 2 ϵẼ ¼ 0 directly. Instead, scattering problems, where Eq. (5) can be exploited, are solved. We call the described approach, which combines Riesz projections and direct differentiation (DD), the Riesz projection DD method. Equation (4) and its numerical implementation exploiting Eq. (5) are the main results of this work and represent the difference from previous works on Riesz projections; cf. Zschiedrich et al. 26 .
Note that the Riesz projection DD method is not limited to the field of nanophotonics, but can be applied to other eigenproblems as well. Maxwell's equations can be replaced by another partial differential equation, and then instead of the analytical continuation of the electric field E, the analytical continuation of another quantity is evaluated for the contour integration.

Application
Eigenfrequency sensitivities of a nanophotonic resonator. We investigate an example from the literature, a dielectric nanoresonator of cylindrical shape placed on a three-layer substrate, where constructive and destructive eigenmode interference has been used to engineer a bound state in the continuum (BIC) 32 . The nanoresonator has been designed taking into account various parameters to suppress radiation losses: The radius, the layer thicknesses, and the layer materials have been chosen to obtain a high-Q resonance. The nanoresonator is made of the high-index material aluminum gallium arsenide (AlGaAs) with 20% aluminum. A silicon dioxide (SiO 2 ) spacer is placed between the nanoresonator and a film of indium tin oxide (ITO) on a SiO 2 substrate. A sketch of the designed system is shown in Fig. 2a. For this specific configuration, a high-Q resonance with a Q-factor of Q = 188 ± 5 has been experimentally observed, and numerical simulations have resulted in Q = 197, where the real part of the resonance wavelength is in the telecommunication wavelength regime, close to 1600 nm. The nanophotonic resonator has been exploited as a nanoantenna for nonlinear nanophotonics 32 .
In the following simulations, we consider the constant relative permittivities ϵ r = 10.81 and ϵ r = 2.084 for AlGaAs and for SiO 2 , respectively, which are extracted from experimental data 32,33 . For the ITO layer, the Drude model ϵ r ðω 0 Þ ¼ ϵ inf À ω 2 p =ðω 2 0 þ iω 0 γÞ is chosen, where ϵ inf ¼ 3:8813, ω p = 3.0305 × 10 15 s −1 , and γ = 1.2781 × 10 14 s −1 . This Drude model is obtained by a rational fit 34 to experimental data 32 and describes the material dispersion of the system. We further exploit the rotational symmetry of the geometry. On the one hand, this reduces the computational effort and, on the other hand, the eigenmodes can be easily distinguished by their azimuthal quantum numbers m, which correspond to the number of oscillations in the radial and axial directions. When the light sources used for computing Riesz projections are not rotationally symmetric, such as oblique incident plane waves, the source fields can be expanded into Fourier modes in the angular direction. Considering Fourier modes with certain quantum numbers, only the eigenmodes, eigenfrequencies, and corresponding sensitivities associated with these quantum numbers are accessed.
We start with computing a Riesz projection to obtain the eigenfrequencyω of the high-Q resonance. Figure 2b shows the complex frequency plane with the calculated eigenfrequency, ω ¼ ð1:17309 À 0:00296iÞ 10 15 s À1 , and the corresponding circular integration contourC for the computation of the Riesz projection. The center and the radius of the contour are selected based on a-priori knowledge from Koshelev et al. 32 . Alternatively, without a-priori knowledge, a larger integration contour can be used 25 . The simulations are performed using eight integration points on the contourC, where a sufficient accuracy with respect to the integration points is ensured. The computations are based on a FEM mesh consisting of 306 triangles. To compare the size of the contour with the distances between the eigenfrequencies within the spectrum of the nanoresonator, the two eigenfrequencies which are closest toω are also shown. We obtain a Q-factor of Q = 198 for the high-Q resonance, which is in good agreement with the experimental and numerical results from Koshelev et al. 32 . The corresponding electric field intensity jẼj 2 is shown in Fig. 2c. The eigenmodeẼ has the quantum number m = 0 and is strongly localized in the vicinity of the nanoresonator.
Next, the eigenfrequency sensitivities ∂ω=∂p i with respect to the parameters p 1 , p 2 , …, p 5 sketched in Fig. 2a are computed. In order to validate the approach, a convergence study for the polynomial degree d of the FEM ansatz functions is performed. Figure 2d, e shows the relative errors for the real and imaginary parts, respectively. Exponential convergence can be observed for all sensitivities with increasing d. The computed sensitivities for d = 5 are shown in Table 1. Exemplary source code for the Riesz projection DD method and simulation results are presented in Binkowski et al. 35 .
Performance benchmark. The computational effort of the numerical realization of the Riesz projection DD method is compared with the computational effort of the finite difference method. We choose the central difference scheme ∂ω=∂p i % ωðp i þ δp i Þ Àωðp i À δp i Þ À Á = 2δp i À Á for the comparison. Computing central differences is more computationally expensive than computing forward or backward differences. However, more accurate results can be achieved as the error decreases with ðδp i Þ 2 . To achieve an adequate accuracy, sufficiently small step sizes δp i are selected. For example, for the radius of the nanoresonator, we choose δp 1 = 0.1 nm. Note that, also for the finite difference method, we compute the eigenfrequencies by using the contourintegral-based formula in Eq. (3).
We increase the degrees of freedom of the system shown in Fig. 2a by deforming the cylindrical nanoresonator to an ellipsoidal nanoresonator. This breaks the rotational symmetry yielding a full three-dimensional system with new parameters, the radius of the nanoresonator in x direction and the radius in y direction. Figure 3 shows, for the three-dimensional implementation and for the rotational symmetric implementation, the normalized computational effort for different numbers of computed sensitivities. We compute the eigenfrequencyω and then we add the sensitivities, starting with ∂ω=∂p 1 , one after the other. It can be observed that the Riesz projection DD method requires less computational effort than the finite difference method, for any number of computed sensitivities, i.e., for all N ≥ 1. In the case of using finite differences, the computational effort has a slope of about 200% because for each sensitivity two additional problems with typically the same dimension as the unperturbed problem have to be solved. In the three-dimensional case, a linear regression for the computational effort gives a slope of about 4% for the Riesz projection DD method. The computational effort needed for the LU-decomposition is significant compared to the matrix assembly and to the other solution steps, so the possibility of exploiting Eq. (5) gives a great benefit for the Riesz projection DD method. For N = 5, the CPU time required to solve the linear system of equations, which includes the LU-decomposition, takes 81% of the accumulated CPU time. In the rotational symmetric case, the time for solving the linear system is negligible. However, the trend is the same for the three-dimensional and for the computationally cheaper rotational symmetric case: The advantage of using Riesz projections significantly increases with an increasing number of computed sensitivities.
Note that contour integral methods are well suited for parallelization because the scattering problems can be solved in parallel on the integration contour. However, as total CPU times are considered for Fig. 3, this is not reflected by the time measurements.
Q-factor optimization. The Riesz projection DD method is applied to further optimize the Q-factor of the high-Q resonance of the nanophotonic resonator from Koshelev et al. 32 shown in Fig. 2a. A rotational symmetric nanoresonator is considered because simulations show that an ellipsoidal shape does not lead to a significant increase of the Q-factor. We use a Bayesian optimization algorithm 36 with the incorporation of sensitivity information. This global optimization algorithm is well suited for problems with computationally expensive objective functions and benchmarks show that providing sensitivities can significantly reduce computational effort 37 . However, other optimization approaches could be used as well.
For the optimization, we choose the parameter ranges 435 nm ≤ p 1 ≤ 495 nm, 575 nm ≤ p 2 ≤ 695 nm, 150 nm ≤ p 3 ≤ 550 nm, 100 nm ≤ p 4 ≤ 500 nm, and 60 ∘ ≤ p 5 ≤ 90 ∘ . To ensure that the optimized nanoresonator can also be used as nanoantenna in the telecommunication wavelength regime, like the original system, we add the constraint that the optimized eigenfrequency must lie in the circular contour with the center ω 0 = 2πc/(1600 nm) and the radius r 0 = 4 × 10 13 s −1 . In each optimization step, the Riesz  projection DD method is used to compute the eigenfrequency with a quantum number of m = 0 lying inside the contour and to calculate the corresponding sensitivities.
A nanoresonator with a Q-factor of Q = 292 is obtained after 61 iterations of the optimizer yielding an increase of about 47.5% over the original resonator. More iterations yield only a negligible increase of the Q-factor. The optimized nanoresonator with a sketch of the electric field intensity of its high-Q resonance and the values for all underlying parameters are shown in Fig. 4. The corresponding eigenfrequency is given byω opt ¼ ð1:176897À 0:002015iÞ 10 15 s À1 . Note that, in the optimization domain, the average sensitivity of the Q-factor with respect to the ITO layer thickness p 4 is negligible.

Conclusions
An approach for computing eigenfrequency sensitivities of resonance problems was presented. The numerical realization of the Riesz projection DD method relies on computing scattering solutions and their sensitivities by solving Maxwell's equations with a source term, i.e., solving linear systems of equations. This enables direct differentiation for the efficient calculation of eigenfrequency sensitivities. Although sensitivities of resonances are computed, no eigenproblems have to be solved directly. The performance of the approach was demonstrated by a comparison with the finite difference method. The Riesz projection DD method was incorporated into a gradient-based optimization algorithm to maximize the Q-factor of a nanophotonic resonator.
The savings in computational effort are particularly significant for optimization with respect to several parameters, which is a common task in nanophotonics. Therefore, we expect the approach to prove especially useful when many sensitivities are to be calculated. The Riesz projection DD method can not only be applied to problems in nanophotonics, but to any resonance problem.

Data availability
All relevant data generated or analyzed during this study are included in this published article. Tabulated data files are included in a corresponding data publication 35 . Fig. 4 Optimization of a nanophotonic resonator. The optimized nanophotonic resonator with a sketch of the electric field intensity jẼj 2 corresponding to the high-Q resonance is shown. The high-Q resonance has a Q-factor of Q = 292. The materials of the nanoresonator are the same as for the reference structure in Fig. 2a.