General Theory of Spontaneous Emission Near Exceptional Points

We present a general theory of spontaneous emission at exceptional points (EPs)---exotic degeneracies in non-Hermitian systems. Our theory extends beyond spontaneous emission to any light--matter interaction described by the local density of states (e.g., absorption, thermal emission, and nonlinear frequency conversion). Whereas traditional spontaneous-emission theories imply infinite enhancement factors at EPs, we derive finite bounds on the enhancement, proving maximum enhancement of 4 in passive systems with second-order EPs and significantly larger enhancements (exceeding $400\times$) in gain-aided and higher-order EP systems. In contrast to non-degenerate resonances, EPs are associated with non-Lorentzian emission curves, leading to enhancements that scale polynomially with the resonance quality factor.


Introduction
Electromagnetic resonances enable precise control and enhancement of spontaneous emission and other light-matter interactions. While it is well known that resonances can enhance spontaneous emission rates via the celebrated Purcell effect [1][2][3] by confining light to small volumes for long times, recent work [4][5][6] suggests that giant enhancements can occur via the less familiar Petermann effect [7][8][9][10][11]. The Petermann enhancement factor is a measure of non-orthogonality of the modes in non-Hermitian systems and it appears to diverge when two modes coalesce at an exceptional point (EP)-an exotic degeneracy in which two modes share the same frequency and mode profile [12,13]. In recent years, there has been an explosion of interest in EPs due to the many interesting and counter-intuitive phenomena associated with them, e.g., unidirectional reflection and transmission [14][15][16], topological mode switching [17][18][19], intrinsic single-mode lasing [20,21], and lasers with unconventional pump dependence [22][23][24]. An understanding of spontaneous emission at EPs is essential for their implementation in optical devices, but the existing theory is limited to one-dimensional [25,26] or discrete oscillator systems [27].
In this paper, we present a general theory of spontaneous emission near EPs. Our theory extends beyond spontaneous emission to any light-matter interaction described by the local density of states (LDOS) or, more precisely, any situation in which one analyzes the contribution of a given resonance to the emission of a source, such as narrowband thermal emission [28][29][30][31], absorption [32], perfect coherent absorption [33][34][35], and nonlinear harmonic generation [36]. Whereas traditional theories of spontaneous emission imply infinite enhancement factors at EPs (since the Petermann factor diverges), we use a modified Jordan-form-based perturbation theory to derive (finite) bounds on the enhancement at second-and higher-order EP systems. We show that line narrowing leads to a maximum enhancement of 4 in passive systems with second-order EPs and significantly larger enhancements (exceeding 400×) in gain-aided and higher-order EP systems. Our analytical results are presented in Sec. 2, where we express the emission rate at an EP in terms of the degenerate mode and its corresponding Jordan vector. This derivation assumes negligible dispersion, but we show in appendix B that the effect of dispersion amounts to merely modifying the normalization of the resonant modes, changing the results quantitatively but not qualitatively. Then, in Sec. 3, we demonstrate the implications of our theory via a concrete numerical example of coupled plasmonic resonators. Motivated by the fact that an EP is associated with a double pole in the Green's function, we find specific locations where the emission lineshape becomes a squared Lorentzian, with peak amplitude scaling as Q 2 , where Q is the resonance quality factor (a dimensionless measure of the resonance lifetime). We show that the enhancement at the EP is thus, potentially, much larger than the Purcell factor, which scales linearly with Q. Then, in Sec. 4, we derive bounds on the maximal enhancement at an EP, and we explore these bounds using a periodic system, which allows us to tune gain, loss, and degeneracy independently. Our theory provides a quantitative prescription for achieving large enhancements in practical optical systems, which is applicable to arbitrary geometries and materials and can be implemented with the recent experimental realizations of EPs [16-18, 20-22, 37, 38].
Traditional enhancement formulas fail at EPs since they are based on non-degenerate perturbation theory, which is invalid at EPs. Standard perturbation theory relies on Taylor expansions of differentiable functions while, near EPs, eigenvalues change non-analytically in response to small matrix perturbations. Instead, one needs to use a Jordan-form-based perturbative expansion [39,40]. Although Jordan-vector perturbation theory is well known in linear algebra, along with related results on resolvent operators, these algebraic facts have not previously been applied to analyze Purcell/Petermann enhancement or LDOS in a general EP setting. By using such a modified expansion, we obtain a quantitative formula for the LDOS, which is a measure of how much power a dipole source can radiate [41] or, equivalently, a measure of its spontaneous emission rate. Note that similar expansion methods were previously used to evaluate the Green's functions at EPs [42][43][44]; however, these works were limited to one-dimensional and paraxial systems, and were not applied to study spontaneous emission. An alternative semi-analytic approach was presented in [25]. In this work, the authors applied a scattering matrix formulation to model a simple one-dimensional system and analyzed its emission properties under PT -symmetry conditions. With proper modifications, such an analysis could be generalized to handle more complicated one-dimensional structures (e.g., continuously varying media or complex layered media). However, our general formulas [Eqs. (9) and (10)] can be directly applied to any system with an EP (e.g., three-dimensional photonic or plasmonic structures). Our results demonstrate that the unique spectral properties at EPs are general and do not rely on certain symmetry or dimensionality. Moreover, our theory enables modeling complex experimental apparatuses, performing numerical optimization and design, and deriving bounds on the enhancement, thereby clarifying the usefulness and limitations of EPs for enhancing light matter interactions.
Formally, the Petermann factor is inversely proportional to the "unconjugated norm" of the resonant mode ∫ dx ε E 2 n , where E n is the mode profile and ε is the dielectric permittivity (with modifications to this "norm" when treating dispersive media [45].) At an EP, the unconjugated norm vanishes, ∫ dx ε E 2 n = 0 [13] (a property also called "self-orthogonality" [12]), and the Petermann factor diverges. In fact, the Petermann factor can only diverge at an EP. This is because the Petermann factor is proportional to the sensitivity of an eigenvalue to perturbations [46] (its "condition number"), and that sensitivity can only diverge when two eigenvectors coalesce (i.e., at an EP) by the Bauer-Fike theorem [46]. This implies that our theory is applicable to any system exhibiting a giant Petermann factor. Specifically, in any laser and optical parametric oscillator (OPO) system where a giant Petermann factor was identified [47][48][49], there must have been a nearby "hidden" EP.

