Quasinormal modes and Purcell factors of coupled loss-gain resonators and index-modulated ring resonators near exceptional points

We first present a quasinormal mode (QNM) theory for coupled loss-gain resonators working near an exceptional point. Assuming linear media, which can be fully quantified using the complex pole properties of the QNMs, we show how the QNMs yield a quantitatively good model to a full dipole spontaneous emission response in Maxwell's equations at various spatial positions and frequencies (linear response). We also develop a highly accurate and intuitive QNM coupled-mode theory, which can be used to rigorously model such systems using only the QNMs of the bare resonators, where the hybrid QNMs of the complete system are automatically obtained. Near a lossy exceptional point, we analytically show how the QNMs yield a Lorentzian-like and a Lorentzian-squared-like response for the spontaneous emission lineshape, consistent with other works. However, using rigorous analytical and numerical solutions for microdisk resonators, we demonstrate that the general lineshapes are far richer than what has been previously predicted. Indeed, the classical picture of spontaneous emission can take on a wide range of positive and negative Purcell factors from the hybrid modes of the coupled loss-gain system. These negative Purcell factors are unphysical and signal a clear breakdown of the classical dipole picture of spontaneous emission in such media, though the negative local density of states is correct. We also show the rich spectral features of the Green function propagators, which can be used to model various physical observables. Second, we present a QNM approach to model index modulated ring resonators working near an exceptional point and show unusual chiral power flow from linearly polarized emitters, in agreement with recent experiments, which is quantitatively explained without invoking the interpretation of a missing dimension (the Jordan vector) and a decoupling from the cavity eigenmodes.


