Computing eigenfrequency sensitivities near exceptional points

Exceptional points are spectral degeneracies of non-Hermitian systems where both eigenfrequencies and eigenmodes coalesce. The eigenfrequency sensitivities near an exceptional point are significantly enhanced, whereby they diverge directly at the exceptional point. Capturing this enhanced sensitivity is crucial for the investigation and optimization of exceptional-point-based applications, such as optical sensors. We present a numerical framework, based on contour integration and algorithmic differentiation, to accurately and efficiently compute eigenfrequency sensitivities near exceptional points. We demonstrate the framework to an optical microdisk cavity and derive a semi-analytical solution to validate the numerical results. The computed eigenfrequency sensitivities are used to track the exceptional point along an exceptional surface in the parameter space. The presented framework can be applied to any kind of resonance problem, e.g., with arbitrary geometry or with exceptional points of arbitrary order.


I. INTRODUCTION
Resonance phenomena play a crucial role in the field of photonics.The interaction of light with photonic nanoresonators leads to a strong increase of the electromagnetic fields.This effect can be used for, e.g., probing single molecules with ultrahigh sensitivity [1], designing nanoantennas with a tailored directivity and large spontaneous emission rate [2], or realizing efficient singlephoton sources [3].Resonances are characterized by loss mechanisms, such as damping or open boundaries, yielding non-Hermitian systems with complex-valued eigenfrequencies [4,5].In general, the resonances are numerically computed by solving the source-free Maxwell's equations [6][7][8].For the investigation of photonic systems, not only the resonances are of interest, but also their partial derivatives, the so-called sensitivities, with respect to certain system parameters.These enable a better understanding of the underlying physical effects [9] and an efficient optimization of corresponding photonic devices [10,11].
The eigenfrequencies near an EP show a characteristic complex-root topology [13].Consequently, near the EP, small parametric changes of the system are amplified to a strong response of the eigenfrequencies.This EP-enhanced eigenfrequency sensitivity is challenging in experiments and may also spoil the accuracy of numerical simulations.On the other hand, calculation of the eigenfrequency sensitivities near an EP is crucial for applications such as sensors or optimization schemes.Hence, accurate and efficient numerical methods for photonic systems operating near an EP are of fundamental interest in research and engineering.
In this work, we address this challenge by combining contour integration with algorithmic differentiation (AD).Contour integration gives access to eigenfrequencies in non-Hermitian systems by solving Maxwell's equations with source terms at frequencies on a contour in the complex frequency plane.For the corresponding scattering problems, AD is exploited so that the calculation of eigenfrequency sensitivities is numerically accurate and efficient.Although the eigenfrequency sensitivities diverge near an EP, which is sketched in Fig. 1, the presented framework allows to capture the EP-enhanced eigenfrequency sensitivity.We apply the framework to an optical microdisk cavity with an EP of second order.The two coalescing eigenfrequencies and the corresponding sensitivities are computed near the EP.The sensitivities are further used for an optimization scheme to track the EP along a surface in the underlying parameter space.Such exceptional surfaces are used to combine the sensitivity of EPs with robustness against fabrication tolerances [39][40][41].

