Hierarchical construction of higher-order exceptional points

The realization of higher-order exceptional points (HOEPs) can lead to orders of magnitude enhancement in light-matter interactions beyond the current fundamental limits. Unfortunately, implementing HOEPs in the existing schemes is a rather difficult task, due to the complexity and sensitivity to fabrication imperfections. Here we introduce a hierarchical approach for engineering photonic structures having HOEPs that are easier to build and more resilient to experimental uncertainties. We demonstrate our technique by an example that involves parity-time (PT) symmetric optical microring resonators with chiral coupling among the internal optical modes of each resonator. Interestingly, we find that the uniform coupling profile is not required to achieve HOEPs in this system -- a feature that implies the emergence of HOEPs from disorder and provides resilience against some fabrication errors. Our results are confirmed by using full-wave simulations based on Maxwell's equation in realistic optical material systems.


Introduction-
The notion of non-Hermitian photonics [1][2][3][4], together with its central concept of exceptional points (EPs) [5] have attracted considerable attention in the past few years. At an EP, two or more eigenvalues and their corresponding eigenstates/eigenvectors coalesce. As a result, the dimensionality of the eigenspace is reduced. An interesting feature of systems operating at EPs is their strong (nonlinear) response to small perturbation [6,7]. More specifically, for an EP of order N (or EPN for short) formed by the coalescence of N different eigenvalues and associated eigenstates, the response scales as N √ , where is the perturbation strength. Clearly, for N > 1, this response 'amplifies' weak effects, which could be potentially useful for sensing applications [8][9][10][11], wireless energy transfer [12] and enhancing lightmatter interactions [13]. However, current schemes for achieving this requires the implementation of complex networks and tuning of considerable number of design parameters [10,[14][15][16][17].
In this work, we introduce a general approach for constructing tight-binding networks supporting higher-order exceptional points (HOEPs) out of arrangements that support only lower-order EPs. As an example, we apply our scheme for constructing photonic networks with HOEPs based on judicious engineering of the gain/loss profile and coupling coefficients in both real and synthetic dimensions. More specifically, this strategy allows us to start with a parity-time (PT) symmetric arrangement having an N th order EP and modify it to obtain a 2N th order EP in a straightforward fashion. Importantly, the resultant configuration requires tuning a small number of design parameters than those of conventional PT symmetric arrays with the same EP order, and hence is more robust against fabrication errors and experimental uncertainties. * Corresponding author: ganainy@mtu.edu General procedure for doubling the order of EPs-To this end, we start with two different discrete systems described by two matrix Hamiltonians H a and H b , each with dimension N × N and exhibiting an EP of order N . Moreover, we assume that both EPs correspond to the same eigenvalue. We will denote the eigenvalue and eigenvector at the EPs by (µ, a) and (µ, b), respectively, where the N -dimensional column vectors a and b are exceptional eigenvectors (i.e. the coalesced eigenvectors at the EP) of H a and H b . Let us now consider the following mathematical construction: where K is an N × N matrix, and0 is an N × N null matrix. Note that H has a dimension 2N × 2N . The eigenvalues and eigenvectors of H can be written as where the superscript T indicates matrix transpose; and c 1,2 are each N -dimensional column vectors (thus the column vector c has a dimension of 2N ). The eigenvalue problem λI − H can be then rearranged into: Eq. (2a) admits two possible solution: a nontrivial solution with an eigenvalue λ = µ and an eigenvector c 1 = a; and a trivial one with c 1 = 0. Let us first consider the trivial solution. By substituting in Eq. (2b), we obtain (λI − H b ) c 2 = 0, which has the solution λ = µ and  Fig. 1.
Having introduced the general framework, we now consider the following special case where: H a = H b ≡ H are discrete PT symmetric Hamiltonians with an eigenvector a = b = v, K ≡ diag(k 1 , k 2 , ..., k N ) is a diagonal matrix whose elements represent a unidirectional coupling from the different elements of H a to the corresponding elements in H b , as shown in Fig. 2. Without any loss of generality, we also assume that µ = 0. Under these conditions, Eq. (2b) reduces to H c 2 = −K v. If this latter equation has one solution, then the spectrum of the combined system will contain two EPs of order N each. On the other hand, if there is no solution, the spectrum will consist of only one EP of order 2N . To investigate these two situations in more detail, we express H in its Jordan canonical form as defined by the similarity transformation H J = SHS −1 with S being the mapping matrix. The explicit form of H J is: In these bases K J = SKS −1 , c J 2 = S c 2 and v J = S v. Accordingly: If the diagonal matrix elements of K are identical and equal to k, the combined structure will respect PT symmetry, and we obtain K J = K = kI where I is the N ×N identity matrix. It follows that Eq. (4) is undetermined with a family of solutions given by c J 2 = [x, −k, 0, ..., 0] T , where x is a free parameter. Clearly, the spectrum of H in this case does not contain an EP of order 2N . A more general argument applies when [S, K] = 0. If, however, one introduces a disorder to the coupling matrix K such that the N th element of the vector K J v J is not zero, then Eq. (4) is inconsistent and has no solution, and thus H would contain an EP of order 2N .
Application to PT J x arrays-In order to illustrate the power of the scheme presented above, we demonstrate its application to J x photonic array [18][19][20] that posses PT symmetry [14]. A tight binding network that realize an EP of order N using PT J x arrays, will consist of N sites and its Hamiltonian H N is given by [14,21]: In the above equation, g n = (2n − N − 1)γ and κ n = n(N − n)γ with n = 1, 2, ..., N , and γ is a non-Hermitian parameter (in optics it represents gain or loss depending on its sign). Obviously, an implementation of an EP of order 2N can be directly obtained by scaling the above system according to N → 2N to obtain H 2N . This however introduces more non-uniformity in the coupling and gain/loss profiles, which poses practical limitations on realizing these arrangements. On the other hand, our scheme relaxes some of these constraints by enabling the construction of an EP of order 2N out of two exact copies of H N according to Eq. (1) by substituting H a = H b = H N , as demonstrated schematically in Fig. 2. In what follows, we take K to a be a diagonal matrix and we will show that this choice can lead to an EP of order 2N . Before we proceed, we highlight the interesting observation that both H 2N and H 2N are connected via a similarity transformation. This can be demonstrated by noting that both H 2N and H 2N have the same Jordan form which we denote by H J 2N , i.e.
for some mapping matrices R and S. This connection can be also expressed in terms of discrete supersymmetry [22,23]: by setting Before we proceed further, we would like to mention that few studies have presented different routes for implementing higher order EPs. For instance, Nada et al. proposed a scheme that leads to high order EP without relying on gain/loss distribution [17]. Our work goes beyond previous studies by introducing a straightforward way of modifying an already implemented system that has N -th order EP to double the order of the EP to 2N , with the minimal number of additional design parameters, and by showing for the first time that higher-order EPs can arise due to disorder, which completely goes against the conventional wisdom in the field.
Implementation of H 2N using synthetic dimensions-In principle, the structure shown in Fig. 2 can be implemented by using 2N resonators together with optical isolators to introduce the unidirectional coupling. This however does not provide any advantage in terms of scalability or fabrication. Alternatively, we explore a different route that relies on synthetic dimensions [24][25][26][27][28]. Particularly, we consider a PT symmetric J x array made of N microresonators. Time reversal symmetry implies that each resonator supports two traveling waves at each resonant frequency, clockwise(CW) and counterclockwise (CCW). Taken together, these 2N modes can act as bases for implementing an EP of order 2N . Figure 3(a) illustrates one possible realization of the above photonic network when N = 2. The unidirectional coupling here is introduced via an evanescently coupled waveguide with an end mirror -a geometry that was recently proposed for building robust EP sensors [29], amplifiers [30], and directional absorbers [31]. Figure 3(b) depicts the equivalent arrangement in real space. Interestingly, if a mirror is introduced on each waveguide, the system will not exhibit a HOEP. In other words, the formation of the HOEP in this structure can be achieved only by breaking the PT symmetry -a rather remarkable observation.
Applications in enhanced light-matter interaction-Before we present a concrete design for the above structure with realistic material systems, we first confirm that it can be indeed used to enhance light-matter interaction. As a prototype example, we focus again on the implementation of H 4 as shown in Fig. 3 and we assume that a nanoparticle exists in the near field region of one of the resonators as shown in Fig. 3. By assuming that the nanoparticle introduces a perturbation strength , it follows from the perturbation theory of non-Hermitian operators [6,7] that an eigenfrequency splitting scaling should take place. To confirm that this is indeed the case for our structure, we consider its Hamiltonian as derived by using temporal coupled mode theory [32], which   FIG. 3. (a) A PT symmetric system consisting of two microring resonators having equal loss (−iγ), gain (iγ) and coupling to implement an EP of order two. By introducing an evanescent coupling between any of the rings (the lossy ring here) and a waveguide with an end mirror, the system realizes the Hamiltonian H4 of Eq. (6). The nanoparticle close to the gain ring is used as a small perturbation to confirm that the combined system indeed has an EP of order four in the later simulation. (b) Equivalent system in (a) in real space, detailing the loss/gain distribution and the coupling profile between the various CW/CCW components.
takes the form H 4 = H 4 + H p , with: Here ω 0 is the resonant frequency of each mode, ±γ is the gain/loss coefficients of the modes in the gain and loss resonators, which is also taken to be identical to the coupling coefficient between the two resonators in order to implement a PT symmetric dimer (when k = 0). Additionally, H p is the perturbation Hamiltonian, and = | | exp (iφ ). We remark that the above Hamiltonian is expressed in the bases [a cw , b ccw , a ccw , b cw ] T (see Fig. 3). When = 0, the above system has an EP4 with an eigenvalue ω = ω 0 . Assuming a small perturbation | | |γ|, |k|, we obtain the following asymptotic expression (see Appendix A for more details) for the eigenvalues ω n : where χ n = i n−1 4 exp(iφ )kγ 2 , and n = 1, 2, 3, 4. In other words, by just evanescently coupling the resonators to waveguides and terminating one waveguide with a mirror, the eigenvalue splitting changes from | | to 4 | |.
Similarly, one can apply the same strategy to the PT system that supports EP of order three, which has been implemented experimentally in [10], in order to obtain a modified structure supporting an EP of order six. Full-wave simulations-Finally, we confirm our results by performing full-wave analysis using finite element  7), while the dots represent actual values obtained from fullwave numerical simulations using FEM. As expected, the perturbation due to the particle causes the unperturbed eigenstate to split into four different eigenmodes. From the figures, we observe good agreement between theory and simulations. We note however that even in the absence of the particle, our simulations indicate a finite splitting between the eigenfrequencies. This is due to the fact that the waveguides themselves can introduce back-scattering between the modes. This means that our system operates in the vicinity, rather than exactly at, the EP. As we will see below, this does not have a significant impact on the ability of the structure to enhance light-matter interactions. Figure 4(c) plots Re[ω 1 − ω 0 ] as a function of using log scale. A curve fitting of the data produces the relation Re[ω 1 − ω 0 ] = A 0.24 , corresponding to a line with a slope 0.24 on the log scale, which is very close to the 0.25 value expected for an ideal EP of order four. Moreover, a comparison between the numerical and analytical values for the splitting amplitudes (i.e. the value of the constant A in the formula above and χ 1 in Eq. (7)) confirms the excellent agreement with only 4% error. For completeness, we discuss the modal profiles associated with the resonant modes of the structure investigated above in Appendix E.
Conclusions-In summary, we have introduced a new, generic approach for constructing tight-binding Hamiltonians with HOEPs out of initial arrangements with lower-order EPs. As an illustrative example, we have demonstrated the detailed application of our scheme to PT symmetric J x arrays. By focusing on PT symmetric dimer and utilizing the concept of synthetic dimensions in multimodal systems, we have presented an elegant and more robust (compared to conventional configurations) implementation of optical ring dimer that exhibits a fourth-order EP. We should note here that many interesting features of EPs, such as a system's response to perturbations as in sensors, stem from the order of EPs: the higher the better, regardless of whether it is an odd or even order EP. Mathematically speaking, there is no bound on the order of EP that can be engineered using our approach. Limitations originate from the maturity and scalability of the fabrication process of the underlined physical platform. For example, in electronics where fabrication is highly scalable and fabrication errors are minimal, one can construct EP of any even order EP for telemetry and wireless energy transfer. Scalability in photonics has always been a critical problem despite decent progress in the past few years, and in general building photonic systems with HOEP is challenging (i.e., none beyond order 3 yet). Our approach provides a methodological approach and a realistic route to build photonic systems with HOEP, alleviating some of the difficulties encountered. This in turn may have an impact on building better sensing devices (though this topic is still subject to debate, see for instance [34][35][36][37][38]), as well as photonic components with enhanced nonlinear interactions.