I. INTRODUCTION
Lossless photonic systems (such as closed resonators with no material absorption) can be formulated as a Hermitian eigenvalue problem, which yields real eigenfrequencies from the source-free Helmholz equation, and corresponding normal modes (NMs). This is also true for periodic systems with bound modes, e.g., lossless waveguide modes. However, real cavity structures (resonators) with open boundary conditions yield finite loss (or gain) and produce complex eigenfrequencies. Thus most optical systems are naturally dissipative via material absorption and/or radiation.
To have a detailed understanding of such systems, and also for connecting to new applications in quantum optics, it is desirable to have an accurate model of coupled loss-gain resonators at the level of a rigorous and intuitive mode theory, which can allow one to describe lightmatter interactions at various spatial points and frequencies. From a theoretical perspective, temporal coupled mode theory (CMT) has proved to be an efficient approach to investigate coupled resonator systems [54,55], where only the solution from the bare systems are required, and one assumes that the coupled modes can be represented by a linear superposition of the bare modes (weak coupling regime). However, the coupling coefficients are typically used as heuristic parameters in that they are usually extracted from fitting the full solution of the coupled system properties [37,[56][57][58], or they are mainly used to explain the basic physics of coupling.
Another potential problem with such approaches is arXiv:2101.07633v2 [physics.optics] 24 Mar 2021 that the underlying modes of the bare resonators are assumed to be NMs (Hermitian system), and a finite decay rate, to account for real losses, is added phenomenologically. This finite decay rate already relates to a non-Hermitian problem, so it is natural to model the bare resonators also with a non-Hermitian theory, where the correct eigenvalues and modes would then be obtained for a more general CMT for loss-gain systems. This is not only a more correct approach, but it changes some of the fundamental coupling regimes and constraints significantly, and open up much richer light-matter interaction regimes.
In this work, we first present a QNM theory for general media containing both lossy resonators and gain resonators. We describe a rigorous and intuitive CMT based on the Green function solution for the coupled loss-gain QNMs, where we analytically obtain the hybrid modes from only the QNMs of the bare resonators (gain or loss). We demonstrate the extremely high accuracy of the analytical theory by comparing with full numerical dipole simulations for whispering-gallery modes (WGMs) of microdisk resonators, and show excellent agreement for various designs and spatial positions without using any fitting parameters. A QNM approach also allows one to justify the the underlying assumptions of a linear gain medium, which requires an analysis of the complex poles. Without such an analysis, then composite systems may not even constitute a physically meaningful solution with gain treated at the level of a material response in Maxwell's equations. Indeed, as far as we are aware, this is likely the only approach to confirm this directly (in the absence of having an analytical expression for the medium Green function).
After obtaining the physically meaningful QNMs for loss-gain resonator systems, we apply our theory to study the unusual Purcell factors and Green function propagators at various spatial positions for eigenfrequencies close to an EP. For our numerical example, the considered coupled loss-gain disk resonators are shown in Fig. 1(a), where one microdisk has material loss and the other has material gain. Such systems with balanced gain and loss form a general optical system to investigate PT symmetry and EP physics [37,38,45]. In our model, we choose a gain coefficient that is slightly less than the loss coefficient, which can support hybridized QNMs with finite Figure 1. Schematic diagram of a coupled loss and gain microdisk resonator system. The radius of both resonators is R = 5 µm. The refractive index of the lossy (gain) resonator is n loss (ngain), in a background medium with nB = 1 (free space). The gap distance between the resonators is dgap, and the dipole (emitter) is placed within the gap, at several possible positions: r1, r2, r3, r4, and r5 (see text). The origin of the coordinate system is at the gap center.
loss, which is also a requirement for assuming linear gain media [1].
Next, we apply our general QNM approach to model emitters coupled to index-modulated ring resonators, which were recently experimentally studied [81]; in this work, the authors found the following interesting observation, quoting from their abstract: "We find a chiralityreversal phenomenon in a ring cavity where the radiation field reveals the missing dimension of the Hilbert space, known as the Jordan vector. This phenomenon demonstrates that the radiation field of an emitter can become fully decoupled from the eigenstates of its environment." In contrast to this interpretation, using similar resonator designs, we show that the chiral reversal behaviour is quantitatively explained by the emitter coupling to the two fundamental QNMs of the resonator structure. We explain this net chiral power flow analytically in terms of only the QNM properties, and show excellent agreement with full dipole simulations from Maxwell's equations.
The rest of our paper is organized as follows: In Sec. II, we introduce optical QNMs, Green functions in terms of QNMs, and show how these relate to SE decay and Purcell factors. In Sec. III, we present a detailed CMT, using both NMs and QNMs, to obtain analytical insight into coupled loss-gain resonators. We subsequently use these to explain when an EP may form, and how things change when one uses a QNM theory. We obtain explicit expres-sions for the hybrid modes using only the QNMs of the bare loss or gain resonators. We explain the limits and failures of the usual CMT for such systems. Section IV discusses Green functions and Purcell factors at the EP, and shows how a Lorentzian-like or Lorentzian-squaredlike lineshape forms [82][83][84]. Section V presents detailed numerical results for coupled microdisk resonators, and confirms the excellent agreement with our analytical CMT and full dipole solutions, for various gap distances between the resonators. We then study various Purcell factor regimes at different dipole positions as a function of frequency, and show highly unusual and rich spectral lineshapes, including negative Purcell factors, and discuss the essential role of the QNM phase; again, all of these show quantitatively good agreement with full dipole calculations. Negative total purcell factors are not physical and motivate the need for a corrected derivation of the accepted Fermi's golden rule for such media, which is described elsewhere [85]. We also show the Green function propagators, which connects to various observables outside the resonators, which also yield rich non-Lorentzian lineshapes. Finally, in Sec. VI, we study EP physics for emitters coupled to index-modulated ring resonators, and connect to recent experiments on chiral emission for linear dipoles [81] using the QNM theory. We give our conclusions in Sec. VII.
In addition to the main text, we also present several Appendices. Appendix A discusses the numerical QNM normalization approaches using COMSOL, where we show three different approaches yielding the same normalized QNMs within numerical precision. Appendix B discusses why we only need to consider one QNM of the miscrodisk resonator for the dipole locations we study, which is constructed from a symmetric linear combination of clockwise and counter clockwise WGMs. Full dipole calculations in COMSOL are also discussed in Appendix C, which are used to check the validity of the QNM results. To compare with the coupled cavity systems, the results for single loss and single gain cavities are shown in Appendix D, which also confirms the extremely high accuracy of the single QNM approximation for these resonators. Naturally, one can also solve the coupled system with a QNM approach directly, instead of using CMT from the bare solutions; thus Appendix E shows the direct QNM approach for a coupled resonators, where quantitatively good agreement with our analytical CMT results and full dipole results are obtained. In addition to the loss-gain cavities shown in the main text, we also show two more loss-gain examples in Appendix F, with different gain coefficients.

II. QUASINORMAL MODES AND SEMICLASSICAL THEORY OF SPONTANEOUS EMISSION AND PURCELL FACTORS
We first introduce the electric-field QNMs,f µ (r), which are solutions to the Helmholtz equation, where c is the vacuum speed of light,ω µ = ω µ − iγ µ is the complex eigenfrequency of each QNM, and (r,ω µ ) is the dielectric function, which is in general complex and dispersive, though for our numerical examples below, we will assume this is a constant complex value in the frequency regime of interest (this is not a model restriction). The open boundary conditions ensure the Silver-Müller radiation condition [86]. It is also worth noting that this boundary condition leads to quite different asymptotic behaviour of gain QNMs and loss QNMs due to the change of sign for γ µ . Namely, the lossy QNMs diverge in space but converge in time, while the gain QNMs converge in space, but diverge in time. For the composite system, then the hybrid modes must converge in time, forcing the complex poles to have loss. These subtleties are clearly missing and overlooked in heuristic theories of coupled loss and gain resonators, but are essential to get correct for a physically meaningful model. Using n to represent the refractive index and α the loss or gain coefficient, for a lossy resonator with permittivity L = (n L +iα L ) 2 , we assume a dominant QNM resonancẽ ω L ≈ ω L − iγ L , where α L , γ L > 0. Similarly, for the gain resonator, G = (n G + iα G ) 2 , we have a dominant QNM resonanceω G ≈ ω L − iγ G , where now α G , γ G < 0. The quality factor is defined from Q µ = ω µ /(2|γ µ |). Since we treat the gain amplifier in terms of a linear amplifying medium (e.g., we neglect saturation effects in the gain medium), the composite system must have γ +/− µ >0 for the hybrid modes [1]. We will refer to the gain and loss QNMs asf G andf L , respectively, and the hybrid modes (i.e., in the presence of coupling) asf ± .
To connect to a general definition of the SE in an arbitrary medium, we seek to obtain the Green function, defined through with corresponding radiation conditions, and 1 is the unit tensor. The normalized QNMs can be used to define the Green function for locations near (or within) the scattering geometry through [60,87] with A µ (ω) = ω/2(ω µ − ω). Expanding the Green function with the QNMs can easily be used to compute the SE rate and the Purcell factor. We stress again that the medium must meet the condition for linear amplifying media, which coincides with a causal Green function in the sense of linear response theory (Kramers-Krönig relations). Here, the Green function must be analytic in the upper half complex plane to fulfill the Kramers-Krönig relations, and this can be rigorously justified by using a QNM approach. Indeed, without such an approach, it is not known whether the model represents a physically meaningful solution for Maxwell's equations. Considering a dipole emitter, d = d 0 n d , at location r d , then the classical SE rate is [66] and the generalized Purcell factor reads [66,88] where Γ 0 (r d , ω) = 2d · Im{G B (r d , r d , ω)} · d/( 0 ), and G B is the Green function for a homogeneous medium (known analytically). For a 2D TM (TE) dipole, Im{G B (r d , r d , ω)} = ω 2 /4c 2 (ω 2 /8c 2 ). The factor of 1 appears naturally for dipole positions outside the resonator [89].
For an arbitrary photonic cavity medium, the QNMs for both the bare resonators (i.e., without coupling), and also for the coupled system, can be obtained from an efficient dipole scattering approach in complex frequency [67], described in more detail in Appendix A. The total Green function can also be obtained numerically from the full dipole response (i.e., without any modal approximations), which we carry out in COMSOL to check the accuracy of the QNM expansion form. Although the hybrid QNMs can be obtained numerically as well, it is far more insightful to develop a coupled-mode formalism to describe the coupling geometry.

A. Wave equation and normal modes
Before developing a QNM CMT, here we present a NM approach and also connect to the common literature for describing when EPs can occur for coupled loss-gain cavity modes.
To simplify the equations and terminology, we introduce shorthand notation, and define the wave equation: where the fields are assumed to have a harmonic frequency dependence, e −iωt , L = c 2 ∇ × ∇×, and t (r) is the total dielectric constant that we will assume is nondispersive. The operatorˆ t is defined as r|ˆ t |r = To construct a Green function solution, we consider a situation where we start with cavity 1, and then add cavity 2. The dielectric constant defining cavity 1 isˆ 1 = b +V 1 , so we can also write the wave equation as whereV 2 defines the dielectric constant change after adding in cavity 2, andˆ b defines the entire background without either cavity. Naturally, we can also start from cavity 2 and add in cavity 1, and the end Green function that includes both cavities must be the same.
B. Coupled mode theory and lossless exceptional points using normal modes Exploiting the fact that L is a linear self-adjoint operator over space, the homogeneous part of Eq. (7) defines an orthogonal set of eigenstates on a single cavity. It follows that where ω k are the eignfrequencies of the eigenstates f k . These states are also complete and orthogonal [90], so and the sum includes all modes, physical and unphysical. Next, we can formulate a scattering problem, based on some (arbitrary) reference input field, so that where G(r, r ) = r|Ĝ|r is the total Green function of the system (including both cavities), and is defined from For weakly coupled resonators, we can expandĜ in terms of a restricted set of carefully chosen basis states, which will be the dominant modes of the individual cavity systems [90]. Thus one obtains the Green function expansionĜ where both sums extend over all states of interest. If we obtain a solution for B α,β , then the scattering problem is solved as we have the total Green function, including the new poles of the coupled resonator system. In the absence of cavity 2, we define the cavity mode of the bare cavity 1 (i.e., no cavity mode 2 yet), from with the normalization f 1 |ˆ 1 |f 1 = 1. For simplicity, we assume the cavity supports a dominant single mode in the frequency of interest, but this can easily be generalized to allow for N modes per cavity. Similarly, we can define the solution of cavity 2, from with f 2 |ˆ 2 |f 2 = 1.
Subsequently, we substitute the mode expansions into the main Green function Eq. (11), to obtain where i, j refer to any basis states, andˆ α can refer tô 1 orˆ 2 . Equation (15) defines a matrix equation whose poles correspond to the new eigenfrequencies of the composite system. To proceed, we exploit the fact that the modes are only weakly coupled to each other, and so which can typically be easily checked numerically (else this approximation can also be relaxed, which just yields a more complex matrix to be solved). The matrix defined from Eq. (15), namely M BT = T , has the solution B α,β = [M −1 ] α,β , with elements: Thus, the matrix M is and we define the inter-mode coupling rate: so that Matrix inversion can be solved without approximations, however it is appropriate to obtain an easier form within a rotating-wave approximation. Using and we obtain an explicit solution for the Green function expansion coefficients where the pole frequencies are Finally, we note that for closed cavity systems, unitarity of a Hermitian system also requires that κ 12 = κ * 21 , and the pole frequencies are simply This concludes the derivation of the NM Green function with weakly coupled cavities. With regards to EPs, if we now consider the case with two cavity systems, one with loss: −γ (ω 1 = ω 0 − iγ), and one with a loss compensating gain: +γ (ω 2 = ω 0 + iγ), then one might be tempted to predict a situation where ω ± → ω 0 , if |κ 12 | = ±γ. The problem with this argument is that the original cavity modes here do not satisfy a Hermitian eigenvalue problem (assuming they are open and/or contain some loss), and thus the above coupled mode Green function solutions are not valid. Strictly, they are only valid for real eigenfrequency cavity modes. For very high Q cavities, however, the theory may be approximately correct, but the question and definition of a true EP still becomes questionable.
C. Coupled mode theory and lossy exceptional points using quasinormal modes Since we are interested in open cavities with loss and gain, we now adopt a more rigorous and appropriate resonator approach, using QNMs. One form of the QNM normalization can be defined from: where some coordinate transform has been applied to regularize the outgoing surface fields, e.g., through perfectly matched layers [65], but such terms will not be needed in the region for CMT overlap integrals as discussed below. Alternative QNM normalizations are discussed in Appendix A. These fields are now solutions to the eigenvalue problem with complex frequencies, and virtually all of the previous equations apply, with some simple replacements: (i) The eigenfrequencies become complex and formally discrete 1 : 1 Although the NMs are also assumed to be discrete for resonator problems, formally they yield continuous eigenfrequencies.
(ii) The QNM Green function expansion iŝ (iii) The completeness relation becomes which is assumed to be valid for spatial regions near or inside the scattering geometry.
(iv) For regions far outside the scattering geometry, one can use regularized QNMs (non divergent), such that where |F α is obtained from a Dyson solution using the original QNM [87], or using near-field to far-field transformations [91]. For high Q cavities however, using the QNMs in the perturbative cavity region is in practice extremely accurate, as we also confirm later.
With these replacements, we can use the same approach as before. Assuming again that the solution is first solved for cavity 1, and then we add in cavity 2, we derive: which notably now uses an unconjugated norm in the QNM overlap integrals. For the QNM formalism, notẽ κ 12 =κ * 21 , in contrast to NM theory. In QNM theory, these off-diagonal terms are not the complex conjugates of each other, since the open cavity system does not obey Hermiticity. Indeed, using such a relation with open cavities is simply ill-defined.
The full QNM Green function solution (using Eq. (30) and Eq. (27)) can be written as follows: with the two new QNM pole frequencies for the composite cavity system, In the limit of no coupling,κ 12 =κ 21 = 0, andω ± → ω 1,2 , and we recover the original Green function expansion for two separated resonators.
Finally, we briefly discuss the condition for EPs. Equation (33) shows that the exceptional point occurs when and in general one might obtain a lossy EP (at best), since the resonators are open, and the coupled system must also yield a finite loss for any linear amplifying medium [1]. The influence of dissipation on EPs and CMT is often argued heuristically. For example (defining 32, Miri and Alù discussed two limits where (i) the coupling is real, i.e., B = 0, so the EP condition would become A = ±2γ 0 , with γ 1 = γ 0 and γ 2 = −γ 0 (ω 1 = ω 2 = ω 0 ); and (ii) assume the coupling is purely imaginary, i.e., A = 0; then the EPs condition would become B = ±∆, where ω 1 = ω 0 + ∆ and ω 2 = ω 0 (γ 1 = γ 2 ). In practise it would be very difficult to reach exactly this regime (Eq.(34)), though one can likely come close, and certainly find signatures of EP behavior. Experimentally, in Refs. 37, 38, and 92, the spectra below and above EPs are observed, as well as the associated unidirectional transmission, but an exact EP was not demonstrated. Note, however, that these systems involve waveguides, so the CMT theory should be modified to connect more closely to such works. The theory in Ref. 45 also finds non-ideal EPs if one shows a zoom-in of the EP region in the complex frequency plane.
For a high Q resonator, one can assume that an approximate EP may be obtained when Re(κ 12 ) = γ 0 , but in general this condition is an approximate one for open resonators (and QNMs). Moreover, the lineshapes associated with the QNMs are far richer than with two coupled Lorentzian lineshapes, as we will demonstrate in more detail below, even for very high Q modes (Q ≈ 10 5 ).
D. Hybrid quasinormal modes using the couple mode theory Next, we introduce a model for obtaining the hybridized QNMs, which can be obtained analytically in terms of the uncoupled QNMs; this considerably simplifies the numerical solutions and helps to identify the underlying physics of how the hybridized modes are formed.
In terms of the hybrized eigenfrequenciesω ± , defined in Eq. (33), we obtain: or Assumingκ =κ 12 ≈κ 21 , then Subsequently, we also obtain the new Green function: which gives the same results as Eq. (32) 2 , but is now in diagonalized form.

IV. GREEN FUNCTIONS AND PURCELL FACTORS AT THE EXCEPTIONAL POINT
The SE rates (and generalized Purcell factors) can be significantly modified close to EPs [82][83][84][93][94][95], where a squared Lorentzian contribution has been emphasized as well as signatures of linewidth narrowing. With higher-order EPs, these linewidths may be reduced further [83], e.g., a cubic Lorentzian lineshape has been predicted with third-order EPs [94].
Below, we focus on the more general second-order EPs and first briefly comment on previous theoretical predictions about the modified lineshapes. In Ref. 84, the frequency dependent response of SE has the following form: which is a complex Lorentzian squared and a single Lorentzian, if very near the EP. A similar regime was predicted and shown numerically in Ref. 83. The main Green function response was predicted to have the following form: where C and D are connected with the Jordan vectors. They also give an approximate coupled mode theory expression for the Green function expansion.
Next, we will first show how our QNM Green function is fully consistent with the above predictions, and then show why one can find a much richer range of complex lineshapes, which we also demonstrate explicitly in the numerical results section.
We define the EP resonance frequency fromω ± = ω EP = (ω 1 +ω 2 )/2, which occurs when∆ 2 12 which is identical to what we have done already but this form allows us to expand the Green function near the EP resonances. For example, theĜ 1,1 term (expanded in terms of the bare mode from resonator 1) is: where we clearly see the separation of a Lorentzian and a Lorentzian squared contribution, in agreement with Refs. 82-84.
We stress that our spectral forms are not actually Lorentzian or Lorentzian-squared because of the QNM phase. Indeed, given the appropriate QNM phase, these terms can contribute negatively, a feature that is already known with coupled QNMs yielding a Fano resonance [91,[96][97][98]. It is also useful to compare with the response of the single cavity solution, which is simplŷ which has a single Lorentzian-like feature, again modified by the QNM phase. Thus the Lorentzian-squared feature is caused by the EP coupling regime. Sinceω EP = (ω 1 +ω 2 )/2, we can also write the QNM Green function solution to the coupled resonator EP regime aŝ Note that any divergences from the hybrid modes at the EP are avoided here, since we use an expansion in terms of the bare modes and CMT. Interestingly, we have obtained this known (and highly unusual) form without having to perform any Jordan expansion around the EP pole [83,84,93]. Assuming the CMT is accurate (and we show later that it can be quantitatively accurate), this is clearly a more convenient form to work with. It is also important to note that the hybrid QNM modes, e.g., in Eq. (37), diverge at the EP, since the QNMs at exactly this point are ill-defined and selforthogonal. In practise, however, this is not a restriction, as most solutions will deviate from precisely this point, where the two Green function responses become identical and well defined; thus one can use either the bare modes (non-diagonal form) or the hybrid mode solutions (diagonal form). This is a significant advantage that benefits from the constructed CMT.
In the next section, we will highlight much more general forms of the spectral lineshapes near the EP, which are fully verified by numerically exact solutions (within numerical precision). Our calculations also point out a fundamental problem with defining a SE rate in media with coupled loss and gain resonators, even when the total Green function is analytic in the upper complex half plane [1].

V. NUMERICAL RESULTS FOR COUPLED LOSS-GAIN MICRODISK RESONATORS
We consider two coupled loss-gain microdisk resonators, both with a disk radius of R = 5 µm (ct. Fig. 1). The refractive index of the lossy (gain) resonator is n loss = 2 + 10 −5 i (n gain = 2 − 5 × 10 −6 i), unless stated otherwise. The background medium is free space with n B = 1. The gap distance between the resonators is d gap (around 1120 ∼ 1200 nm). The dipole (out-ofplane line current, a point in 2D, shown as red dot in Fig. 1) is placed within the gap, at several possible positions: r d = r 1 , r 2 , r 3 , r 4 , and r 5 , where l = 10 nm and h = (d − 2l)/4 (r 3 is at gap center).

A. Quasinormal modes for single loss and single gain resonators
Before investigating the coupled resonator regime, we first show the QNMs for the single loss or gain resonators, which will be used as input for the CMT in the next subsection.
A single 2D microdisk resonator supports WGMs [99,100], which generally can be described by three integers: radial mode number q (=1,2,3...), azimuthal mode number m, and polarization p (TM or TE). Note that for the same q, m, and p, there are two degenerate counterpropagating modes [60,[101][102][103]: with counter clockwise (ccw) direction. for a TM mode, there is only the z component for electric fields), which share the same eigenvalues ( for a TE mode, the z component of the magnetic field has similar properties). A linear combinations of these modes result in degenerate standing waves [101][102][103], such as a symmetric standing mode E s = E cw + E ccw ∝ cos (mφ), and an antisymmetric standing wave mode E as = E cw − E ccw ∝ sin (mφ).
In general, the mode with q = 1 and m 1 has a high Q and strong field confinement. In this work, we focus on a TM mode (H x , H y , E z ) with q = 1, and m = 37 (yielding a resonant wavelength near the telecommunication band, around 1487 nm). To compute the normalized QNMs, we employ an efficient dipole scattering approach to obtain the QNMs in complex frequency [67], where an out-of-plane line current (a point in 2D, z-polarized) is placed at 10 nm away from the 2D microdisk (the dipole is at r d = r 1 for the single lossy cavity and at r d = r 5 for the single gain cavity). For more details, see Appendix A 1.
Conveniently, for our TM dipole location along the xaxis, we excite only one of the two degenerate standing wave modes, as shown in Fig. 2(a) for the single lossy cavity, which has the form cos(mφ) if one considers φ = 0 at the positive x-axis. The orthogonal and generate QNM has the form sin(mφ), constructed from an asymmetric linear combination of the clockwise and counter clockwise modes, though we only need to consider the symmetric QNM for our chosen dipole locations below. More details can be found in Appendix B. Also note that this working QNM dominates in the frequency region of interest below, because the two closest modes are q = 4, m = 25 and q = 3, m = 29, and the angular frequency spacing between them and the working mode are ∆ω = 3.43 THz and ∆ω = 19.6 THz, which are much larger than the FWHM (∼ 12.52 GHz for the single lossy cavity, and ∼ 6.26 GHz for the single gain cavity). Also note, the free spectral range for modes with q = 1 is also much larger than these FWHM values, e.g., the angular frequency spacing between q = 1, m = 36 (or q = 1, m = 38) with the working mode (q = 1, m = 37), would be around ∆ω = 31.6 THz. Thus, we can adopt a single QNM approximation for the mode of interest for each resonator, which we will also verify below by performing full dipole calculations with no mode approximations.
Numerically, the complex angular eigenfrequency for single lossy resonator is found to beω L = ω L − iγ L ≡ ω 0 − iγ 0 = 1.266666 × 10 15 − 6.260269 × 10 9 i rad/s, with a quality factor Q ∼ 10 5 ; here ω 0 (γ 0 ) is the real part (opposite imaginary part) of the complex eigenfrequency. The corresponding QNM field distribution (real part Re[f L z ]), is shown in Fig. 2 (a); the imaginary part is much smaller than the real part, and thus is not shown.
To better understand the overlap integrals for use in CMT, the QNMs at the second resonator region are also shown (1155 nm away, i.e., d gap = 1155 nm), which is hardly seen on a linear scale because they are very small compared with fields close to the resonator. This also indicates that the lossy QNMs show no sign of a spatial divergence in this region and thus they can be accurately used for input to CMT, which is a consequence of the high Q.
Similarly, the complex angular eigenfrequency for the single gain resonator is found atω Fig. 2(b). In addition, the Purcell factors from the dipole location (10 nm away from the loss and gain resonator), using a single QNM contribution, (Eq. (5)) agree quantitatively well with full dipole results (see Appendix D). We also highlight that the Purcell factors (as defined in a semi-classical model)  35)), for a separation distance dgap = 1155 nm; this regime is close to the EP as shown in (e) and (f). (e-f) Complex QNM eigenfrequencies of the coupled microdisk resonators from a direct numerical eigenfrequency solver in COMSOL (numerical, black circles) and the analytical CMT (Eq. (33), red stars forω+, green stars forω−), as a function of gap separation dgap. Note using the eigenfrequency solver, there are two bare degenerate modes per resonator, and four new eigenfrequencies for coupled resonators (four black circles for each gap distance), but with two pair of degenerate modes. For the dipole QNM technique, only one of degenerate standing-wave QNMs are obtained and used for single resonator (see text and Appendix A); thus, the analytical CMT produced two new eigenfrequencies in the presence of finite coupling (one red star and one green star for each gap distance). The arrows show two solutions near an EP. Note that ω0 and γ0 are for uncoupled single lossy resonator with n loss = 2 + 10 −5 i (ωL = ωL − iγL = ω0 − iγ0); while for the uncoupled single gain cavity, with n loss = 2 − 5 × 10 −6 i, the resonance is around at ωG ≈ ω0 + i0.5γ0.
a The cavity at the right is shown for clarity and for use in the CMT, but of course is not there for the single resonator calculations.
are net negative for single gain cavity only (since the field is being amplified rather than dissipated), a regime that will also be shown below for the coupled resonator system.

B. Hybrid quasinormal modes for coupled loss-gain resonators
We next study the hybrid QNMs formed from the two coupled microdisk resonators, using the QNM CMT and also using full numerical solutions to confirm the accu-racy of our semi-analytical results.
Using only the QNMs from single lossy/gain cavities as input, the properties for the coupled modes are obtained analytically. First, the new angular eigenfrequenciesω ± are computed from Eq. (33), where the coupling coefficientsκ 12 andκ 21 (both complex, Eq. (31)) are related to the overlap QNM integrals for different gap distances d gap . Note the input variablesω 1/2 andf 1/2 in theory parts are changed toω L/G andf L/G in the numerical parts accordingly. As shown in Figs. 2(e) and (f), the analytical eigenfrequencies (Eq. (33)) versus full numerical solution in COMSOL (eigenfrequency solver) show excellent agreement, and at all gap separations. This quantitative level of agreement is also obtained for the QNMs and the QNM Green functions as we will show in more detail below.
Note, with a direct QNM eigenfrequency solver, there would be four new eigenfrequencies for coupled resonators, because there are two degenerate standing modes per resonator (see Appendix B); thus, there are four black circles for each gap distance in Fig. 2 (5) or Eq. (A8))) for the hybrid QNM contributions as well as their sum. Also shown is the full dipole simulation with no approximations (red marker, Eq. (C1)), and the QNM contribution from a single lossy resonator (multiplied by a constant for clarity). The dipole is placed at r d = r1 (10 nm away from the lossy resonator). (b) Similar to (a), but with a different dipole position r d = r5, 10 nm away from the gain resonator.
Purcell factors are obtained analytically from Eq. (5), using the QNM Green function (Eq. (38)). As shown in Fig. 1, we consider five potential dipole positions along the x-axis, and study the Purcell factors as a function of frequency in each case.
For a dipole at the position r d = r 1 (10 nm away from the lossy cavity), the Purcell factors are shown in Fig. 3  all frequencies. We stress there are no fitting parameters in the QNM solutions. Moreover, from our theory, the contribution from the hybrid QNMsf + andf − can be shown separately, as indicated by the green and blue solid curves (Eq. (38), the first and second term). As a reference, we also show the Purcell factors for a single lossy cavity (orange dashed curve), which is net positive and multiplied by 8 for a better graphical comparison. The gain resonator clearly acts to suppress the broadening and enhance the overall Purcell factors.
To help explain these unusual lineshapes, we can study the contributions to the QNM Green functions from the phases for the two hybrid QNMs. Definingf + z (r) = |f + z (r)|e iφ + (r) andf − z (r) = |f − z (r)|e iφ − (r) , the QNM Green function can be expressed as from which we extract the imaginary part for use in Purcell's formula [98]: where we introduced the normalized Lorentzian lineshapes, . Similar expressions have been used to explain Fano resonances formed by elastic QNMs in coupled cavity beams [98]. For comparison, the Green functions for the singlemode bare QNMs are also given as with normalized Lorentzian lineshapes for the loss and gain resonator modes, and we have redefined the QNMs asf L z (r) = |f L z (r)|e iφ L (r) andf G z (r) = |f G z (r)|e iφ G (r) . In Fig. 3(a), we show the total Purcell factors for d gap = 1155 nm, as well as the hybrid mode contributions, and the lossy mode result on its own for comparison. First, the single mode case (orange dashed curve, single lossy resonator) shows a typical Lorentzian peak, which can be explained by the QNM phase contributions: cos [2φ L (r 1 )] = 0.999999 and sin [2φ L (r 1 )] = −0.001244 (also shown in Table I). Applying these to Eq. (48) results in a typical Lorentzian lineshape.
Next we focus on the hybrid modes (coupled resonator case), where now: cos [2φ + (r 1 )] = 0.1514, sin [2φ + (r 1 )] = Table I. Phase contributions in the imaginary part of QNM Green function, which is used to explain the lineshapes of the corresponding Purcell factors shown in Fig. 3 and Fig. 4.  Table I), which explain the non-Lorentzian lineshapes for the separate contributions (two terms in Eq. (46), green solid curve and blue solid curve shown in Fig. 3 (a)). Then combining the weights for each term (from |f + z (r 1 )| 2 and |f − z (r 1 )| 2 as shown in Eq. (46)), we obtain negative Purcell factors in a wide range of frequencies (black solid curve in Fig. 3 (a)).
Furthermore, these QNM phases will result in positiondependent lineshapes for the Purcell factors. Thus when we change the dipole position from r d = r 1 to r d = r 5 (10 nm away from the gain cavity), the Purcell factor lineshapes can change significantly, as shown in Fig. 3 (b) (again for d gap = 1155 nm), where the total contribution (black solid curve) are found to be mainly net positive in a wide range of frequency, and also show excellent agreement with the full dipole method (red circles). Once again, the contribution fromf + andf − can be given separately, as shown with the green and blue solid curves. In addition, the Purcell factors for the single gain cavity is shown as an orange dashed curve, which is net negative (and multiplied by 4 for a better graphical comparison).
Similar to case with a gap distance d gap = 1155 nm, the Purcell factors with d gap = 1160 nm for a dipole at r 1 are shown in Fig. 4(a). Their separate contributions from QNM + (green solid curve) and QNM − (blue solid curve) show Lorentzian lineshapes, but with different linewidth since γ + = γ − (see Fig. 2 (f)) and different signs due to cos [2φ + (r 1 )] = −1.0000 and cos [2φ − (r 1 )] = 1.0000. Combing these contributions, the total contribution are found to be net negative in a wide range of frequencies. When the dipole is at r 5 , the Purcell factors are shown in Fig. 4 (b), where the Lorentzian lineshapes for separate contributions could also be explained by the QNM phases cos [2φ + (r 5 )] = 0.9999 and cos [2φ − (r 5 )] = −0.9999. Considering their separate contributions, the total Purcell factors are now mostly positive. More details about the QNM phases are shown in Table I.
These QNM phases have a significant effect on the spectral lineshapes of the Purcell factors, which can be further seen in Fig. 5(a) (for d gap = 1155 nm) and (b) (for d gap = 1160 nm), where the total Purcell factors are given for five different dipole positions as shown in Fig. 1. For better comparison, the Purcell factors at r 2 , r 3 and r 4 are multiplied by 10, 300 and 10. Note, we do not show the full dipole results here, as they are basically indistinguishable from the the two QNM solutions, similar to the previous comparisons.
We stress that these negative Purcell factors are unphysical because the classical Purcell factor formulas, and Fermi's golden rule, is no longer working with a gain medium [85], though the local density of states (LDOS, proportional to Im[G(r d , r d , ω)]) is correct. Reference 85 discusses the breakdown of classical formulas in more detail, and shows that the correct SE rates must be described fully quantum mechanically, in which case the quantum Purcell factors (modified SE decay rayes) are indeed net positive.