II. THEORETICAL BACKGROUND
In photonics, in the steady-state regime, light scattering by an open system can be described by the time- where E(r, ω) ∈ C 3 is the electric field, J(r) ∈ C 3 is a source term corresponding to an optical source, ω ∈ C is the angular frequency, and r ∈ R 3 is the spatial position.For optical frequencies, the permeability µ(r, ω) = µ r (r, ω)µ 0 typically equals the vacuum permeability µ 0 .The permittivity ϵ(r, ω) = ϵ r (r, ω)ϵ 0 , where ϵ r (r, ω) is the relative permittivity and ϵ 0 is the vacuum permittivity, describes the material dispersion and the spatial distribution of material.Problems given by Eq. ( 1) are called scattering problems.Resonance problems are given by Eq. ( 1) but without a source term.Note that, in the following, we consider non-Hermitian systems based on open boundaries.However, the inclusion of damping as loss channel is also possible.
A. Contour integration for computing eigenfrequency sensitivities near exceptional points Electromagnetic quantities q(E(r, ω)) ∈ C, such as reflection or transmission coefficients, are measured for real excitation frequencies ω ∈ R. To compute complexvalued eigenfrequencies and their sensitivities in non-Hermitian systems, we consider the analytical continuation of q(E(r, ω ∈ R)) into the complex frequency plane ω ∈ C, which we denote by q(ω) as a short notation of q(E(r, ω ∈ C)) [42].The L eigenfrequencies inside a chosen contour C are given by the eigenvalues ω l of the generalized eigenproblem [43][44][45] where Ω = diag(ω 1 , . . ., ω L ) is a diagonal matrix containing the eigenvalues, the columns of the matrix X ∈ C L×L are the eigenvectors, and are Hankel matrices with the contour-integral-based elements The second equality results from applying Cauchy's residue theorem, where a l are the residues of q(ω) corresponding to the eigenfrequencies ω l , and it holds for simple eigenfrequencies.Once the elements s k are computed, there is not only access to the eigenfrequencies by solving the eigenproblem given in Eq. ( 2), but also to the eigenfrequency sensitivities ∂ω l /∂p, where p is some parameter.Computing the derivatives in Eq. ( 3) yields the linear system of equations for the unknowns ∂ω l /∂p and ∂a l /∂p, where ω l and a l are known due to solving Eq. ( 2) and using diag(a 1 , . . ., a L ) = X T HX(V T X) −2 with the Vandermonde matrix [42] .
The eigenfrequencies at EPs are not simple.However, experimental and numerical realizations show always an eigenfrequency splitting due to fabrication inaccuracies [39] or numerical approximation [46].For this practical reason, we can consider Eq. ( 4) to compute eigenfrequency sensitivities in EP-based systems.
The elements s k of the Hankel matrices and their sensitivities ∂s k /∂p are obtained by numerical quadrature [47] of the contour integrals in Eq. (3) and Eq. ( 4), respectively, where q(ω) and ∂q(ω)/∂p are computed for complex frequencies on the integration contour C. The quantities q(ω) and ∂q(ω)/∂p only have to be computed once for each integration point and then all contour integrals can be evaluated.The integrands differ only in the weight functions ω k .Information on the numerical realization of the contour integration can be found in the Refs.[48,49].