Appendix B: Geometric and material parameters
For the fullwave simulations, we used the following parameters: each of the rings and the waveguides has a width w = 0.25 µm, and refractive indices of 3.47 in a background refractive index n 0 = 1.44 (relevant to SiN on silica platforms); the radii of the rings are identical and taken to be R = 4.75 µm each; the edge-to-edge separation between each ring and its neighbor waveguide is h 1 = 0.25 µm while that between the two rings is h 2 = 0.4 µm; the distance between the center of the nanoparticle and the edge of the ring resonator is 0.1 µm; and the distance between the ring-waveguide junction and the mirror is L = 5.11 µm. Moreover, the mirror is implemented by using a thin layer of silver with a thickness of 100 nm, and the loss and gain are modeled by including an imaginary part to the refractive index of the two microrings, i.e. n 1,2 = 3.47 ± in i . Finally, the refractive index of the nano-particle is n p = 3.47. In order to minimize the effect of the surface roughness along the ring resonator due to numerical discretization, we set up the average mesh size equal to λ/50 (here λ = λ f /n where λ f = 1550 nm is some free space reference wavelength and n is the refractive index of the materials) and use a finer mesh size of 1 nm around the nanoparticle region.
Appendix C: Optical parameters of the structure In order to extract the numerical values of the various optical parameters, we first simulate a single microring resonator and compute its resonance frequency ω 0 , or equivalently the corresponding free space wavelength of 1548 nm.
Next, we consider a coupled resonators configuration without any waveguides. In this case, the coupling coefficient between the resonators can be obtained by calculating the resonant frequencies of the even/odd supermodes in the absence of any gain or loss. By doing so  Fig. 3 without the waveguide. The EP is located at ni = 2.9051 × 10 −4 . The coupling coefficient between the two resonators can be obtained from the data at ni = 0, which is found to be κ = 98.8 GHz.
using COMSOL software package, we find κ = 98.8 GHz. Next, we include the gain and loss in the two resonators by adding an imaginary part ±n i to their refractive indices in order to create a PT symmetric arrangement and we plot the resonant frequency as a function of n i as shown in Fig. 5. From these plots, we find that the EP occurs when n i = 2.9051 × 10 −4 , which indicate that the field amplitude gain and loss at this point is equal to the coupling between the resonators.
Finally we evaluate the decay rate from the microring to the waveguide γ w by constructing a symmetric adddrop ring resonator device and measuring the bandwidth of the transmission between the input and the output port as a function of frequency ( Fig. 6(a)). In this case, for a negligible material absorption and radiation loss, the bandwidth is given by BW = 4γ w . The FEM simulations gives the value γ w = 43.25 GHz. Thus, for an all-pass ring-waveguide structure, the resonant frequency is complex and is given by ω 0 − iγ w .
Appendix D: Perturbation coefficient due to the nano-particle Here we characterize the perturbation introduced by the nanoparticle as a function of its radius R p . To do so, we first note that a particle located in the near field of the microresonators will introduce a coupling between the CW and CCW modes, forming standing wave patterns. These new supermodes are distributed such that the particle is located at the node of the antisymmetric mode, and at the maximum of the symmetric as shown in Fig. 7(a) for a 10 nm particle. Here symmetric and antisymmetric refer to the field distribution with respect to a horizontal line that crosses the nano-particle's center. Consequently, the antisymmetric mode is not affected much by the particle while symmetric mode experiences a frequency shift of 2| |. The perturbation Hamiltonian H p in the main text, which is written in the CW/CCW bases, capture this behavior. Finally, Fig. 7(b) plots the value of | | (the frequency is red-shifted because the particle's refractive index is larger than that of the surrounding [39]) as a function of R p when the latter varies from 2 to 20 nm. In obtaining this plot, we fit the results obtained from full-wave simulations to those obtained using the Hamiltonian H test = ω 0 + ω 0 + .
Appendix E: Modal profile Figure 8 presents an example of the spatial field distribution (amplitude of the electric field) of the four eigenmodes of the structure for a particle of radius 10 nm. We observe that both eigenmodes 1 and 3 have equal power distribution in both rings, which explain their dominant real frequency splitting as demonstrated in Fig. 4(a). On the other hand, modes 2 and 4 have more intensity in the gain/loss microring, respectively which explain their dominant imaginary frequency shift as can be seen from Fig. 4(b). These numerical results thus confirm the analytical predictions and demonstrate the feasibility of our approach for practical implementations. FIG. 7. A nanoparticle in the vicinity of a microring resonator (left panel of (a)) will introduce a coupling between the CW and CCW modes. The new supermodes of the ring will have symmetric/asymmetric field distribution with respect to the particle location. These modes are depicted in the right panel of (a) for a nanoparticle of radius 10-nm located at the node of the asymmetric mode and anti-node of the symmetric mode, respectively (the rest of the parameters as identical to those used in the text). (b) By using full-wave simulations to compute the frequency splitting due to nanoparticles and fitting these results to those obtained from the Hamiltonian Htest, we obtained the perturbation | | as a function of the particle's radii, Rp. (b) plots such a relation when Rp varies from 2 to 20 nm.