D. Green function Propagators
This next section presents example Green function propagator results, namely |G(R, r d , ω)|, which can be useful to relate to a range of experimental observables. For example, this function is required to model the spectrum at a point detector located at R, emitted from a dipole at r d . Unlike the (classical) Purcell factor, it is a well defined quantity for use in both classical and quantum field theory, and is used frequently in the exploration of light-matter coupling regimes [104][105][106]. For example, considering an excitation dipole, d, then the spectrum is S(R, ω) ∝ |G(R, r d , ω) · d| 2 .
A far-field detection point is first chosen at R =   Fig. 6(a), for d gap = 1160 nm. For comparison, the propagator for a single lossy (gain) cavity is also shown as green (red) dashed curve, where the dipole is placed at r d = r 1 (r 5 ). Furthermore, for the far field detected point at R = R 2 = [14.09, 0] µm, the propagators are shown in Fig. 6(b).
To better explain the physics of the propagator lineshape, we once more write out the two QNM expanded Green function in terms of the hybrid modes, where again we can recognize non-Lorentzian features from the QNM phase terms, now from both the dipole position and the detector position. The detailed phases are shown in Table II. Note, the two-space point Green function has a range of uses for describing light-matter interactions, including the description of photon transport [104,107,108], and can be used to model collective effects with multiple emitters and dipoles.  (50)).