B. Structure exploiting algorithmic differentiation
To evaluate q(ω) and ∂q(ω)/∂p on the integration contour, the electric field E(r, ω) and the corresponding sensitivity ∂E(r, ω)/∂p have to be computed.For this, we solve Maxwell's equation given in Eq. ( 1) using the finite element method (FEM) with the implementation of the software package JCMsuite [50].The spatial discretization of Eq. (1) yields the linear system of equations AE = f , where A ∈ C n×n is the FEM system matrix, E ∈ C n is the scattered electric field in a finitedimensional FEM basis, and f ∈ C n realizes the source term.
For the computation of ∂E(r, ω)/∂p, an approach based on directly using the FEM system matrix is applied [11,51,52].With this direct differentiation approach, the sensitivities of scattering solutions can be computed by Instead of directly computing A −1 , an LU -decomposition of A is calculated to efficiently solve the linear system AE = f .The LU -decomposition is used to obtain E and also ∂E/∂p.In the FEM context, an LU -decomposition is usually a computationally expensive step, so reusing it significantly reduces computational effort [11].
Several approaches could be used to evaluate the sensitivities of the system matrix, ∂A/∂p, and of the source term, ∂f /∂p.The required derivative information can be provided analytically by means of corresponding computer algebra systems if the considered functional dependence is not too complicated.Alternatively, one may employ finite differences to approximate the gradient information.However, the resulting computational effort scales linearly with the number of parameters p [11].Furthermore, finite differences only yield an approximation of the required derivative information, a fact that may cause problems if the accuracy of the derivatives is essential.For this reason, we propose AD to provide the required sensitivities within working accuracy in an efficient way [53].Suppose a function F : R n → R m , y = F (p), is given in a computer language like C or C++.Then, the evaluation of F (p) can be decomposed into socalled elementary functions, the derivatives of which are well known and easy to implement in a software package.The basic differentiation rules, such as the product rule, the quotient rule, etc., can be applied to each statement of the given code segment to calculate the overall derivatives.Hence, exploiting the chain rule yields the derivatives of the whole sequence of statements with respect to the input variables.
One distinguishes two basic modes of AD, namely the forward mode and the reverse mode.The former one evaluates derivatives together with the function evaluation.In mathematical terms, one obtains, for a given direction ṗ, the Jacobian-vector product ẏ = F ′ (p) ṗ.Hence, for a unit vector ṗ = e i ∈ R n , one obtains the ith column of the Jacobian ∇F (p).When applying the reverse mode of AD, one propagates derivative information from the dependents y to the independents p after the evaluation of the function F (p).This yields, for a given weight vector ȳ, the vector-Jacobian product p = ȳ⊤ F ′ (p).Similar to the forward mode, for a unit vector ȳ = e j ∈ R m , one obtains the jth row of the Jacobian ∇F (p).
Over the past decades, extensive research activities led to a thorough understanding and analysis of these two basic modes of AD, where the complexity results with respect to the required runtime are based on the operation count O F , i.e., the number of floating point operations required to evaluate F (p).The forward mode of AD yields one column of the Jacobian ∇F for no more than three times O F [54,55].Using the reverse mode of AD, one row of ∇F is obtained for no more than four times O F , also see Ref. [54].It is important to note that this bound for the reverse mode is completely independent of the number n of input variables.Hence, if m = 1, the gradient of the then scalar-valued function F can be calculated for four function evaluations.This observation is called cheap gradient result and is used extensively for derivative-based optimization approaches.However, the reverse mode requires the knowledge of intermediate results computed during the function evaluation.Therefore, the basic implementation of the reverse mode leads to a memory requirement that is proportional to the operation count O F .For a considerable amount of applications, this fact does not cause any problems.For problems of larger scale, checkpointing strategies have been developed, see, e.g., Ref. [56].Here, instead of storing all intermediates, only a few of them are recorded.Subsequently, the missing intermediate values are recomputed using the data stored in the checkpoints.Hence, these checkpointing strategies seek for a compromise between storing and recomputing data.Besides the theoretical foundation, numerous AD tools have been developed, e.g., CppAD [57], ADOL-C [58], and Tapenade [59].Due to the language features, the implementation of the AD packages is based on source transformation for Fortran codes and operator overloading for C or C++ codes.
When applying these AD-tools in a black box fashion to large-scale simulation codes, one usually does not observe the theoretical runtime factors mentioned above, e.g., due to memory issues.Therefore, a structure exploiting approach, where AD is only applied to relevant parts of the code, as used in this work, is very often beneficial for the efficient calculation of the derivatives.This includes also the appropriate choice of the mode of AD.For the application considered here, where, e.g., derivatives of a matrix with respect to a parameter vector are required, the forward mode is the method of choice.

III. APPLICATION TO A MICRODISK CAVITY
We apply the presented framework to a twodimensional optical microdisk cavity with two concentric layers of different refractive indices embedded in an open environment [60].The system is sketched in Fig. 2(a).The interest lies in the electric fields which are perpendicular to the cavity plane (TM polarization).This simplifies the time-harmonic Maxwell's equations (1) to the scalar Helmholtz equation where E z m,l (x, y) is the z-component of the l-th eigenmode with the azimuthal mode number m ∈ Z, the refractive index n(x, y) = √ ϵ r µ r describes the material in the region of interest, ω m,l ∈ C is the eigenfrequency, and c is the speed of light.
Computing eigenfrequencies and their sensitivities near the EP.(a) Complex frequency plane with integration contour C. The center of the contour is at ω EP and the radius is given by 0.01 × Re(ω EP ).Four eigenfrequencies from Fig. 3(a) are shown, where ω A 8,1 and ω A 8,2 or ω B 8,1 and ω B 8,2 are the result of a contour integration with a different R1 each case.(b) Relative error |ω8,1 −ωqex|/|ωqex| with respect to the number of integration points on the contour C, where ωqex is the quasiexact solution computed with 64 integration points.The maximum side length of the FEM mesh is 0.05 × (c) Relative error with respect to maximum side length of the FEM mesh, where ωqex is the quasiexact solution computed with a maximum side length 0.0125 × R. The number of integration points is 16.Legend as before.