Local density of states and Green's function expansions
In the following section, we give some background on LDOS calculations in non-degenerate systems, i.e., systems without EPs (Sec. 2.A), and then we review perturbation theory for systems with EPs (Sec. 2.B). Finally, in Sec. 2.C, we present a condensed derivation of our key analytical result-a formula for the LDOS at an EP [Eqs. (9) and (10)].

2.A. LDOS formula for non-degenerate resonances
The spontaneous emission rate of a dipolar emitter, oriented along the directionê µ , is proportional to the local density of states (LDOS) [50][51][52], which can be related to the dyadic Green's function G via [41,50,53] Here, G is defined as the response field to a point source J = δ(x − x )ê µ at frequency ω. More generally, currents and fields are related via Maxwell's frequency-domain partial differential equation, (∇ × ∇ × − ω 2 ε)E = iωJ, where ε is the dielectric permittivity of the medium. Throughout the paper, we use bold letters for vectors, carets for unit vectors, and Greek letters for vector components. Moreover, we set the speed of light to be one (c = 1). Computationally, one can directly invert Maxwell's equations to find G and evaluate Eq. (1), but this provides little intuitive understanding. A modal expansion of the Green's function, when applied properly, can be more insightful. Away from an EP, one can use the standard modal expansion formula for non-dispersive media [54]: (We review the derivation of this formula for non-dispersive media in appendix A and treat dispersion effects in appendix B). Here, E R n is a solution to the source-free Maxwell's equation with outgoing boundary conditions or, more explicitly, is a right eigenvector of the eigenvalue problem: In reciprocal media ε = ε T , and one can easily derive a simple relation between left and right eigenvectors: E L n = εE R n . Right and left modes which correspond to different eigenvalues are orthogonal under the unconjugated "inner product" [56][57][58]. [The convergence of the denominator (E L n , E R n ) is proven in appendix C.] Due to the outgoing boundary condition, the modes solve a non-Hermitian eigenvalue problem whose eigenvalues (ω 2 n ) are generally complex, with the imaginary part indicating the decay of modal energy in time (in accordance with our intuition that typical resonances have finite lifetimes). From Eq. (2), it follows that the eigenfrequencies, ±ω n , are poles of the Green's function-a key concept in the mathematical analysis of resonances [59]. When considering dispersive media, the denominator in Eq. (2) changes to , as shown in appendix B. In many cases of interest, one can get a fairly accurate approximation for the LDOS by including only low-loss resonances in the Green's function expansion [Eq. (2)] since only those contribute substantially to the emission spectrum. Under this approximation (i.e., considering only resonances ω n = Ω n − iγ n which lie close to the real axis in the complex plane, with γ n Ω n ), the spectral lineshape of the LDOS reduces to a sum of Lorentzian functions, weighted by the local field intensity: Near the resonant frequencies, ω ≈ Ω n , the peak of the LDOS scales linearly with the resonance quality factor Q n ≡ Ω n 2γ n , leading to the celebrated Purcell enhancement factor [1]. On the other hand, the "unconjugated norm," (E L n , E R n ), which appears in the denominator of Eq. (3) leads to the Petermann enhancement factor, defined as [5]. In non-Hermitian systems, the mode profiles (E n ) are complex and the Petermann factor is, generally, greater than one. At the extreme case of an EP, the unconjugated norm in the denominator vanishes and the enhancement factor diverges. However, this divergence does not properly describe LDOS or spontaneous emission at EPs since Eq. (3) is invalid at the EP. That is because the derivation of Eq.
(2) assumes that the set of eigenvectors of the operatorÂ spans the Hilbert space, but this assumption breaks down at the EP. In order to complete the set of eigenvectors ofÂ into a basis and obtain a valid expansion for the Green's function and the LDOS at the EP, we introduce in the following section additional Jordan vectors [39,46,60].

2.B. Jordan vectors and perturbation theory near EPs
At a (second order) EP, the operatorÂ 0 is defective-it does not have a complete basis of eigenvectors and is, therefore, not diagonalizable. However, one can find an eigenvector (E R 0 ) and an associated Jordan vector (J R 0 ), which satisfy the chain relations [39]: where λ EP = ω 2 EP is the degenerate eigenvalue. Equivalent expressions can be written for the left eigenvector E L 0 and Jordan vector J L 0 . In order to uniquely define the Jordan chain vectors, we need to specify two normalization conditions, which we choose to be (E L 0 , J R 0 ) = 1 and (J L 0 , J R 0 ) = 0. Near the EP, on can find a pair of nearly degenerate eigenvectors and eigenvalues that satisfŷ where p 1 represents a small deviation from the EP. [More explicitly,Â(p) ≡ 1 ε(p) ∇ × ∇× = A 0 +Â 1 p+ O(p 2 ), withÂ 0 being defective]. In order to obtain consistent perturbative expansions for E ± and λ ± near the EP, one can use alternating Puiseux series [39]: Substituting Eq. (6) into Eq. (5) and using the additional normalization condition (J L 0 , E R ± ) = 1, one finds that the leading-order terms in the series are where ∆ = In the next section, we use these results to derive a formula for the LDOS at the EP.

