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

### 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 steady-state regime, light–matter interactions can be described by the time-harmonic Maxwell’s equations in second-order form,

$$\nabla \times {\mu }_{0}^{-1}\nabla \times {{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}},{\omega }_{0})-{\omega }_{0}^{2}\epsilon ({{{{{{{\bf{r}}}}}}}},{\omega }_{0}){{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}},{\omega }_{0})={{i}}{\omega }_{0}{{{{{{{\bf{J}}}}}}}}({{{{{{{\bf{r}}}}}}}}),$$

(1)

where \({{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}},{\omega }_{0})\in {{\mathbb{C}}}^{3}\) is the electric field, \({{{{{{{\bf{r}}}}}}}}\in {{\mathbb{R}}}^{3}\) is the position, \({\omega }_{0}\in {\mathbb{R}}\) is the angular frequency, and \({{{{{{{\bf{J}}}}}}}}({{{{{{{\bf{r}}}}}}}})\in {{\mathbb{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 \(\tilde{{{{{{{{\bf{E}}}}}}}}}({{{{{{{\bf{r}}}}}}}})\in {{\mathbb{C}}}^{3}\) and the corresponding complex-valued eigenfrequency by \(\tilde{\omega }\in {\mathbb{C}}\). The *Q*-factor of a resonance is defined by

$$Q=\frac{{{{{{{{\rm{Re}}}}}}}}(\tilde{\omega })}{-2{{{{{{{\rm{Im}}}}}}}}(\tilde{\omega })}$$

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 \(\tilde{{{{{{{{\bf{E}}}}}}}}}\) and in the corresponding eigenfrequency \(\tilde{\omega }\), which describes the sensitivity of \(\tilde{{{{{{{{\bf{E}}}}}}}}}\) and \(\tilde{\omega }\) with respect to the parameter *p*. To compute the eigenfrequency sensitivity, we introduce a contour-integral-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 Riesz-projection-based approach for computing eigenfrequency sensitivities, which are the partial derivatives of the eigenfrequency, we consider the electric field \({{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}},{\omega }_{0}\in {\mathbb{R}})\) as a solution of Eq. (1) and \({{{{{{{\bf{E}}}}}}}}({{{{{{{\bf{r}}}}}}}},\omega \in {\mathbb{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 \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\) be a physical observable, where \({{{{{{{\mathcal{L}}}}}}}}:{{\mathbb{C}}}^{3}\to {\mathbb{C}}\) is a linear functional, and \(\tilde{C}\) be a contour enclosing the pole \(\tilde{\omega }\) of the order *m* and no other poles. Then, the Laurent expansion of \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\) about \(\tilde{\omega }\) is given by

$${{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})=\, \mathop{\sum }\limits_{k=-m}^{\infty }{a}_{k}{(\omega -\tilde{\omega })}^{k},{{{{{{{\rm{where}}}}}}}}\\ {a}_{k}(\tilde{\omega })=\, \frac{1}{2\pi {{{{{\rm{i}}}}}}}\oint \limits_{\tilde{C}}\frac{{{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}}(\omega ))}{{(\omega -\tilde{\omega })}^{k+1}}\,{{{{{\rm{d}}}}}}\omega \in {\mathbb{C}}.$$

(2)

The coefficient \({a}_{-1}(\tilde{\omega })\) is the so-called residue of \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\) at \(\tilde{\omega }\). Using Eq. (2) with the assumption that \(\tilde{\omega }\) has the order *m* = 1 and applying Cauchy’s integral formula yield

$$\oint \limits_{\tilde{C}}\omega {{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega =\oint \limits_{\tilde{C}}\frac{\omega }{\omega -\tilde{\omega }}{a}_{-1}(\tilde{\omega })\,{{{{{\rm{d}}}}}}\omega =\tilde{\omega }\oint \limits_{\tilde{C}}{{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega ,$$

where, due to the closed integral in the complex plane, the regular terms in the expansion vanish. With this, the eigenfrequency \(\tilde{\omega }\) is given by

$$\tilde{\omega }=\frac{\oint \limits_{\tilde{C}}\omega {{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega }{\oint \limits_{\tilde{C}}{{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega }.$$

(3)

The contour integrals in this equation are essentially Riesz projections for \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\) and \(\tilde{C}\)^{24}. Partial differentiation with respect to a parameter *p* directly gives the desired expression for the eigenfrequency sensitivity,

$$\frac{\partial \tilde{\omega }}{\partial p}=\, \left(\frac{\partial u}{\partial p}v-u\frac{\partial v}{\partial p}\right)\frac{1}{{v}^{2}},\;{{{{{{{\rm{where}}}}}}}}\\ u=\, \oint \limits_{\tilde{C}}\omega {{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega ,v=\oint \limits_{\tilde{C}}{{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{E}}}}}}}})\,{{{{{\rm{d}}}}}}\omega ,\\ \frac{\partial u}{\partial p}=\, \oint \limits_{\tilde{C}}\omega {{{{{{{\mathcal{L}}}}}}}}\left(\frac{\partial {{{{{{{\bf{E}}}}}}}}}{\partial p}\right)\,{{{{{\rm{d}}}}}}\omega ,\frac{\partial v}{\partial p}=\oint \limits_{\tilde{C}}{{{{{{{\mathcal{L}}}}}}}}\left(\frac{\partial {{{{{{{\bf{E}}}}}}}}}{\partial p}\right)\,{{{{{\rm{d}}}}}}\omega .$$

(4)

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 \(\tilde{{{{{{{{\bf{E}}}}}}}}}\) and its sensitivity \(\partial \tilde{{{{{{{{\bf{E}}}}}}}}}/\partial p\) can be represented by the contour integrals

$$\tilde{{{{{{{{\bf{E}}}}}}}}}=\oint \limits_{\tilde{C}}{{{{{{{\bf{E}}}}}}}}\,{{{{{\rm{d}}}}}}\omega \;{{{{{{{\rm{and}}}}}}}}\;\frac{\partial \tilde{{{{{{{{\bf{E}}}}}}}}}}{\partial p}=\oint \limits_{\tilde{C}}\frac{\partial {{{{{{{\bf{E}}}}}}}}}{\partial p}\,{{{{{\rm{d}}}}}}\omega ,$$

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 *A**E* = *f*, where \(A\in {{\mathbb{C}}}^{n\times n}\) is the system matrix, \(E\in {{\mathbb{C}}}^{n}\) is the scattered electric field in a finite-dimensional FEM basis, and \(f\in {{\mathbb{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 time-domain method, the Fourier modal method, or the boundary element method^{12,28}.

In order to calculate eigenfrequencies \(\tilde{\omega }\) and their sensitivities \(\partial \tilde{\omega }/\partial {p}_{i}\) with respect to parameters *p*_{i}, the electric fields **E** and their sensitivities ∂**E**/∂*p*_{i} are computed for complex frequencies \(\omega \in {\mathbb{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

$$\frac{\partial E}{\partial {p}_{i}}={A}^{-1}\left(\frac{\partial f}{\partial {p}_{i}}-\frac{\partial A}{\partial {p}_{i}}E\right).$$

(5)

In a first step, instead of directly computing *A*^{−1}, an *L**U*-decomposition of *A*, which can be seen as the matrix variant of Gaussian elimination, is computed to efficiently solve the linear system *A**E* = *f*. In the FEM context, this step is usually a computationally expensive step in solving scattering problems, so reusing an *L**U*-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* = *L**U*, *E*, ∂*A*/∂*p*_{i}, and ∂*f*/∂*p*_{i} are used to compute ∂*E*/∂*p*_{i} in Eq. (5). The *L**U*-decomposition 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 \({{{{{{{\mathcal{L}}}}}}}}({{{{{{{\bf{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 \(\tilde{\omega }\) and its sensitivity \(\partial \tilde{\omega }/\partial {p}_{i}\) can be calculated without solving resonance problems \(\nabla \times {\mu }^{-1}\nabla \times \tilde{{{{{{{{\bf{E}}}}}}}}}-{\tilde{\omega }}^{2}\epsilon \tilde{{{{{{{{\bf{E}}}}}}}}}=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 \({\epsilon }_{{{{{{{{\rm{r}}}}}}}}}({\omega }_{0})={\epsilon }_{\inf }-{\omega }_{{{{{{{{\rm{p}}}}}}}}}^{2}/({\omega }_{0}^{2}+{{{{{\rm{i}}}}}}{\omega }_{0}\gamma )\) is chosen, where \({\epsilon }_{\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 \(\tilde{\omega }\) of the high-*Q* resonance. Figure 2b shows the complex frequency plane with the calculated eigenfrequency, \(\tilde{\omega }=(1.17309-0.00296i)\times 1{0}^{15}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}\), and the corresponding circular integration contour \(\tilde{C}\) 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 contour \(\tilde{C}\), 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 \(\tilde{\omega }\) 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 \(|{\tilde{{{{{{{{\bf{E}}}}}}}}}}|^{2}\) is shown in Fig. 2c. The eigenmode \(\tilde{{{{{{{{\bf{E}}}}}}}}}\) has the quantum number *m* = 0 and is strongly localized in the vicinity of the nanoresonator.

Next, the eigenfrequency sensitivities \(\partial \tilde{\omega }/\partial {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 \(\partial \tilde{\omega }/\partial {p}_{i}\approx \left(\tilde{\omega }({p}_{i}+\delta {p}_{i})-\tilde{\omega }({p}_{i}-\delta {p}_{i})\right)/\left(2\delta {p}_{i}\right)\) 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 \({(\delta {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 contour-integral-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 \(\tilde{\omega }\) and then we add the sensitivities, starting with \(\partial \tilde{\omega }/\partial {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 *L**U*-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 *L**U*-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 \({\tilde{\omega }}_{{{{{{{{\rm{opt}}}}}}}}}=(1.176897-0.002015i)\,\times \,1{0}^{15}\,{{{{{{{{\rm{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.