A. Semi-analytical solution
Semi-analytical results for Maxwell's equations can be obtained in the case of two concentric dielectric disks [61].As the system is rotational invariant, the eigenfrequencies are characterized by an azimuthal mode number m ∈ Z.Therefore, combining the continuity conditions at the dielectric interfaces with the outgoing wave condition results in a conditional equation S m (ω m,l ) = 0 for the eigenfrequency ω m,l .In the case of TM polarization, the function S m (ω) is given by where k j = n j ω/c and J m , H m , H m are Besseland Hankel functions of first and second kind and order m, respectively.The complex roots of S m (ω) need to be determined numerically, e.g., by Newton's method, where an initial guess can be taken from a single disk with modified refractive index [61].With the implicit function theorem, the derivative of the eigenfrequency ω m,l with respect to a parameter p ∈ {n 1 , n 2 , n 3 , R 1 , R} can be calculated as The partial derivatives ∂S m /∂ω and ∂S m /∂p are inde-pendent of the eigenfrequency.Therefore, they can be calculated analytically with a common computer algebra system from Eq. ( 6), see also Ref. [49].

B. Numerical solution
To numerically solve Eq. ( 5), the software RPExpand [48] is applied.The eigenfrequencies ω m,l and the corresponding sensitivities ∂ω m,l /∂p with respect to a parameter p are computed by using the contourintegral-based framework.Here, the quantity q(ω) is the z-component of the electric field at the position For each point in the parameter space, the parameter n2 is shifted by 0.0025.
[x, y] = [0, 0.9 × R] resulting from an z-polarized plane wave travelling in y-direction.RPExpand contains an interface to the software JCMsuite.The approach of AD utilized within JCMsuite.Convergence with respect to the FEM parameters is ensured by refining the spatial mesh, where the degree of the polynomial ansatz functions is set to four.For computational efficiency, the mirror symmetries of the system are exploited.Further information on the numerical implementation can be found in Ref. [49].
We study an EP of second order [60].Two eigenmodes with an azimuthal mode number m = 8 and the corresponding eigenfrequencies coalesce at the parameter combination n 1 = 3.1239791, n 2 = 1.5, and R 1 = 0.4970147 × R. The electric field intensity of a corresponding eigenmode E z 8,1 is shown in Fig. 2(b).The eigenfrequency is given by ω 8,1 ≈ ω 8,2 = (6.96185− 0.089761i)×c/R.In the following, for a simpler notation, we denote this specific eigenfrequency by ω EP .
The eigenfrequencies ω 8,l and their sensitivities ∂ω 8,l /∂R 1 near the EP are investigated.Figure 3(a) shows the trajectories of the real parts of the eigenfrequencies ω 8,l when the parameter R 1 is varied and n 1 and n 2 are fixed.For each radius R 1 , a contour integral is computed leading to two eigenfrequencies.The superscript letters A and B describe two points in the parameter space that differ in the radius R 1 .When R 1 is increased, i.e., going from ω A 8,1 to ω B 8,1 or ω A 8,2 to ω B 8,2 , then the two eigenfrequencies converge.Figure 3(b) shows the imaginary parts and Fig. 3(c) and Fig. 3(d) show the corresponding sensitivities.The results can be compared with the semi-analytical solution ω an .For ω B 8,1 , the relative error |Re(ω B 8,1 −ω an) /Re(ω an) | is smaller than 7×10 −7 and, for the corresponding imaginary part, the relative error is smaller than 5 × 10 −5 .The relative error for the underlying sensitivities is smaller than 2 × 10 −3 for real and imaginary part.It can be observed that the sensitivities diverge in the vicinity of the EP.The divergence can be explained by the square root dependence of the eigenfrequencies on perturbations when the system is near the EP [13].
For all calculations corresponding to Fig. 3, the contour C shown in Fig. 4(a) with 16 integration points and a maximum side length of the FEM mesh of 0.05 × R are used.Figure 4(a) also sketches the locations of the eigenfrequencies ω A 8,1 , ω A 8,2 , ω B 8,1 , and ω B 8,2 with the trajectories when R 1 is varied in the direction of the parameter for the EP.To demonstrate the accuracy of the approach, Fig. 4 4(c) shows the relative error with respect to the maximum side length of the FEM mesh.It can be observed that the relative error for the sensitivities is smaller for ω A 8,1 , which lies further away from the EP.The closer the two calculated eigenfrequencies lie to each other, the less accurate the results are.With only 16 integration points, it is possible to obtain a maximum relative error of smaller than 10 −5 for the eigenfrequency sensitivities ∂ω 8,1 /∂R 1 for both points A and B in the parameter space.Using a maximum side length of the FEM mesh of 0.05 × R is sufficient to achieve a maximum relative error of smaller than 3 × 10 −3 .