2.C. LDOS formula at exceptional points
Near the EP (i.e., for small but non-zero p), one can use the non-degenerate expansion formula Eq. (2) to compute G. In order to compute G at the EP, we substitute the perturbative expansions for λ ± and E ± [Eq. (7)] into Eq. (2) and take the limit of p going to zero, namely: The denominators in Eq. (8) vanish in the limit of p → 0 since (E L ± , E R ± ) = ±2 p 1 2 λ 1 + O(p) [as follows from Eq. (6)]. Most importantly, however, the opposite signs of the denominators lead to cancellation of the divergences and a finite value remains, leading to A formula for the LDOS at the EP is obtained by taking the imaginary part of Eq. (9). Considering again low-loss resonances (at which the enhancement is largest), and evaluating the LDOS near the EP frequency (ω ≈ ω EP ), we find Crucially, Eqs. (9) and (10) yield finite results at the EP resonance frequency Ω EP . Moreover, a squared Lorentzian prefactor, Ω n 2π ( γ n (ω−Ω n ) 2 +γ 2 n ) 2 , replaces the traditional Lorentzian prefactors near non-degenerate resonances, 1 π γ n (ω−Ω n ) 2 +γ 2 n [compare with Eq. (3)]. This unique spectral feature follows directly from the existence of a double pole in the modified expansion formula for G [the first term in Eq. (9)]. As shown in the following section, a squared Lorentzian lineshape implies a narrower emission peak and greater resonant enhancement in comparison with a non-degenerate resonance at the same complex frequency. The strength of our formulation [Eqs. (9) and (10)] is that it is applicable to arbitrary structures and can therefore be used to design and optimize complex three-dimensional geometries with EPs. Moreover, it clarifies the conditions under which LDOS enhancement is maximal; essentially, determined by the relative magnitude of the two terms in Eq. (9), depending on the location of the emitter.

Properties of the LDOS at EPs
In this section, we explore the spectral emission properties at EPs [following Eqs. (9) and (10)] via a numerical model system of coupled plasmonic resonators.

3.A. Example: EPs in plasmonic systems
A convenient approach for tailoring the LDOS in practice is by using plasmonic resonances in metallic structures, which can enable ultra-high LDOS enhancements and are widely used in various applications [61]. In this subsection, we numerically explore a system of two plasmonic resonators with approximate parity-time (PT ) symmetry. PT symmetric systems in optics are characterized by having balanced distributions of gain and loss, and are known to possess EPs at critical gain/loss values at which the mode profiles undergo spontaneous symmetry breaking [62]. Our numerical setup is shown on the left-hand side of Fig. 1(a). It includes two rods with metallic cores (ε = −2.3 + 0.0001i) and a silica coating (ε = 2.1316) surrounded by air (ε = 1). The dimensions of the structure (b = 0.643a, t 1 = 0.536a, t 2 = 0.16a, and t 3 = 0.268a) were tuned in order to have two nearly degenerate resonances that are spectrally separated from the neighboring resonances of the structure. We implement outgoing boundary conditions using perfectly matched layers (PML) [41]. In order to preserve the approximate PT symmetry of the system, gain and loss are added symmetrically to the outer parts of the coating. By brute-force optimization, we find that an EP occurs when the gain/loss value is |Im ε| ≈ 0.0002551 and the background permittivity of the upper half of the silica coating is slightly shifted to ε ≈ 2.129. The right-hand side of Fig. 1(a) depicts the trajectories of the two eigenvalues (red and blue curves) that merge at the EP (orange dot) upon varying the gain/loss parameter |Im ε|. More details on the numerical optimization procedure are given in appendix D.  We discretize Maxwell's equations using a finite-difference frequency-domain (FDFD) approach [63] and compute the LDOS by taking imaginary part of the Green's function. The Green's function is found via three methods: (i) directly inverting Maxwell's partial differential equation, (ii) using the non-degenerate expansion formula, Eq. (2), which is valid away from the EP, and (iii) using the degenerate formula, Eq. (9), at the EP. In principle, one can use the non-degenerate formula, Eq. (2), very close to the EP to compute the LDOS, relying on cancellation between the terms to yield a finite result. However, non-degenerate perturbation theory obfuscates the finite enhancement at the EP, confusing previous estimates of LDOS, whereas our formulation, based on degenerate perturbation theory, naturally captures the finite enhancement.
Our numerical results are shown in Fig. 1(b). The plot on the left shows the Petermann factor, which diverges at the EP, while the right-hand side plot shows the finite LDOS. We excite TM modes and, therefore, compute the LDOS for x-polarized modes (denoted by LDOS x ). The red and blue curves show the two terms that contribute to the sum in Eq. (2). Upon approaching the EP, the individual contributions to the sum diverge with opposite signs, while their sum (black curve) remains finite. The red and blue curves are scaled by 10 −3 for ease of presentation.

