Finding exceptional points in realistic systems using full-wave simulations

Here, we present an approach to finding exceptional points using the finite-element method. Using this method, we demonstrated exceptional points in 2D dimers of infinite cylinders and infinite parallelepipeds. The results agree well with the analytical coupled-dipole model, however a deviation due to the contribution of higher multipoles, is present. Our approach can be applied to three-dimensional particles as well.


Introduction
Adding gain and loss to optical systems leads to many counterintuitive effects, such as degeneration of field distributions of two or more optical modes into one at an exceptional point (EP) [1]. This leads to unique effects, such as unidirectional reflectionless propagation, asymmetric backscattering and directional emission, non-reciprocal light transmission and many others (see, e.g. [1], [2] and references therein). Moreover, the unique nature of EP degeneracy lifting found its application in sensing [3]. In optical systems, EPs have so far been intensively studied in systems larger than wavelength, such as waveguides [4], resonators supporting whispering gallery modes [5,3] and layered structures [2]. Due to the high Q-factor of these structures, their optical modes are coupled through evanescent fields, which allows to treat such (a) Gain Loss   [3]. Microparticles, in contrast, have very large radiative decay that leads to a mixed coupling of the eigenmodes consisting of an evanescent and radiative parts, which severely complicates the analysis [6,7]. Nevertheless, realizing the exceptional point regime in a subwavelength structure or a particle is promising for creation of metasurfaces with unique properties.
Full-wave simulations are the most convenient method of finding the eigenmodes in a structure consisting of complex shapes which cannot be treated analytically, such as finite cylinders, pyramids, cubes or droplets. However, at EP the field distributions of two or more modes become degenerate, as do the corresponding frequencies. This makes a simple eigenmode analysis relying on full-wave simulations unapplicable to finding the system parameters corresponding to the EP: even if we somehow guess the right parameter values, a full-wave simulation without proper modifications would only find one mode without any proof that it is degenerated.
Here, we propose an algorithm relying on eigenmode analysis using finite-element method (FEM), that allows to find the exceptional point in a 2D space of system parameters (e.g., material gain/loss and distance between the resonators) assuming that a single second-order EP exists in the range of interest. We demonstrate our results on 2D dimers consisting of infinite cylinders and parallelepipeds, and compare the results to an analytic model.

Exceptional point search with full-wave simulations
To find the exceptional point, we utilize the specific topology that eigenfrequency sheets display in the parametric space around the EP, which is shown in Fig. 1. We introduce the following   Fig. 1, or possibly an odd number of crossings. If there is an even number of crossings, or no crossings at all, there is no EP inside the circle.
Using this criterion, we implement a recursive algorithm to enclose EP within the desired accuracy. First, we select a point A in the parametric space and find the two eigenmodes which are supposed to coalesce, and the corresponding eigenfrequencies. Then we choose an initial grid cell ABCD, which contains the exceptional point. To achieve that, we encircle the cell twice, and if the crossing number criterion indicates the abscence of EP, we extend the cell. Afterwards, the initial cell is split into four subcells AKON , KBLO, OLCM and N OM D, as shown in Fig. 2. Then we encircle each of the subcells twice and use the crossing number criterion to check which of the subcells contains the EP. Then we take this subcell as the initial cell and continue the procedure recursively until the desired accuracy is met.
Let us now consider the EP encircling procedure in more detail. The eigenfrequencies (ω 1A , ω 2A ) at the start point A in the parametric space are already known from the initial the eigenmode analysis. At a point A which lies close to the point A, we use FEM to compute a single eigenmode using a known shift value for the eigenfrequency, searching for the mode with the eigenfrequency closest to the shift value. Using ω 1A as the shift value produces the new eigenfrequency ω 1A , and using ω 2A produces ω 2A . The set (ω 1A , ω 2A ), which corresponds to the point A , can in turn be used to find the eigenfrequencies at A and so on. Provided that the points A, A , A , . . . lie close enough to each other, the frequencies ω 1A , ω 1A , ω 1A , . . . all correspond to mode 1, and the frequencies ω 2A , ω 2A , ω 2A , . . . all correspond to mode 2. Our studies showed that splitting each of the segments AB, BC, CD, AD into 50 points was sufficient to prevent the modes from mixing up.

Results
First, we consider a dimer of infinite cylinders shown in Fig. 2a which can be treated analytically. In our previous work [7], we have shown using the dipole approximation that EP can be achieved at ω EP = 89.9 · 10 6 s −1 in a dimer made of cylinders with radius 1 m, center-to-center distance 8.28 m and refractive indices n 1 = 7.86 − 0.27i and n 2 = 7.86 − 0.088i in TE polarization. We set n 2 = n + in , where n and n become the varied system parameters, i.e., form the 2D parameter space. Using n 2 = 7.76 − 0.19i as the start point, we utilize the recursive algorithm with the exit condition |ω 1 − ω 2 | < 10 5 s −1 . Fig 2b,c show the real and imaginary parts of eigenfrequencies ω 1 and ω 2 obtained by our algorithm. A crossing-to-anticrossing transition can be observed as we go from the upper left to the lower right quadrant through the exceptional point at n 2 = 7.97 − 0.044i. The obtained value for n 2 is rather close to the analytical result n 2 = 7.86 − 0.088i. The EP eigenfrequency ω EP = (89.4 + 0.3) · 10 6 s −1 obtained with FEM has a small imaginary part, as opposed to the purely real value predicted by the coupled-dipole model. We attribute the observed deviation to the contribution of higher multipoles. Now let us consider a dimer of infinite parallelepipeds with a 1 m square cross-section separated by an 1 m gap, as shown in the inset of Fig. 3. We choose TE polarization and set n 1 = 7.86 − 0.088i and n 2 = n + in , where n and n are the varied parameters. Fig 3a,b show the real and imaginary parts of eigenfrequencies obtained by our algorithm for this dimer. As in the previous case, a crossing-to-anticrossing transition can be observed as we go from the upper left to the lower right quadrant through the exceptional point at n 2 = 7.81 − 0.14i. The EP eigenfrequency is ω EP = [167.0 − 0.3i] · 10 6 s −1 .

Conclusion
We have developed an approach to finding the system parameters corresponding to an exceptional point using the finite-element method. Using this method, we demonstrated exceptional points in 2D dimers of infinite cylinders and infinite parallelepipeds. The results for infinite cylinders agree well with the analytical coupled-dipole model, however a deviation due to the contribution of higher multipoles, is present. In this work we considered only two-dimensional resonators for simplicity, however our approach can be applied without any modification to three-dimensional particles as well.