C. Exceptional surface
We track the EP in the parameter space by using Newton's method with the sensitivities ∂ω 8,l /∂n 1 and ∂ω 8,l /∂R 1 .We are looking for the zeros of the squared eigenfrequency splitting (ω 8,1 − ω 8,2 ) 2 .The reason for choosing the square of the splitting is the square root dependency of the eigenfrequencies near the EP.The starting point is the parameter combination corresponding to ω B 8,1 and ω B 8,2 shown in Fig. 3(a), given by n 1 = 3.1239791, n 2 = 1.5, and R 1 = 0.497004557 × R. With a fixed n 2 = 1.5, a certain tolerance, and a few iterations of Newton's method, the two eigenfrequencies ω EP 8,1 = (6.961993− 0.089638i) × c/R and ω EP 8,2 = (6.961996− 0.089642i) × c/R with the parameters n 1 = 3.123979246 and R 1 = 0.497014753 × R are obtained.Then, the parameter n 2 is shifted by a constant value of 0.0025 and Newton's method is restarted with the shifted n 2 together with the previously calculated values for the parameters n 1 and R 1 , and with a new integration contour centered at the previously computed eigenfrequency.The results of this procedure are shown in Fig. 5, where Fig. 5(a) contains the trajectory of ω EP and Fig. 5(b) shows the corresponding parameter combinations.
Although the sensitivities diverge near the EP, the numerical calculations of the sensitivities are still accurate enough for Newton's method to work very well.Therefore, the exceptional surface is identified as a curve in the parameter space spanned by R 1 and n 1 .Figure 5 also contains the semi-analytical solution and it agrees with the numerical results.

IV. CONCLUSION
Eigenfrequency sensitivities near EPs were investigated.For this purpose, we developed a framework based on contour integration with AD.We applied a numerical implementation of the framework to an optical microdisk cavity with an EP of second order.It was shown that the eigenfrequency sensitivities near the EP can be captured accurately and efficiently.The sensitivities were applied to track the EP in the parameter space.A semianalytical solution for the system was derived and used to validate the numerical results.
The presented contour-integral-based framework cannot only be applied to a benchmark problem as studied in this work, but also to any resonant system, e.g., with complex geometry or with an EP of arbitrary order.The combination of contour integration with AD makes the framework very general since the approach is essentially based on solving scattering problems.The treatment of the corresponding linear systems of equations is a stan-dard task for state-of-the-art software packages used in the field of computational physics.
Computing eigenfrequency sensitivities near EPs can aid in optimizing setups for EP-based sensors.For example, exceptional surfaces [62] in higher dimensional parameter spaces can be identified, which can then be used to optimize an EP along such a surface to maximize sensitivity to a targeted perturbation while minimizing the response to other perturbations that may arise from fabrication tolerances or noise.For such an optimization scheme, the numerically accurate and efficient computation of eigenfrequency sensitivities with respect to specific system parameters is crucial.

DATA AND CODE AVAILABILITY
Supplementary data tables and source code for the numerical experiments for this work can be found in the open access data publication [49].Links to further AD tools and comprehensive information on AD can be found at www.autodiff.org.

FIG. 1 .
FIG.1.Eigenfrequencies ω l and their sensitivities ∂ω l /∂p in the complex frequency plane with respect to a parameter p of an optical microdisk cavity.A specific change in the parameter cause the eigenfrequencies and the associated eigenmodes to coalesce at the EP.The sensitivities of the eigenfrequencies diverge near the EP.

5 FIG. 5 .
FIG. 5. Tracking the EP in the parameter space.(a) Numerical and semi-analytical results for the trajectory of the eigenfrequency ω EP .The arrow indicates the eigenfrequency ω EP from Ref. [60].(b) Corresponding parameter combinations.For each point in the parameter space, the parameter n2 is shifted by 0.0025.
(b) shows the relative error of the absolute value of ω A 8,1 , ∂ω A 8,1 /∂R 1 , ω B 8,1 , and ∂ω B 8,1 /∂R 1 with respect to the number of integration points on the contour C. Figure