3.B. Simplified model for the LDOS
Although our general formula for the LDOS at the EP can be directly applied to the plasmonic system of Fig. 1, it is useful to introduce a simplified model to interpret the results. First, we project Maxwell's operator,Â = ε −1 ∇ × ∇×, onto the two-dimensional subspace spanned by the nearly degenerate eigenvectors near the EP, thus producing a reduced 2 × 2 matrix,Â (which gives a meaningful description of the system as long as the emission spectrum is dominated by the two coalescing resonances). Second, we employ an approach similar in spirit to coupled-mode theory (CMT) [64,65], which involves expressing the modes of the two-rod system (the coupled system) in terms of modes of two reduced systems (the uncoupled systems), containing only one or the other rod. Such an approach is valid whenever the rod-rod separation (t 1 ) in the coupled system is large enough so that the frequency splitting induced by the coupling [κ depicted in Fig. 1(a)] is smaller than the uncoupled resonance frequencies. Denoting by U and V the matrices whose columns are the right and left eigenvectors of the uncoupled system (E R 1,2 and E L 1,2 corresponding to rod 1 and 2 respectively), we find in appendix E that the reduced operator iŝ Here, ω EP ≡ Ω EP − iγ is the degenerate eigenvalue (Ω EP is the resonant frequency while γ denotes the decay rate), and κ is the near-field coupling between the rods. In general, κ is complex; however, in low-loss systems, κ is almost real and, in the current analysis, we neglect its imaginary part entirely for ease of discussion.
is the imaginary frequency shift induced by the gain and loss (Im ε) in the coating. (This definition of η follows from perturbation theory for small gain/loss [55].) This approach is closely related to the recent 2 × 2 formalism used in PT -symmetry works [15,16,[20][21][22]. In fact, the formulations are equivalent for low-loss resonances, γ Ω EP , which is the regime considered in this section. We note that the EP occurs at a complex frequency (i.e., below the real axis in the complex plane) due to outgoing boundary conditions.
Next, we obtain a simplified formula for the LDOS at the EP. The reduced matrixÂ has an EP at the critical gain/loss value: η = κ. Denoting the defective matrix byÂ 0 and the identity opertor by 1 1, we can relate the full Green's function at the EP toÂ 0 via: A formula for the LDOS is obtained by taking the imaginary part of Eq. (12), which allows expressing the LDOS in terms of entries of the 2 × 2 resolvent operator, (Â 0 − ω 2 1 1) −1 , weighted by a product of the left and right local fields E R n E L m (with n, m = 1, 2 denoting the two resonances of the uncoupled system).
The advantage of this formulation becomes apparent when evaluating the LDOS in close proximity to one of the resonators. For instance, near the gain region [e.g., at r 0 in Fig. 1(a)], the lossy mode intensity, ∝ |E R 2 E L 2 |, is negligible, so it follows from Eq. (12) that the LDOS is proportional to the first diagonal entry of the resolvent Moreover, since we consider low-loss resonances, we can normalize the mode profiles so that they are mainly real (Re[E 1 ] ≈ E 1 ), and we find The inset in Fig. 2(b) demonstrates the nearly perfect agreement between this simplified CMTbased LDOS formula (red dashed curve) and brute-force inversion of Maxwell's equation, discretized via FDFD (black curve).

3.C. Quadratic scaling and linewidth narrowing
In this subsection, we apply our CMT-based simplified formulas [Eqs. (12)- (14)] to evaluate the LDOS x near the upper rod [see Fig. 1(a)]. As can be seen from Eq. (14), the LDOS peak value scales as When the resonance-decay rate γ is much smaller than the mode-coupling rate κ, the lineshape approaches a squared Lorentzian curve and M EP scales quadratically with the quality factor Q ≡ Ω EP 2γ [55]. On the other hand, when γ κ, the LDOS peak scales linearly with Q (similar to isolated resonances, as predicted by Purcell [1]). Figure 2(a) demonstrates this key feature. To this end, we computed the LDOS peak, M EP , for several resonance decay rates γ (corresponding to several Q values), varied by introducing homogeneous background gain in the coating. We repeated this procedure for five rod-rod separations (corresponding to five κ values, ranging from 0.0015 to 0.0274), and verified the scaling laws of Eq. (15). The inset in Fig. 2(a) presents the blue and red data points from the main plot on a log-log scale, providing additional confirmation for the quadratic scaling at γ κ (blue) and the linear scaling in the opposite limit of γ κ (red).
Another consequence of a squared Lorentzian lineshape is a narrower emission peak, compared to that of a standard Lorentzian spectrum. The full-width half maximum (FWHM) of a standard Lorentzian curve, 1 π γ n (ω−Ω n ) 2 +γ 2 n , is Γ 0 = 2γ (where Ω 0 ± Γ 0 /2 is the frequency at which the Lorentzian drops to half of its maximal value). On the other hand, following Eq. (14), the FWHM of the LDOS near an EP with Im[ω EP ] = γ is As shown in Fig. 2(b), the FWHM at the EP is always smaller than Γ 0 , approaching a value of √ 5 − 2Γ 0 ≈ 0.48Γ 0 in the limit of κ/γ → ∞. We computed the FWHM directly from the FDFD-discretization of Maxwell's equations (dots) and by using the CMT-based simplified expression, Eq. (14) (black solid line), proving very good agreement between the two methods.

LDOS enhancement at EPs
In the previous section, we showed that a squared Lorentzian lineshape can lead to enhanced emission rates and reduced linewidth. In this section, we quantify the enhancement at the EP and study its bounds. To demonstrate our results, we employ another numerical example: a periodic waveguide (Sec. 4.A). This structure allows us to independently tune gain/loss and degeneracy and, therefore, demonstrate the impacts of the two effects separately. We treat both passive (Sec. 4.B) and active (Sec. 4.D) systems (i.e., systems without and with gain respectively), and then generalize our results (in Sec. 4.C) to higher-order EPs, which form at the coalescence of multiple resonance.

4.A. Example: EPs in periodic structures
To demonstrate our results in this section, we numerically explore the periodic structure shown in Fig. 3(a). The modes of a periodic system are Bloch wavefunctions of the form E(r) = E k (r)e ik·r , where E k (r) is a periodic function and k is the wavevector [66]. At each k, the mode E k (r) solves an eigenvalue problem of the form [64,65]:Â(k)E nk = ω 2 nk E nk , wherê A(k) ≡ ε −1 (∇ + ik) × (∇ + ik)× and n labels the band. The figure of merit for spontaneous emission in periodic structures is the LDOS per wavevector k and field component µ, which is a measure of the power expended by a Bloch-periodic dipole source with a particular k-vector, in the presence of an electromagnetic field polarized along direction µ. We abbreviate it as LDOS k (also called the mutual density of states [67]). The LDOS k can be integrated over k to obtain the LDOS from an isolated (non-periodic) point source in the periodic structure [68]. However, the effect of the EP is much clearer in the integrand (LDOS k ) than in the integral (LDOS), and so we focus here on the former for illustration purposes, exploiting the fact that k allows us to control how close we are to an EP without altering losses.
Our example system consists of a waveguide with periodic index modulation along its central axis,x [ Fig. 3(a)]. PML are used to truncate the transverse (ŷ) direction. The design parameters are: ε 1 = 12, ε 2 = 13.137, d = 0.51a, and t = 0.25a. These parameters were chosen so that the corresponding one-dimensional system (with the same ε 1 , ε 2 , d and infinite thickness, t → ∞) has nearly degenerate second and third frequency bands near k x = 0 [guaranteed by choosing parameters close to the quarter-wave plate condition: We force an EP for TM-polarized modes (E z , H x , H y ) by fine-tuning the wavevector k x and the permittivity contrast δε ≡ ε 2 − ε 1 (see appendix D, Fig. 6). Figure 3(a) also depicts the real and imaginary parts of the coalescing eigenvalues (blue and red curves respectively). At k = 0 (also called the Γ point), Maxwell's eigenvalue problem is complex-symmetric and, consequently, the eigenvectors E R 1 and E R 2 are orthogonal (see mode profiles in the lower-left corner, having even/odd symmetry when x is flipped to −x). At the EP, the eigenmodes E R 1,2 merge into a single degenerate mode E R 0 (upper-right panel). We also compute the Jordan vector (J R 0 ) with our novel iterative method [69], and use it to verify our formula for the LDOS at the EP [Eq. (10)]. Figure 3(b) depicts the LDOS k near the periodic waveguide at several wavevectors (k EP , k 1 , k 2 , k ∞ marked on the lower plot). Far from the EP (at k ∞ ), the LDOS k is a sum of two non-overlapping Lorentzian curves. As k approaches the EP (e.g., at k 1 and k 2 ), the resonance peaks begin to overlap and the LDOS k peaks increase due to the growing Petermann factor. Most importantly, for k values near but not equal to k EP , the traditional modal expansion formula [Eq. (2)] approaches the limiting Jordan-form-based formula [Eq. (9)]. Physically, this means that structures with nearby EPs can be approximated by truly defective structures, making this analysis useful for experimental systems, which can only be close to but not exactly at an EP due to fabrication and calibration imperfections. Computationally, this implies that as long as Maxwell's operator A(k) is not exactly defective, one can use Eq. (2) to evaluate the LDOS k , when properly canceling the divergent terms in the sum. The lower plot in Fig. 3(b) compares the enhancement, M k /M ∞ , for structures with varying quality factors, where we introduce the definitions M k ≡ γ · max ω [LDOS k (ω)] and M ∞ ≡ γ · max ω [LDOS ∞ (ω)]. We change Q by modifying the permittivity contrast δε: Radiation losses decrease with decreasing index contrast of the periodic modulation [70], with the limit of zero contrast corresponding to infinite Q. By plotting the normalized LDOS k peak (M k /M ∞ ) vs. deviation from the EP (∆k ≡ k − k EP ), we find that the enhancement at the EP (M EP /M ∞ ) is always four-fold, regardless of Q. This follows from a sum rule which states that the spectrally integrated LDOS (and, therefore, also the LDOS k ) is a constant [71]. It implies that the maximum LDOS at an EP [i.e., the peak of a square Lorentzian M EP γ 3 [(ω−Ω EP ) 2 +γ 2 ] 2 ] is four times larger than the maximum LDOS at a non-degenerate resonance [i.e., the peak of a simple Lorentzian

4.C. Passive structures with higher order EPs
Motivated by recent interest in higher-order EPs [72][73][74][75], we generalize our results from the previous section to nth order degeneracies (i.e., EPs that form when n degenerate eigenvectors merge). In this case, we define the enhancement factor as the ratio of the LDOS peak at the EP and at a reference point with n non-degenerate resonances (generalizing our earlier definition for M EP /M ∞ ). Following the arguments from Sec. 3.B, we expect to find a squared Lorentzian emission curve at second-order EPs, a cubic Lorentzian curve at third-order EPs, and a Lorentzian to the nth power, L n (ω) = M n γ 2n−1 [(ω−Ω EP ) 2 +γ 2 ] n , at nth order EPs. (This is essentially equivalent to a known result on the rate of divergence of the norm of the resolvent matrix as an nth order EP is approached [46].) From the sum rule mentioned above [71], the spectrally integrated LDOS at an nth order EP, , is equal to the integrated LDOS before the n resonances merge (here, Γ[n] is the gamma function). Realizing that the enhancement at the EP is maximal when merging n identical resonances, a bound can be computed by solving S n = nS 1 . We find that the enhancement at the EP is at most M n /M 1 =
These results imply that higher order EPs could potentially provide a new route for achieving order-of-magnitude enhancement of the LDOS and order-of-magnitude narrowing of the emission linewidth (in contrast to traditional methods, which typically aim to maximize the traditional Purcell factor by increasing the quality factor and reducing the mode volume [1]). However, this result does not necessarily mean that higher-order EPs will yield a larger LDOS than the best lower-order EP or single resonance. The reason is that the degrees of freedom that one would use to bring n resonances together might otherwise be employed to enhance the Q of an individual resonance. In practice, there will be a tradeoff between the quality of individual resonances and the order of the EP.

4.D. Active structures
We show in this section that much greater enhancements can be achieved in active systems, i.e., by introducing gain. Figure 4 compares four periodic waveguides with different index contrasts (δε), corresponding to different passive quality factors, Q p . Gain and loss are added to each of the waveguides in order to force EPs, all of which share the same active quality factor Q a . The empty (filled) markers in Fig. 4(a) are the EP resonances (ω EP ) in the complex plane before (after) adding gain/loss. As shown in Fig. 4(b), the structure with the smallest passive quality factor Q p has the largest LDOS k enhancement at the EP since it requires more gain in order to attain the same Q a . Structures with smaller initial Q p values can lead to even greater enhancements, and we present such a case in appendix F. In principle, the relative enhancement at the EP, M EP /M ∞ , is not bounded in this computational model. However, in practice, giant relative enhancements do not necessarily imply giant absolute LDOS values (since the value of the LDOS at the reference point may be very small). Moreover, when adding enough gain to bring the system to the lasing threshold, stimulated emission eventually limits line narrowing and LDOS enhancement at the EP.
Similar to the analysis of the plasmonic example in Sec. 3, we can introduce a simplified 2 × 2 model to interpret the results [analogous to Eq. (11)]. We project the full Maxwell's operator A onto the subspace spanned by the two modes at k x = 0 [shown in the left-lower corner of Fig. 3

(a)] and we obtainÂ
where v g ≡ ∂ω ∂k k=0 is the the group velocity (similar to the model used in [37]). Prior to adding gain/loss, an EP occurs when v g k = γ p , and the degenerate frequency is ω EP = Ω EP − iγ p . After adding gain/loss, ±iη, the EP frequency moves vertically in the complex plane, becoming Ω EP − iγ a , while the EP condition remains the same: v g k = γ p . In analogy to Eq. (15), this model implies that the maximal enhancement at the EP scales as M EP ∝ γ p γ 2 a + 1 γ a . The quadratic scaling of M EP with Q a in the limit of γ a γ p is demonstrated in appendix F [but this result is essentially another demonstration of the spectral property shown in Fig. 2(b)].

Discussion
The theory presented in this work provides a quantitative prescription for achieving large spontaneous emission rates using EPs, potentially exceeding by orders of magnitude those attained with standard non-degenerate resonances. Such enhancements could be useful for various applications, including fluorescent and Raman sensing [77], high-power low-coherence light sources [78], and sources with tunable coherence [79]. Moreover, by extending the current linear theory to account for EPs in nonlinear systems, our approach could be applied to study the properties of lasers at EPs-a topic that has recently drawn great attention in the optics community [20][21][22]80].
Our formulation for the LDOS at the EP [Eqs. (9) and (10)] establishes that the enhancement generally consists of two terms, one that scales linearly and one that scales quadratically with the quality factor. In Secs. 3,4, we verified this scaling argument via two numerical examples, and we employed simplified 2 × 2 models to estimate the coefficients of the quadratic terms in the LDOS formula. (In the plasmonic system, we found that the coefficient was the coupling constant κ, while in the periodic example, the coefficient was the passive decay rate γ p .) More generally, we show in appendix G that for arbitrary low-loss systems, the quadratic coefficient is bounded by This result provides an easy-to-evaluate upper bound on the maximal enhancement in complicated geometries, depending only on the maximal gain/loss of the constituent materials. More explicitly, in active systems, we find that the maximal enhancement at the EP is bounded by M EP M ∞ = 2(1 + κ Im[ω EP ] ) ≤ 2 1 + 2Q| max |Im ε| . Note that in our plasmonic example, the quadratic coefficient is within 10% of the bound.
Finally, our theory extends beyond spontaneous emission, and can be applied to a broader class of phenomena described by the LDOS. For example, we anticipate similar enhancements in higher-harmonic generation rates in nonlinear media (e.g., Kerr media). In that case, the input lower-harmonic field (multiplied by the nonlinear susceptibility) acts as a source to the higher-harmonic field. To lowest order in the nonlinearity, the converted power is found by convolving the Green's function with the input signal. This is analogous to our formulation of spontaneous emission (which involved convolving a dipole source with the Green's function), and would therefore result in similar enhancements [81]. More generally, a similar treatment could produce enhancements in related quantities at EPs in other areas of physics (e.g., exciton-polariton [38] mechanical systems [18], and also leak-wave antennas both at radio-frequency and visible frequencies [82,83]). Finally, our theory can be generalized to study scattering and extinction problems. We find, however, that even though EPs can produce special spectral features in scattering cross-sections, they do not give rise to giant scattering enhancements, since the scattered intensity is bounded by the incoming intensity (and the bound of perfect scattering can be easily achieved with a single non-degenerate resonance [84]).