VI. QUASINORMAL MODES FOR INDEX-MODULATED RING RESONATORS NEAR EXCEPTIONAL POINTS
In order to test our general theory even further, and motivated by the recent experiments in Ref. [81], we next study QNMs for EP-like behaviour from indexmodulated ring resonators. In these experiments, the authors reported a novel reversible chiral emission for linearly polarized dipoles coupled near the ring resonators, using both optical and acoustic systems. One of the novel conclusions was that the emitter becomes disentangled from the system eigenmodes, and coupled to a "missing dimension", such as the Jordan vector. Here, we study a similar structure, and show that the same features naturally emerge and are fully expected using a two QNM picture, which is again due to the special properties of the QNM phases and how they couple and hybridize.
In the following, we will keep n 0 = 3.0 and δn Re = 0.001, while δ Im is changed in the range of 0.0005 ∼ 0.0015.
The field distributions Re[f s ], withf =f xx +f yŷ = f sŝ +f φφ , of QNM 1 and QNM 2 for the ring structures, with δn Re = δn Im = 0.001, are shown in Fig. 7(b) and (d). The black square in Fig. 7(b) is enlarged and shown in Fig. 7(c), where the detailed periodic refractive index is displayed. In addition to the general (x,ŷ) basis, we also utilize a polar (ŝ,φ) basis. Their unit vectors are related fromŝ =x cos φ +ŷ sin φ andφ = −x sin(φ) +ŷ cos φ (the positive direction of polar angle φ is along CCW direction). As mentioned before, with the dipole QNM technique, we use a linearφ dipole at φ = 0 (the bottom one d 1 in Fig. 7(c)). Once we obtain the QNMs, then we can easily couple them to any in-plane linear dipoles with any polarization; we schematically give two example dipoles in Fig. 7 (c), one is a linearŝ dipole at φ = φ 2 = π/4m (the center dipole, d 2 ), and the other one is a linear φ dipole at φ = φ 3 = 3π/4m (the top dipole, d 3 ).

A. Purcell factors versus frequency for different dipole locations
To confirm the accuracy of our two QNM description for the index-modulated ring resonators, we calculate the Purcell factors analytically from the QNM properties (Eq. (5)) for aφ dipole at φ = 0 (equal toŷ at φ = 0, d 1 in Fig. 7 (c)), and compare them with the full dipole method (Eq. (C1)). For all configurations investigated, an excellent agreement is shown as a function of frequency as demonstrated in Fig. 8, where the results for six cases are investigated with the analytical two QNM expansion are studied, including δn Im = 0.0005, 0.0008, 0.001, 0.0011, 0.0012, 0.0015 while keeping δn Re = 0.001. Their eigenfrequencies are marked in Fig. 7(e). Although there are no material gain regions in these resonators, we also note the appearance of modal negative Purcell factors. This level of agreement is unusual in the sense that the eigenmodes in Ref. [81] were reported to be completely decoupled, which was used to argue a chiral emission, an effect that is not expected nor . In addition to general unit vectorsx andŷ, we also utilizeŝ =x cos φ +ŷ sin φ andφ = −x sin φ +ŷ cos φ. Three example dipoles are shown, where one (d1) is a linearφ dipole at φ = φ1 = 0 (i.e.,ŷ dipole at φ = 0), one (d2) is a linearŝ dipole at φ = φ2 = π/4m, and the left one (d3) is a linearφ dipole at φ = φ2 = 3π/4m. (e) Eigenfrequencies from QNMs (ω1 = ω1 − iγ1 andω2 = ω2 − iγ2, red and green stars) and approximate COMSOL solver (black circles) for the ring resonator as a function of δnIm when keeping δnRe = 0.001. The insert in the bottom figure shows the zoom in region close to δnIm = 0.0005. One will find that there are two QNMs when, e.g., δnRe = δnIm = 0.001. Specifically, we will focus on six cases (including δnIm = 0.0005, 0.0015 and four more cases that the arrows points to) and compare the results from QNMs and the full dipole simulation, where an excellent agreement is shown for the Purcell factors (cf. Fig. 8). Note that the eigenfrequency of QNM 1 for δnIm = 0.001 is ω0 − iγ0 ≈ 2.50543552 × 10 15 − i3.768337 × 10 11 (rad/s), i.e., ω0 and γ0 are the real part and the opposite imaginary part of this pole eigenfrequency. With the increase of δnIm, the decay rates of these QNMs increase.   -p))φ dipoles at φ = π/4m, φ = 3π/4m, φ = 5π/4m, and φ = 7π/4m (the radial distance between the dipole and the ring surface is 10 nm (δnRe = δnIm = 0.001). The red dot with double black arrows schematically describe the position and polarization of the linear dipoles. The Poynting vectors are labelled as the magenta (for CW direction) and blue (for CCW direction) arrows. The distributions are normalized to S Poyn φ (for those from full dipole method, they are with arb. units). These power flows are calculated at a example real frequency ω = 2.50543 × 10 15 (rad/s) (close to QNM resonance). In the ring region, one will find S Poyn φ < 0 for CW propagation and S Poyn φ > 0 for CCW propagation; the net power flow from QNMs approaches are in excellent agreement with the full dipole simulations. obtained for a regular ring resonator. We will connect to the power flow emission below, and show that it is also well explained in terms of the underlying QNMs.

B. Chiral power flow from linearly polarized dipoles
To directly connect with the experimentally measured chirality in Ref. 81 from linear dipole emitters, here we show chiral power flow using only the QNM propagators in a similar microring with refractive index modulation. We can easily obtain the power flow by computing the scattered fields and the Poynting vector from QNMs and the QNM propagator, and thus the calculations are basically instantaneous.
Once the QNMs are known, the scattered field (in real frequency space) at any point r from a dipole at r d is simply where d = d 0 n d , and we will consider bothŝ dipoles andφ dipoles (i.e., n d =ŝ or n d =φ). The magnetic scattered field is then and thus the Poynting vector is simply where S Poyn φ (r, ω) is the projection alongφ, and one will find S Poyn φ (r, ω) < 0 (> 0) when it goes along the CW (CCW) direction. Using the QNMs, we stress that we can compute the direction of the power flow and scattered fields analytically, for any dipole position and orientation. This further demonstrates the power of having an analytical expression for the Green function in terms of only QNM properties.
In this section, we focus on the case with δn Re = δn Im = 0.001, which is closest to the EP as seen from Fig. 7(e) and Fig. 8(c). Figure 9(a-h) shows the chiral power flow (Eq. (54)) from our QNMs model, forŝ andφ dipole at several specific positions, including φ = π/4m (center of two lossy sections, see Fig. 7(c)), φ = 3π/4m (center of two lossless sections), φ = 5π/4m (center of two lossy sections), and φ = 7π/4m (center of two lossless sections) (the radial distance between the dipole and the ring surface is 10 nm). The red dot with a double black arrows in Fig. 9 schematically shows the position and polarization of these linear dipoles (though this is not to scale). Also note that these power flows presented in Fig. 9 are calculated at a example real frequency ω = 2.50543 × 10 15 (rad/s) (close to QNM resonance), and we stress that the chirality will remain the same in a range around ω 0 ± 1.5γ 0 .
The first thing to observe, is that when a linearŝ dipole is placed at φ = π/4m ( Fig. 9 (a)), the net power flow (S Poyn , Eq. (54)), indicate by the magenta arrows, goes along the CW direction. The distribution shown is normalized to the projection S Poyn φ . We also find that S Poyn φ < 0 in the ring region, which further confirms the net energy flow goes along the CW direction. In addition, the azimuthal period of such phenomenon is π/m, i.e., when a linearŝ dipole is located at φ = π 4m + p π m (center of the two lossy regions), where p = 0, 1, 2, ..., the net power will go along the CW direction, such as at φ = 5π/4m shown in Fig. 9 (c).
To further justify that our QNM picture is both correct and rigorously accurate for describing this chiral emission, we have confirmed these unusual properties with full dipole calculations shown in Fig. 9 ((i) and (k)). This also supports the experimental findings in Ref. 81, where they measured CW propagation when putting a linear dipole at φ = π/4l (center of lossy region) (azimuthal mode number l = 1 in their considered structure, and we are using azimuthal mode number m = 21 for a similar ring structure). However, our interpretation is drastically different. Instead or arguing a decoupling from the eigenmodes and coupling to a missing dimension, we find that our net chiral power flows are fully explained from the underlying QNMs, which are the correct natural eigenmodes of such resonators 3 .
Moreover, we also find that such unusual properties are not limited toŝ dipoles at φ = π 4m + p π m . As shown in Fig. 9(b,d) (QNMs results) and Fig. 9 (j,l) (full dipole results), when a linearŝ dipole is placed at φ = 3π 4m + p π m (center of the two lossless sections, where p = 0, 1, 2, ...), the net power flow is along the CCW direction. In addition, as shown in Fig. 9 (e-h) (QNMs results) and Fig. 9 (m-p) (full dipole results), when a linearφ dipole is placed at φ = π 4m + p π m (φ = 3π 4m + p π m ), the net power flow will go along the CCW (CW) direction. Also note, forŝ dipoles at φ = n π m and φ = 2π 4m + n π m (cf. Fig. 7), respectively, the power flow direction will change from CCW to CW and CW to CCW, around resonance, i.e., at resonance, there is no net flow direction. This is also true forφ dipoles; forφ at φ = n π m (φ = 2π 4m + n π m ), the power flow direction will change from CW (CCW) to CCW (CW). This gives one an external control to change the directionality by simply changing the frequency.
Note that such chirality behaviour is not found in regular ring resonators (without a refractive modulation) for linearŝ andφ dipoles, i.e., no net power flow is obtained. However, by using a circular dipole close to or inside the general ring resonator, chirality will exist, as shown in Ref. 109, where positional dependent chirality for rightor left-handed circular dipoles are demonstrated. This is also similar to how one excites unidirectional propagation in photonic crystal waveguides [110,111]. However, in all these cases, local symmetry breaking is possible by using a circularly polarized dipole.

VII. CONCLUSIONS
We have introduced a powerful and highly accurate QNM approach to coupled loss-gain resonators, and have presented a rigorous and intuitive CMT based on the photonic Green function, which allows one to solve the coupled system efficiently with just the bare QNM solutions from the individual resonators. We have also highlighted the failure of using a NM CMT approach when defining the general conditions for finding EPs.
For the SE response of embedded dipoles in these systems, we have carried out detailed calculations for coupled microdisk resonators; as well as finding Lorentzianlike and Lorentzian-suqared like responses at the EP, consistent with other works, we have shown much richer Purcell factor lineshapes near EPs for various designs and spatial dipole positions, showing excellent agreement with QNM CMT and full dipole calculations. In particular, we have shown how the Purcell factors can also be negative, in loss-gain media, even when the hybrid modes are both lossy (γ µ > 0). This is caused by a breakdown of Fermi's golden rule which incorrectly assumes that the SE rate is propositional to the (projected) LDOS) [85]. In addition, we also showed how the Green function propagators (related to various experimental observables, such as the emitted spectrum) also take on rich non-Lorentzian features, which depend on the values of the QNM phases.
In addition to the coupled loss and gain resonators, we also investigated EP-like resonances formed from indexmodulated ring resonators, where unusual chiral emission from linearly polarized diploes was recently observed [81]. Once again, we showed that the full dipole Purcell factor response is quantitatively well explained in terms of the main two QNMs of this resonator. Moreover, when a linearŝ dipole is located at φ = π/4m, the Poynting vectors from QNMs propagators goes along the CW direction, which supports the experimental findings in Ref. 81 for similar ring resonators. Notably, our explanation is in contrast with the view that the emitter does not couple to the system eigenmodes. In addition, we also showed that such chirality is not limited toŝ dipoles at φ = π 4m + p π m , and we also show the opposite chirality forŝ dipoles at φ = 3π 4m + p π m . There is also similar chirality for linearφ dipoles at these positions. We stress, again, that these net power flows can be well explained and interpreted from the two underling QNMs and the corresponding Green function, where the QNM phases play a decisive and fundamental role on the light emission, without having to invoke any unusual interpretation such as a missing dimension (Jordan vector).
Apart from providing a detailed and intuitive formalism for understanding the classical mode properties of these complex coupled resonator systems, our QNM formalism forms the basis for a rigorous quantum optics approach in media with gain and loss, using new approaches with quantized QNMs [74][75][76], which has already lead to a revision of the usual photonic Fermi's golden rule for coupled loss and gain resonators [85], one in which the net Purcell factor is always a positive quantity.

VIII. ACKNOWLEDGEMENTS
We acknowledge funding from Queen's University, the Canadian Foundation for Innovation, the Natural Sciences and Engineering Research Council of Canada, and CMC Microsystems for the provision of COMSOL Multiphysics. We also acknowledge support from the Alexander von Humboldt Foundation through a Humboldt Research Award. We thank Andreas Knorr for discussions and support, and Thomas Christopoulos for discussions about PML properties in COMSOL.

Appendix A: Quasinormal mode (QNM) normalization
There are several general numerical approaches to obtaining normalized QNMs [71], including a dipole excitation technique in complex frequency [67], PML normalization [65], finite domain normalization with a surface term [72,73,86], finite-difference time-domain methods [89], and Riesz-projection-based techniques [68]. In the main text, we use the dipole technique in complex frequency space [67], and more details are shown in the following subsection A 1. A brief description of two alternative methods is also given in subsections A 2-A 3; numerically, we find all three approaches give the same normalized QNMs for the QNMs used in this work (within numerical precision).

QNM normalization and numerically exact
Green function from a dipole source As shown in main text, with CMT and the Green function theory, only the bare mode solutions (for a single lossy resonator or gain resonator) are required. We employ an efficient dipole scattering approach to obtain the uncoupled QNMs in complex frequency [67], where an out-of-plane line current (a point in 2D, red dot in Fig. 1, for a TM mode; it's a in-plane point dipole for a TE mode) is placed close to the lossy resonator, or the gain resonator, or the index-modulated ring resonators.
We can also use this approach for the coupled resonator problem, which we also do to check the accuracy of the analytical CMT for the hybrid modes (see Appendix E).
The scattered field of this dipole at r 0 is related to the 2D Green function, from where the units of the scattered field E scatt , 2D Green function G 2D , and 2D dipole moment d 2D are, respectively, V/m, m −2 , and C.
Expanding the 2D Green function with one QNM (dominating in the regime of interest), then so that Multiplying Eq. (A3) with d 2D and using r = r 0 , then Substituting this back to Eq. (A3), we obtain the 2D normalized QNM field as a function of spacẽ The above QNM simulations are performed in the commercial COMSOL software [112], where the frequency in Eq. (A5) is set as ω = (1 − 10 −7 ) ×ω c (or (1 − 10 −6 ) ×ω c , (1 − 10 −8 ) ×ω c , adjust with the quality factors), very close to the pole frequency. The computational domain (including PMLs) is around 804 ∼ 814 µm 2 (various gap distance for microdisk resonators; for microring resonators, it's around 72 µm 2 ), where the maximum mesh element sizes are 0.1 nm, 40 nm and 75 nm at the dipole point, inside and outside the 2D microdisks. To minimize boundary reflections, we used 15 layers to form the PMLs with a total thickness of 1.5 µm, which was found to be well converged numerically.

Perfectly matched later (PML) normalization
The PML normalization [65,70,113] approach is another alternative way to get the normalized QNMs, which is given via (for dispersive and nonmagnetic materials) whereh µ (r) is the magnetic field of QNM, V is the whole simulation region, and V PML denotes the PML region.
Outside the PML region, then the fields are zero. For a nondispersive and nonmagnetic material, and extra care is needed for the PML region (second term in Eq. (A10)), though this contribution can be very small for certain problems and geometries. In general, there are several kinds of transformation performed in PMLs to minimize boundary reflections. The first one is using special permittivity PML (r, ω) and permeability µ PML (r, ω) values, which is not always available with some commercial software [70]. A second approach uses a coordinate transformation, which is what we used here, with the built-in stretched-coordinate PML of COMSOL, where the coordinates are transferred from real space to the complex plane [114].
To verify that the dipole normalization technique and PML normalization are consistent with each other, we performed the norm calculation (Eq. (A10)) using the QNM fields from the dipole technique. We obtained f L (r)|f L (r) = 1.0004−0.0011i ∼ 1 for the QNM of interest from the single lossy WGM resonator with n loss = 2 + 10 −5 i, and f G (r)|f G (r) = 0.9991 + 0.0006i ∼ 1 for the QNM of interest from single gain WGM resonator with n gain = 2 − 5 × 10 −6 i. However, note, that for our chosen PML geometry and resonator, the contribution from the PML is practically negligible for this problem of interest (high Q resonance).

Finite-domain normalization with a surface term
A third way to normalize QNMs is through a finite domain normalization with an outgoing surface term [71][72][73]86], is given by and the surface ∂V is the surface of the finite domain V that surrounding the entire resonators, and ∂s is the derivative in the normal direction of the surface ∂V . Now the integral region V has a greater freedom of choice, but at the cost of adding a surface integral (though for 2D, it is a surface integral and line integral). For a nondispersive material, For completeness, we have checked our normalized QNMs with this technique, and obtain f L (r)|f L (r) = 1.0004 − 0.0011i ∼ 1 for the single bare QNM from single lossy WGM resonartor with n loss = 2 + 10 −5 i, and f G (r)|f G (r) = 0.9991 + 0.0006i ∼ 1 for the single bare QNM from single gain WGM resonartor with n gain = 2 − 5 × 10 −6 i. Again, we note for this problem, the contribution from the line integral is negligible, thus we obtain basically the same answer as the PML normalization. As a further check, we also tried above norm with two smaller domain with area of 579 µm 2 and 497 µm 2 , and obtain exactly the same numerical results as above.

Appendix B: Degenerate QNMs for microdisk Resonators
The WGM resonator naturally supports degenerate modes because of circular symmetry. For a specific azimuthal mode number m, the degenerate counterpropagating modes are [101][102][103]: E cw (r, φ) = E(r)e (−imφ) with clockwise (cw) direction, and E ccw (r, φ) = E(r)e (imφ) with counter clockwise (ccw) direction; for a TM mode, we only have the z component for the electric fields. A linear combinations of these fields result in standing waves [101][102][103], such as a symmetric standing mode E s = E cw + E ccw ∝ cos (mφ), and an antisymmetric standing mode E as = E cw − E ccw ∝ sin (mφ). These standing waves could of course also be normalized, and can be expanded to obtain the QNM Green function.
Using the direct eigenfrequency solver (where the boundary condition is PML) in COMSOL, one can obtain a pair of unnormalized standing modes distributioñ f M1/M2 un . This is not used generally as a QNM solver since it is not so robust as a nonlinear eigenvalue solver, but the results appear to be reasonable for the high Q modes of interest in this work.
For example, one can use PML normalization to obtain the relevant QNM from where f Mi un (r)|f Mi un (r) PML takes the form of Eq. A10. For a single lossy WGM resonator with refractive index Fig. 10(b) and (c). To better check their differences, we show them along the edge of the disk (equator) along the small vertical arrow (see Fig. 10(b) and (c)) direction in Fig. 10(d) (green solid and orange solid curves). One can find their antinode positions are off by π 2m (m = 37), corresponding to the difference between cos mφ and sin mφ. However, their absolute phases are offset, since these are related to where φ = 0 is defined in coordinate system (x M , y M ), with the eigenfrequency solver.
In contrast, in the dipole excitation technique [67], only one of the two degenerate standing wave modes is excited, as shown in Fig. 10(a) (which is the same as Fig. 2 (a) in main text), where the dipole is placed at Figure 10. Spatial profile of the degenerate QNMs for a single lossy resonator. (a) QNM obtained with a general QNM approach (Eq. (A5)) (it is the same as in Fig. 2 (a)), where only one of the degenerate standing wave modes is excited, which is with a cos(mφ) dependence, considering φ = 0 at the positive x-axis. r d = r 1 . The field Re[f L z ] along the equator is also shown in Fig. 10(d) and (e) (red solid curve), which is with the form cos(mφ) considering φ = 0 corresponds to the positive x-axis. Moreover, with the dipole technique, one can control the phases accurately by setting the dipole at different locations, as they naturally excite the same standing wave defined from that position. For example, if the dipole is placed at a position corresponding to φ = π/2m (r d = r 1 ), or φ = 2π/2m (r d = r 1 ), or φ = 3π/2m (r d = r 1 ) (the distance between the dipoles and the tangent of the disk surface remains at 10 nm) (see schematic in Fig. 10 (f)), then the same QNM is obtained with a simple shift in phase; these are shown in Fig. 10 (e) (green dashed, black dashed and blue dashed curves), with the phase terms cos(mφ − mφ ) = cos(mφ − π/2), cos(mφ − mφ ) = cos(mφ − 2π/2), and cos(mφ − mφ ) = cos(mφ − 3π/2) (considering φ = 0 is defined at the positive x-axis).

Appendix C: Full dipole simulations
To check the validity of the Purcell factors from the QNMs (for microdisks and microrings) and also the CMT Green function solutions (for microdisks), we compare these directly with the numerical results from a full dipole method (namely, with no modal approximations), which is obtained as follows: F num,2D P (r 0 , ω) = Lcn · S dipole,total (r, ω)dL c Lcn · S dipole,background (r, ω)dL c , (C1) where L c is a small circle (with radius 1 nm) surrounding point current (for TM mode; for 2D TE modes, it's a inplane point dipole) andn is a unit vector normal to L c , pointing outward. The vector S(r, ω) is the Poynting vector at this small circle and the subscript 'total' and 'background' represent the case with and without the resonator.
Appendix D: Checking the accuracy of the QNMs for a single lossy resonator or a single gain resonator In order to compare with the coupled resonator results, and to confirm the accuracy of a single QNM approach for the single resonators, we show the results for single lossy cavity with n loss = 2 + 10 −5 i or single gain cavity with n gain = 2 − 5 × 10 −6 i. Using the dipole technique described above (Eq. (A5)), the QNMsf L z andf G z are shown in Fig. 2 (a) and (b) separately. The dipole is put at 10 nm away from single lossy cavity or the gain cavity.
The corresponding Purcell factors are shown in Fig. 11, which show quantitative agreement with full dipole method. The Purcell factors from single lossy cavity are net-positive and those with single gain cavity are netnegative. Note these results are also shown as dashed curves in Fig. 3 and Fig. 4 in the main text, but multi- plied by constant for better comparison with the results from coupled system.

Appendix E: Direct QNM approach for coupled loss resonators and gain resonators
Naturally, one can also solve the coupled system with a QNM approach directly, instead of using coupled modes theory after the bare solutions are known. In Fig. 12, we show their comparison for the specific example shown in Fig. 3(a)   With the direct QNM approach (when the dipole is located at r d = r1), there are two dominant modes, QNM + and QNM − . Their separate contributions to the total Purcell factors are shown as a dashed magenta curve and a dashed cyan curve, and the total one is shown as a green dotted curve. They show good agreements with CMT results and full dipole results, i.e., the contribution from QNM + (green solid curve, see Fig. 3(a)) and QNM + (dashed magenta curve) agrees very well, as does the contribution from QNM − (blue solid curve) and QNM − (dashed cyan curve); the CMT and two direct QNM totals are practically indistinguishable (black solid curve, green dotted curve; and red circles show the full dipole results). contributions), as well as quantitatively good agreement between the direct QNM results and the full dipole simulations (green dotted curve and red circles).

Appendix F: Additional loss-gain resonator examples with different gain coefficients
Finally, we also show the Purcell factors for two additional loss-gain resonator examples, for different amounts of gain. The refractive index for the lossy resonator is fixed at n loss = 2 + 10 −5 i, as in the main text. When n gain = 2 − 6 × 10 −6 i, the Purcell factors with a gap distance d gap = 1145 nm (close to the lossy EP) for a dipole at r 1 (10 nm away from the lossy cavity) are shown in Fig. 13(a), which show very good agreement with the full dipole method. Again, we see that negative Purcell factors are obtained over a wide frequency range. The separate contributions fromf + andf − are also given. For better comparison, the Purcell factors with single lossy cavity are shown as an orange dashed curve (the dipole is at r 1 ), which is net positive and multiplied by 30 for clarity.
Similarly, the corresponding results for the case with n gain = 2 − 7 × 10 −6 i are shown in Fig. 13(b), where the gap distance is d gap = 1133 nm. Excellent agreement with full dipole results are also obtained. The absolute values of the negative Purcell factors increase here mainly because it is closer to the EP and the gap distance is smaller.