A. Non-degenerate Green's function modal expansion formula
In this section we review the derivation of the Eq. (2) in the main text. Our derivation is similar to standard methods for Hermitian eigenvalue problems [54], with the necessary modifications for treating general non-Hermitian systems (most importantly, using the unconjugated "inner product" [12] between left and right eigenmodes).
Given a non-magnetic medium with dielectric permittivity ε, the fields E and currents J in the medium are related via Maxwell's frequency-domain partial differential equation: whereÂ r ≡ 1 ε ∇ r × ∇ r ×. The response to arbitrary currents can be found by convolving the dyadic Green's function with the current sources: In this section, we expand G in terms of right and left resonant modes (E R n and E L n ), which are outgoing solutions of the partial differential equations When the set of eigenvectors of Eq. (21) forms a complete basis of the Hilbert space, one can introduce the completeness relation, which consists of expanding the Dirac delta function as where ⊗ is the outer/tensor product u ⊗ v = uv T . The question of completeness of eigenmodes in non-Hermitian open systems has not been proven in general, for arbitrary three-dimensional systems. However, since in this work we always evaluate the Green's function in close proximity to the resonators and at frequencies close to the resonance frequencies, the eigenmodes which overlap spectrally and spatially with the emitter give a good approximation for the LDOS, justifying the use of Eq. (22). Similarly, we wish to find an expansion formula for G or, more explicitly, find the coefficients a n (r ) in the series G(r, r ) = n E R n (r) ⊗ a n (r ).
To this end, we substitute Eq. (22) and Eq. (23) into Eq. (20) and obtain n A r − ω 2 E R n (r) ⊗ a n (r ) = − n E R n (r) ⊗ E L n (r ).
Next, we multiply both sides of the equation by [E L m (r)] T . Using the relation: m [E L m (r)] T , integrating over r, and invoking the bi-orthogonality relation for non-Hermitian systems [85]: we obtain a m (r ) = Finally, substituting this result in Eq. (23) we obtain an eigenmode expansion of the dyadic Green's function [which reduces to Eq. (2) in the text]: The integral in the denominator of Eq. (26) is one, but we keep it here for comparison with Eq. (2). Last, we note that in reciprocal media ε = ε T , and there exists a simple relation between left and right eigenvectors: E L n = εE R n . More generally, the left and right eigenvectors of a symmetric generalized eigenvalue problem (EVP): To see this, rewrite the EVP for the left vectors as which shows that B −1 E L n is a right eigenvector. Note that this relation holds also with Bloch-periodic boundary conditions (and the surface-term correction found in [86] does not appear in our formulation). Although the matrix A is no longer symmetric in the k-periodic problem, it satisfies the relation A(k) T = A(−k), and since we are essentially relating k-right eigenvectors to (−k)-left eigenvectors, the relation above remains unchanged.

B. The effect of dispersion on the LDOS formula
In this appendix, we consider the effects of dispersion on the Green's function near and at the EP. In accordance with previous work on quasi-normal modes in dispersive media [87,88], we find that the Green's function has non-diagonal contributions ∝ E L ± ⊗ E R ∓ near the EP. However, exactly at the EP, the Green's function has exactly the same form as the non-dispersive degenerate formula [Eq. (9)], with dispersion affecting only the normalization of the degenerate mode (E 0 ) and Jordan vector (J 0 ).
In the same spirit as our derivation of the non-dispersive formula, we expand the Green's function G in eigenmodes: and use simple algebraic manipulations to express the coefficients α m in terms of the modes.
Recall that G is defined as the solution to the partial differential equation: By multiplying both sides of Eq. (28) from the left by E L n (r) and integrating over dr, we obtain Then, by substituting Eq. (27) into Eq. (29), we find (Note that this result was also derived in [88] by invoking Lorentz reciprocity.) Since we are interested in emission from emitters in close proximity to the resonators and in frequencies near the resonant frequencies, we may assume that a finite set of N eigenvectors adequately describes the system's response [similar to our assumption in the non-dispersive derivation in Eq. (22)]. We introduce the vector s(r) = {E 1 (r), . . . , E N (r)} and the matrix O mn ≡ (E m , E n ), and rewrite Eq. (30) as: s(r 0 ) = O(ω)α(ω, r 0 ) or equivalently α(ω, r 0 ) = O −1 (ω)s(r 0 ). With this notation, Eq. (27) can be rewritten as: where Equation (31) . Next, let us calculate the limit of Eq. (31) when two modes E ± coalesce at an EP. Keeping only the terms corresponding to E ± in Eq. (31), expanding both ω and ω ± in Taylor series around ω 0 and introducing the notation:

and the matrix inverse is
Substituting Eq. (34) into Eq. (31), we find that the Green's function near the EP is Last, we want to take the limit of Eq. (35) as the two modes E ± approach the EP. We expand ω ± and E ± around the EP using degenerate perturbation theory (using the notation of Sec. 2.B): Generalizing our derivation of the non-dispersive formula, we choose the normalization conditions The first two terms in Eq. (35) have the same form as the non-dispersive expansion [Eq. (8)]: . Upon approaching the EP, each contributions diverges with an opposite sign, and a finite contribution remains. The last two terms scale as √ p near the EP, and vanish at the EP. Finally, we obtain Equation (37) implies that the green's function at the EP in dispersive media has the same form as the non-dispersive formula [Eq. (9) in the text], but the normalization of the modes changes.

C. Convergence of the "unconjugated norm"
In this appendix, we show that when perfectly matched layers (PML) are used to implement outgoing boundary conditions, the "unconjugated norm" of a scattering eigenmode, ∫ dx ε(x)E n (x) 2 , converges to a finite result as the PML thickness tends to infinity. Scattering eigenmodes are solutions to Maxwell's eigenvalue problem, ε −1 ∇ × ∇ × E n = ω 2 n E n [55], with outgoing radiation conditions. These solutions (also called "leaky modes" [89]) have complex frequencies ω n which lie in the lower-half of the complex plane (Im[ω n ] < 0) [90] and, consequently, the modal amplitudes (|E n |) grow unboundedly at large distances from the structure. Even though the modal amplitude diverges, we show that the unconjugated inner product is finite and independent of the PML parameters, as long as the PML works effectively (i.e., designed so that outgoing waves are attenuated exponentially inside the PML and are not reflected at the air-PML interface.) We provide a simple proof for a one dimensional geometry. An abstract proof for the independence of the unconjugated norm on the PML parameters appeared in [88] (using a generalized definition of the norm, which includes dispersion). An alternative one-dimensional proof was given in [91], where the authors used analytic continuation of the coordinates, similar to the coordinate stretching method that we use here.
We consider the one-dimensional geometry depicted in Fig. 5(a), which consists of a slab of thickness a and refractive index n, embedded in air. (Generalizations to more complex one-dimensional structures straightforwardly follow.) We truncate the computational cell by placing PML at a finite distance from the slab (|x| = L/2), and we impose metallic boundary conditions at the cell boundary (x = ± N 2 ). The PML can be viewed as an analytic continuation of the spatial coordinate into the complex plane [92] (The latter condition guaranties that fields oscillating at different frequencies ω will be attenuated inside the PML at rate that is independent of ω.) For concreteness, we choose: is the PML thickness. Consequently, the coordinate stretching transformation inside the PML is Scattering solution for this geometry can be written explicitly as A n e −iω n x−(σ 0 /2d)(x+L/2) 2 + B n e iω n x+(σ 0 /2d)(x+L/2) 2 −N/2 < x < −L/2 A n e −iω n x + B n e iω n x −L/2 < x < −a/2 e −iω n nx + e iω n nx −a/2 < x < a/2 A n e iω n x + B n e −iω n x a/2 < x < L/2 A n e iω n x−(σ 0 /2d)(x−L/2) 2 + B n e −iω n x+(σ 0 /2d)(x−L/2) 2 L/2 < x < N/2 , where the coefficients A n , B n and resonant frequency ω n are found by requiring continuity of the field (E n ) and its derivative (dE n /dx) at the interfaces, while imposing the boundary condition: E n = 0 at ±N/2. We obtain where we introduced C n ≡ e iω n N −σ 0 d .
Next, we compute the unconjugated norm, ∫ εE 2 n , and study its convergence. Introducing the antiderivative function the unconjugated norm can be written as We now show that whenever the condition Eq. (41) (below) holds, I is finite and independent of the PML parameters in the limit of d → ∞. Since the coordinate stretching factor is zero at the air-PML boundary [ x( L 2 ) = L 2 ], the antiderivative terms at L/2 cancel (J L 2 = J x( L 2 ) ). Next, consider the boundary term, J x( N 2 ) . Straightforward integration of Eq. (39) yields three terms, all of which decay exponentially as d → ∞, provided that |C n | decays exponentially.
[The first is Introducing ω n = ω n − iω n and α = d/L, one finds from the definition of C n that it decays exponentially whenever Evaluating the remaining terms in Eq.    Fig. 3(a). The EP is found by numerically minimizing the distance between the two eigenvalue sheets (red and blue surfaces). Cyan dots: Eigenvalues for fixed k x = k EP and varying δε. Black dots: Eigenvalues for fixed δε = δε EP and varying k x .

D. Computational details on Forcing EPs in periodic waveguides
In order to force an EP (i.e., an accidental degeneracy of two nearby resonances), we need to satisfy a single complex equation (ω m = ω n ), which can be done by searching two real parameters. In the plasmonic example in Sec. 3, we search the two-dimensional parameter space spanned by the gain/loss parameter (|Im ε|) and the (real part) of the refraction index in the upper-half of the silica coating. In the periodic example from Sec. 4, we search the two-dimensional space spanned by the wavevector k x and the permittivity contrast δε. Figure 6 presents our numerical results for the periodic example. As shown, we find a degenerate resonance at ω EP = 0.3924377 − 0.00004119303i when k EP a/2π ≈ 1.468 × 10 −4 and δε EP ≈ 1.137.  Figure 7 shows our numerical results. At each wavevector, we compute the LDOS k peak value (M k ) vs. deviation from the EP (∆k ≡ k − k EP ) when adding increasing amounts of gain to the waveguide (i.e., fixing Q p while increasing Q a ). The LDOS is evaluated at r 0 [ Fig. 3(a)] by direct inversion of Maxwell's equations. We show here enhancements of ≈ 400, while higher values can easily be obtained; the enhancement is essentially not bounded in this computational model. However, in reality, it is bounded by quantum noise near threshold [93]. The side panels in the figure demonstrate that the LDOS peak scales quadratically with Q a near the EP [ Fig. 7(a)] and linearly with Q a away from the EP [Fig. 7(c)]. Finally, we note that when evaluating the LDOS near the center of the computational cell (i.e., at x ≈ 0), the lineshape changes dramatically, and it actually has a minimum at the resonance frequency and two side peaks, whose amplitude scales as Q 2 a (not shown).

G. Theoretical limit of LDOS enhancement at an EP
In this section, we derive Eq. (18) in the text. Let us first define an effective mode amplitude E, E ≡ ∫ C dx |E| 2 , where C denotes a finite region containing the geometry (e.g., the last scattering surface). In order to obtain an upper bound on the LDOS enhancement at the EP, we need to estimate the quantities E R 0 , E R 0 / J R 0 , J R 0 and E L 0 , E L 0 / J L 0 , J L 0 , which determine the relative magnitude of the two terms in the expansion formula for G at the EP [Eq. (9)]. Let us decompose the complex-symmetric Maxwell's operator into:Â =Â + iÂ , whereÂ ≡Â +Â * 2 andÂ ≡Â −Â * 2i , and the asterisk denotes complex conjugation. In many cases of interest, one can assume Â Â (under an appropriate matrix norm). In such cases, one can use defective perturbation theory (Sec. 2.B) to expand the eigenmodes E R ± and eigenvalues ω ± ofÂ in terms E R 0 , J R 0 and λ EP . Using Eq. (6), we obtain SinceÂ is a real Hermitian operator, it has real and orthogonal eigenvectors E ± . Using Eqs. (47) and (48) and assuming E + , E + ≈ E − , E − and E − , E + = 0, one obtains Substituting the explicit perturbative expansion Eq. (7), we find (with an equivalent expression for the left eigenvector). Next, recall the definition ofÂ ≡ 1 ε ∇×∇×. In the limit of low losses, Im ε/ε 1, we can approximateÂ ≈ − Im ε ε 2 ∇ × ∇× and, consequently where, in going from the first to the second line, we used the relation E R 0 = ε −1 E L 0 (which holds for Maxwell's eigenvalue problem) and, in the following approximations, we used the property E 2 0 ≈ |E 0 | 2 (valid for low-loss systems) and the normalization condition ∫ |E 0 | 2 = 1. This completes the proof of Eq. (18).