Quantized quasinormal-mode description of nonlinear cavity-QED effects from coupled resonators with a Fano-like resonance

We employ a recently developed quantization scheme for quasinormal modes (QNMs) to study a nonpertur-bative open cavity–QED system consisting of a hybrid metal-dielectric resonator coupled to a quantum emitter. This hybrid cavity system allows one to explore the complex coupling between a low-Q (quality factor) resonance and a high-Q resonance, manifesting in a striking Fano resonance, an effect that is not captured by traditional quantization schemes using normal modes or a Jaynes-Cummings (JC)–type model. The QNM quantization approach rigorously includes dissipative coupling between the QNMs and is supplemented with generalized input-output relations for the output electric ﬁeld operator for multiple modes in the system and correlation functions outside the system. The role of the dissipation-induced mode coupling is explored in the strong coupling regime between the photons and emitter beyond the ﬁrst rung of the JC dressed-state ladder. Important differences in the quantum master equation and input-output relations between the QNM quantum model and phenomenological dissipative JC models are found. For the hybridized high-Q cavity mode, we show how the dissipation-induced coupling causes a signiﬁcant reduction in the cavity-emitter coupling rate, and the cavity decay rate, compared to a simpler JC model. In a second step, numerical results for the Fock distributions and system as well as output correlation functions obtained from the quantized QNM model for the hybrid structure are compared with results from a phenomenological approach. We demonstrate explicitly how the quantized QNM model manifests in multiphoton quantum correlations beyond what is predicted by the usual JC models. DOI: 10.1103/PhysRevResearch.2.033456


I. INTRODUCTION
Quantum emitters coupled to photons/plasmons in dissipative nanostructures, such as micropillars [1][2][3] , photonic crystal cavities 4,5 , dielectric microdiscs 6 or metallic nanoparticles [7][8][9][10] , constitute an important field in quantum optics and quantum plasmonics.New classes of nanophotonics structures exhibit enhanced light-matter coupling suitable for studying a range of interesting phenomena and applications, such as non-classical light generation 11,12 , spasing [13][14][15] and quantum information processing 16,17 .While the dielectric cavities have high quality factors with state-of-the-art values 18 around Q ∼ 10 5 , metallic cavities are significantly more lossy (typically Q ∼ 10) due to Ohmic heating, but still provide comparable light-matter coupling regimes thanks to the strong local field confinement below the diffraction limit 19 .Other recent important lossy structures are hybrid metal-dielectric resonators, made from metals and dielectric resonators [20][21][22][23] , which combine the attributes of both resonator types and exhibit Fano interference effects between both systems.
For the theoretical description of quantum light-matter interactions, bound photon states in such systems are often treated using "normal modes" with real eigenfrequencies, typically using Jaynes-Cummings (JC) type models 24,25 .Additionally, since these structures are lossy, dissipation is usually introduced into the model by phenomenologically adding decay rates for the subsystems, typically for the optical modes and quantum emitters 26 .In contrast to this phenomenological approach, so-called quasinormal modes [27][28][29][30][31] (QNMs) with complex eigenfre-quencies (including loss in the imaginary part), intrinsically describe the open/lossy system by solving Maxwell equations with open boundary conditions.This approach allows one to determine useful cavity properties such as the radiative beta factors (quantum efficiencies), quality factors, and effective mode volumes [32][33][34][35][36][37] .With continued developments in computational electromagnetics, the numerical solution of the Helmholtz equation to obtain QNMs is becoming better unraveled and common today 33,36 , but its subsequent quantization using QNMs to retrieve well-known model systems, e.g., a microscopically defined JC model is still highly nontrivial.
The lossy and non-Hermitian character of these open systems prevents the use of a canonical quantization procedure for the discrete modes of interest (cf.Refs.28,  38, and 39).Recently, a general quantization scheme for three-dimensional absorptive and lossy media using QNMs as the basis for the field expansion was presented 40 , based on a Green function quantization approach 41,42 .It was demonstrated that, for more than one QNM, a coupling between the QNMs is induced by the dissipation through proper quantization in the dissipative system.The off-diagonal coupling is especially interesting for mode interference effects in the above mentioned hybrid structures, which is in the semi-classical model a consequence of the complex-valued QNMs 21 , and leads to highly non-Lorentzian line shapes.In fact, Franke et al. 40 demonstrated, in the single-photon limit (weak excitation), that such an interference effect can only be reproduced through the off-diagonal QNM coupling, when starting from a quantized mode approach.This was further confirmed in Ref. 43 for a Fano cavity, by a inde-pendent method and calculation based on introducing a phenomenological mode coupling in a two-mode master equation.
Recently, the approach from Ref. 40 was applied to accurately describe single-photon emission in a single-mode metal resonator 44 and was also used to model the photonic mode quantization for molecular optomechanics in a hybrid metal-dielectric resonator 22 .In the first case, a cavity output field expression for the single-QNM was derived, which is the basis to determine correlation functions and light statistics of the resonator-emitter system, important to simulate experimental situation, such as the Hong-Ou-Mandel 45 or Hanbury-Brown-Twiss 46 setups.Importantly, the QNM quantization scheme allows one to distinguish between radiative and non-radiative decay processes, which both enter naturally into the formalism through the same calculated QNMs.This separation is essential to describe a realistic input-output formalism in cases of absorptive lossy structures, where the output is usually treated in the same way as for systems with only radiative losses [47][48][49] .While in the single-mode case, the results from a more phenomenological dissipative JC model is modified by a loss-induced prefactor, which separates between radiative and non-radiative decay 44 , there are additional changes in the multi-mode case due to offdiagonal mode interaction, which may also effect the output coupling.This may also effect the behavior of multimode systems in higher rungs of the JC ladder, e.g., the change of Poissionian to sub-Possionian light when coupling a resonator-emitter system to another resonator, which is often described with two uncoupled modes 50,51 .
In this paper, we study the nonlinear multiphoton cavity-QED effects of a hybrid cavity structure containing a single quantum emitter, depicted in Fig. 1(a,b), which consists of a metal ellipsoid dimer on top of a photonic crystal beam.In contrast to the hybrid structure used in the Ref. 40, the resonator-TLS system in Fig. 1(a,b) is in the intermediate coupling regime, where the bad cavity approximation is no longer valid.Therefore, processes on higher rungs, i.e., many photon effects, in the anharmonic JC-model are more accessible and we analyse the effects of the inter-mode coupling in this multiphoton regime.To calculate measurable quantities, such as the second-order correlation function, an expression for the electric field operator outside the resonator for multiple QNMs is derived.A second main objective of this work is to compare the results for the Hamiltonian and Liouvillian of the quantized QNM model with a phenomenological dissipative two-mode JC model in the few photon limit.The phenomenological model assumes two uncoupled bosonic modes.It will be shown, that the off-diagonal coupling present in the quantization of the QNMs will induce drastic changes in the master equation and density matrix simulations for higher rungs of the JC ladder, which will further underline the importance of a quantized QNM model.
The rest of our paper is organized as follows: In Section II, we present the theoretical framework for The photonic crystallike mode (violet) is coupled to the plasmonic-like cavity mode (yellow) and the TLS (red) is coupled to both symmetrized QNMs.d) Excitation and input-output scheme.The TLS is driven by a cw-pump field with strength ΩL and decays with a rate γSE, while the QNM subsystem is characterized by the decay rates Γ pl and Γpc, naturally entering the model through the input fields A in pl and A in pc , respectively (see text).
this paper.First, we summarize the phenomenological Green function quantization approach for general spatialinhomogeneous and absorptive media; this includes a coupling between the medium-assisted electromagnetic field and an emitter.Second, the definition and main aspects of the QNM approach will be revisited and clarified.Third, we will briefly recapitulate the QNM quantization from Ref. 40.Then, the equations of motion for the QNM operators and input-output relation for multiple QNMs in the system are derived, and, based on that derivation, the QNM master equation with an additional external pump term is presented.Last, we derive the multi-mode output electric field operator to express the output correlation function in terms of QNM operators.
In Section III, we present three-dimensional numerical results for a metal-dielectric hybrid structure (see Fig. 1(a)) using the two-mode phenomenological dissipative JC model and the rigorous QNM quantum master equation.First, we show results for the Purcell factor of the quantum emitter.Second, we analyse the system properties, including the Fock distributions and populations.Third, we use the derived expressions of the output electric field operators from Section II and show results for the output photon correlation functions.In Section IV, we will summarize our results and give an outlook to future applications of the theory.We complement the main part of this work with five appendices, that contain a more throughout derivation of the QNM input-output relations, details of the QNM parameters of the hybrid structure, a more detailed derivation of the photon correlation functions, as well as discussions about the response of the hybrid cavity to the external driving, and the treatment of the frequency integrals in the case of a few QNM expansion.

A. Green function quantization approach
We start with the Hamiltonian H = H a + H B + H I , where H a describes a two level system (TLS) with a transition frequency ω a , interacting with a electromagnetic field in an absorptive and spatial-inhomogeneous media using the quantization scheme from Refs.41, 42, and 52.The term H B contains the energy of the medium-assisted electromagnetic field, and H I is a dipole-field interaction Hamiltonian (in the rotating wave approximation).The total Hamiltonian thus contains the following contributions: where b ( †) (r, ω) are annihilation (creation) operators acting on the combined Hilbert space of the dissipative medium and the electromagnetic field degrees of freedom, represented by the spatial index r and frequency index ω.The variable ω must be regarded as continuous mode index, rather then a (temporal) Fourier variable.In fact, in the Heisenberg picture the fundamental operators are b ( †) (r, ω, t), and the time evolution is governed by the Heisenberg equation of motion with respect to H. We further note, that b ( †) (r, ω) fulfills canonical commutation relations.The terms σ ± denote lowering and raising operators of the TLS, describing a point-like emitter with dipole moment d a at the position r a .
The TLS interacts with an effective semi-classical excitation field E L (r a , ω), which reflects a contribution of an incident laser field at the quantum emitter position, that is enhanced by the scattering structures in the dielectric medium, and a medium-assisted quantized electromagnetic field Ê(r a , ω), which obeys the quantized Helmholtz equation where k 0 = ω/c, and (r, ω) = R (r, ω) + i I (r, ω) is the dielectric permittivity.The noise current density ĵN (r, ω) = ω 0 I (r, ω)/π b(r, ω) counteracts the absorption, such that the commutation relations between the electromagnetic field operators is spatially preserved 41,53,54 .A formal solution of Eq. ( 4) is the source- where G(r, r , ω) is the photon Green function, solving the usual Helmholtz equation together with suitable boundary conditions.We emphasize again, that the coupling between the total electric field and the TLS is assumed to be below the so-called ultrastrong coupling regime 55,56 , i.e., |d a , which is consistent with our rotating wave approximation.

B. Quantized quasinormal mode approach
After presenting the quantized Maxwell theory for absorptive systems with a continuous set of modes, we now briefly recapitulate the definition and properties of the QNMs (discrete modes, with complex frequencies), and the Green function expansion.We also discuss the regularized QNMs for expanding the quantized mediumassisted electric field operator (outside the system).
The QNMs for open systems can be viewed in a similar way as the normal modes for closed systems; the QNM vector-valued functions fµ (r) are solutions to the Helmholtz equation however, with open boundary conditions, e.g., the Silver-Müller radiation conditions: for |r| → ∞.Here, (r, ωµ ) is the analytical continuation of the permittivity into the complex plane, and we assume in addition to the lossy media (r, ω) a background region with homogeneous refractive index n B , in which the resonator structure is embedded.As a consequence of the open boundary conditions, the QNM eigenfrequencies ωµ = ω µ − iγ µ are complex numbers with a negative imaginary part, i.e., γ µ > 0. In combination with the fact, that the corresponding (classical) QNM electric fields are harmonic solutions of the wave equation, i.e., Ẽµ (r, t) ∝ fµ (r) exp(−iω µ t), this leads to a lossy character of these modes.Moreover, the QNM eigenfunctions behave in the far field as fµ (r) ∝ exp(in B ωµ |r|/c), and because of the complex eigenvalues (with negative imaginary part) they spatially diverge for far field positions.This leads to an nonphysical behavior 58 of the QNMs outside of the resonator, but is, in fact, a property of any solution to a Helmholtz equation for a lossy geometry with open boundary conditions.In spite of this spatial divergence, when properly normalized 32,34,59,60 , the QNMs can be used to expand the photonic Green function in the form [32][33][34]59,61 for positions r, r nearby (or within) the resonator and where is the QNM Green function expansion coefficient.We remark that there exist alternative forms of A µ (ω), which can be obtained by applying a sum rule of the QNMs 59 , i.e., µ fµ (r) fµ (r ) In this way, one can find an equivalent form of A µ (ω) in G ff (r, r , ω), defined as However, it was shown in Ref.31, that, upon using a Riesz projection technique, the latter form seems to lack a pole contribution at ω = 0, which can lead to an unphysical behaviour of the Green function, when it is not complemented by additional contributions.In contrast, the form in Eq. ( 10) is precisely the expansion coefficient for the total (transverse) Green function, when more general Green functions (including poles at ω = 0) are considered.Therefore, in the following, we will use the arguably more general form in Eq. (10).For positions, R, outside the resonator region, we replace fµ (R) with a regularized QNM Fµ (R, ω) via the Dyson equation 40,62 , where G B (R, r, ω) is the background Green function, solving the Helmholtz equation ( 6) for (r, ω) = B (ω); and ∆ (r, ω) = (r, ω) − B (ω) is the permittivity difference, where B (ω) = n 2 B is the homogeneous permittivity of the background region.
Alternatively, as shown in Ref. 63, one can approximate the expression in Eq. ( 13) with a near-field to farfield transformation using the field equivalence principle 64 : where the terms Jµ S (s ) = n × hµ (s ), (15) Mµ S (s ) = −n × fµ (s ), (16)   are the sources on the boundary, and hµ (r ) = ∇ × fµ (r )/(iω µ µ 0 ) is the magentic field of the associated QNM µ and n is the normal vector on the surface S .In contrast to the regularized QNM in Eq. ( 13), the expression Fµ (R, ω) obtained from the near field to far field transformation, requires the QNM source quantities on a surface S surrounding the resonator region (for details, cf.Ref. 63).
Using the Green function expansion in Eq. ( 9) with the regularized QNMs, obtained either from the Dyson approach or the near field to far field transformation, gives an approximated form of the photonic Green function in terms of QNMs for positions inside and outside the resonator.Using Eq. ( 5) together with Eq. ( 9) and the regularization (either Eq. ( 13) or Eq. ( 14)), we can formulate the total electric field Ê(r) = ∞ 0 dω Ê(r, ω)+H.a. as 40 which explicitly expands in terms of the QNMs.Here, r is the position in the resonator region, and we have introduced the QNM operators: where fµ (r) is replaced by the regularized QNM Fµ (r, ω) for positions outside the resonator.These operators fulfill non-bosonic commutation relations, i.e., [ αµ , α † η ] = S µη , where S µη is a dissipation-induced coupling matrix between QNMs µ, η and has the form describes nonradiative loss processes, i.e., Ohmic heating into the lossy medium, and accounts for the radiative loss, where Hµ (s, ω) = ∇ × Fµ (s, ω)/(iωµ 0 ) is the regularized QNM magnetic field and n points outwards of S. Importantly, we have radiative loss even if the material system is not lossy.Numerically, we have verified that the radiative part of S has to be chosen in the far field (i.e., at least half a wavelength away from the resonator), in order to get a convergent value of S rad µµ .This is clear, since otherwise one gets near-field evanescent contributions that can even be negative, and these are not associated with far field propagation decay.Therefore, by choosing S = S ∞ as a far field surface and applying Silver-Müller radiation condition, we obtain the form After applying a symmetrization orthogonalization transformation 40 , we can rewrite the electric field operator from Eq. ( 17) in a symmetrized QNM basis, with the symmetrized QNM functions, 2 ) νµ ω ν /ω µ fν (r), (24)   and a µ = η (S −1/2 ) µη αη and a † µ are the annihilation and creation operators for the symmetrized QNMs.As a consequence of the symmetrization 40 , these operators (a µ and a † µ ) obey bosonic commutation relations and can be used to describe Fock number states of mixed photon/lossy medium excitations.These photon number states are not eigenstates of the photon Hamilton operator, namely Eq. ( 2), since the photon number operators do not commute with the Hamiltonian H B .

C. Multi quasinormal-mode master equation
Next, we provide a more detailed derivation of QNM master equation, introduced in Ref. 40, also taking into account the external laser field.Exploiting the general quantized QNM theory developed in the last subsection, we can determine the time evolution of the QNM annihilation operator, a µ , with respect to the Hamiltonian H in Eqs.(1)- (3).Using the Heisenberg equations of motion, we derive within the Markov approximation (cf.App.A) where we have defined H sys = H em +H a +H em-a +H L as the effective system Hamiltonian in the symmetrized QNM basis, with µη a † µ a η and H em-a = µ gs µ a µ σ + + H.a., and we have introduced the symmetrized QNM-TLS coupling constant gs µ = ω µ /(2 0 )d a • f s µ (r a ) and the coherent mode-coupling term χ 2 ) νη .The pumping term, H L , accounts for the interaction of the dipole with the external laser field (cf.Eq. ( 3)).
We next derive the QNM master equation based on the quantum Langevin equation in Eq. (25).We first assume that the incident laser field reflects a cw excitation with a detuning with respect to the TLS frequency ω a , i.e., such that the effective classical scattered field, i.e., the incident laser field enhanced by the cavity structure, can be approximated as E L (r a , ω, t) ≈ F L (r a )e −iωLt .We subsequently rewrite the dipole-laser interaction Hamiltonian from H I of Eq. ( 3) as where is the Rabi frequency of the enhanced classical laser field.It should be noted, that although we choose here an effective driving of the quantum emitter, there is an equivalence of cavity pumping and quantum emitter pumping, in the sense that the cavity driving leads to an effective quantum emitter driving via the cavity-exciton interaction 66 .Furthermore, we treat the input operators as white noise, and assume that the corresponding input state is the vacuum state, i.e., there are initially zero quanta in the input states, such that a in µ (t)a in † η (t ) = δ µη δ(t−t ) and all other second order correlation functions vanish 67 .In addition, we assume that the eigenfrequencies of Hsys are not degenerate 67,68 .
Applying the Ito-Stratonovich calculus 67 to the Heisenberg equation of motions of the symmetrized QNM operators, and using the procedure from Ref. 40, we obtain the master equation for the symmetrized QNMs and TLS (in a rotating frame with respect to ω L ): where H sys = H a + H em + H L + H em−a is the effective system Hamiltonian (described after Eq. ( 25)) with , and H L = Ω(σ + + σ − ) in the rotating frame with the (laser) detuned TLS frequency ∆ a = ω a − ω L .We stress again, that due to inter-mode coupling terms, the photon number operators N µ = a † µ a µ do not commute with the Hamiltonian, as would be usually the case in a Fock space without dissipation.The QNM Lindblad dissipator is derived as also yielding an off-diagonal coupling via the decay matrix χ (−) µη .We further added the Lindblad dissipator, with the (background) spontaneous emission (SE) rate which accounts for non-cavity decay of the TLS.We highlight that all mode related coupling parameters entering the above QNM master equation are directly obtained from the QNM calculations and quantum emitter properties without any form of phenomenological fitting.We also note that, for convenience, we will refer to the quantized QNM model as QNM-JC model, since the QNM master equation from Eq. ( 27) can be viewed as a generalized and rigorous dissipative JC model.

D. Diagonalization of the Lindblad dissipator
In this subsection, we apply a unitary transformation to the QNM master equation, Eq. ( 27), to diagonalize the decay matrix χ (−) (defined below Eq. ( 25)).This will support the discussion of Section III, and will make the role of off-diagonal QNM coupling more clear, as it will be entirely encoded in the Hamiltonian part of the master equation.Since χ (−) is a semi-positive definite and Hermitian matrix, there exists a unitary transformation U (−) , that diagonalizes χ (−) , such that ν,ν where Γ µ are the eigenvalues of χ (−) .In the new basis, the Lindblad dissipator L em takes the diagonal form where µ are QNM annihilation (creation) operator in the diagonalized dissipator frame with µη is unitary, the bosonic commutation relations of the QMN operators are preserved.The effective system Hamiltonian in the diagonalized basis reads with ∆ µ = Ω µ −ω L , and the coupling constants transform as In Eq. ( 35), we have defined g µη em ≡ χ(+) µη for µ = η as the photon-photon coupling constant between QNM µ and η, and Ω µ ≡ χ(+) µµ as the bare mode frequencies in the diagonalized and symmetrized picture.We note that χ(+) has the same eigenvalues as χ (+) , i.e., the eigenenergies of the full photon Hamiltonian are not changed by the unitary transformation.Furthermore, the photon number operators N µ = A † µ A µ do also not commute with the Hamiltonian in the diagonalized frame.

E. Input-output relations and output electric field operator
Here, we derive the output electric field operator in the far field region and the input-output relations for multiple QNM operators (in the Heisenberg picture).First, we write down the time-reversed quantum Langevin equation of Eq. ( 25): where a out µ = ∞ 0 dωa out µ (ω) is the output operator, which is explicitly given in App.A, Eq. (A27).Next, we subtract Eq. ( 38) from Eq. ( 25) and multiply from the left with χ (−) −1/2 to obtain the input-output relations, We remark that as a consequence of the dissipationinduced coupling, the output and input of a symmetrized QNM µ is related to a linear combination of all QNM operators, which is in contrast to the standard input-output relation 67 , where the input and output channel is connected via a single system operator.In the diagonalized basis, using A µ instead of a µ , we obtain the (diagonalized) input-output relations: One can also formulate input-output relations in ω-space of the forms where, obviously, one obtains the ω-independent rela-tions by integrating over all ω on both sides, respectively.These input-output relations will be used in the following to formulate the output electric field operator, which is important for simulations involving correlation functions at a outside detector.
We start with the full (positive-rotating) electric field operator Ê(+) (R, t) = ∞ 0 dω Ê(R, ω, t) at a position R outside the resonator from the source-field expression from Eq. ( 5) using the QNM Green function together with the field regularization (Eq.( 13),( 14)): where is a regularized QNM (Eq.( 13),( 14)) in the symmetrized basis, and Subsequently, we rewrite Ê(+) (R, t) from Eq. ( 43) as with Now we can use the input-output relations from Eq. ( 41), to obtain the representation Ê are the (cavity) output and input electric field operators.Equation ( 48) yields a general expression for the QNM output/input field for the multi-QNM case.In the following, we concentrate on positions R in the far field, i.e. |R| max(λ µ ), to obtain an approximated form of the output/input fields for numerical calculations in Section III, that connects to the ω-independent QNM system operators a µ .Since Fµ (R, ω) is proportional to the background Green function G B (R, r, ω), where r is either located in the resonator volume (Eq.( 13)) or at the resonator boundary (Eq.( 14)), i.e., |R| |r|, we can approximate the regularized function as where Zµ (R, ω) for the Ansatz in Eq. ( 14), is given as and R = R/|R| is the unit vector in the direction of R. Inserting Eq. ( 49) into the output/input electric field operator from Eq. ( 48), yields where Z s µ (R, ω) is implicitly defined via Eq.( 47) together with Eq. ( 49) and (50).Using the definition of a out/in µ (ω, t) (Eq.(A23) and Eq.(A27)), it follows that ) varies slowly with respect to ω around the QNM frequency ω µ , we apply a resonance approximation to obtain the final expression for the output field operator: where we used again a out/in µ = ∞ 0 dωa out µ (ω).In the diagonalized basis, the output/input field reads with Using the ω-independent input-output relations from Eq. (39) or Eq. ( 40), we can connect the output field to the system QNM operators a µ or A µ .In particular, the output electric field operator at position R and time t is then a linear combination of far-field regularized QNMs Z s µ (R)(Z sU µ (R)) and QNM system operators a µ (A µ ) as well as QNM input operators a in µ (A in µ ) at time t−n B |R|/c.The introduction of the above output electric field operators (Eq.( 52) and ( 53)) allows one to calculate, e.g., second-order photon correlation functions.

III. APPLICATIONS TO COUPLED OPEN RESONATORS
In this Section, we will apply the theory from Section II to a two-QNM, multiphoton system, using first principle calculations for a specific open cavity structure.A typi-cal system to study in terms of two dominant but different QNMs are metal-dielectric hybrid structures, where one mode is photon dominated and one is plasmon dominated 21,47,69 , but with a sufficiently different quality factor.
We will focus on the hybrid structure depicted in Fig. 1, which shows two fundamental QNMs in the optical frequency regime.In particular, we will discuss differences between a phenomenological dissipative JC model and the QNM-JC model with respect to density matrix equations results of the hybrid structure.In subsection III A, we will discuss the Hamiltonian and the dissipator of the master equations on a formal basis.Afterwards, both models (phenomenological dissipative JC and QNM-JC) will be compared in the weak photonemitter coupling regime in subsection III B, and subsequently in the intermediate emitter-photon coupling regime in subsections III C and III D.
A. Two-mode master equations, hybrid cavity and TLS parameters We start with the formulation of the master equation for the hybrid structure, e.g., as shown in Fig. 1.All input QNM and TLS parameters used for numerical evaluation can be found on the fourth row of Tab.I.A more detailed description of the QNM input parameters is given in App.B.
For the two-mode hybrid case, we rewrite the Hamiltonian for the QNM-JC model H QNM sys ≡ H sys from Eq. ( 35) with µ = {pc, pl} as where pl is the annihilation (creation) operator for the plasmon-like mode and A ( †) pc is the annihilation (creation) operator for the PC-like mode in the diagonalized dissipator picture, as illustrated in Fig. 1 (c-d).
In this basis, ∆ pl(pc) = Ω pl(pc) − ω L is the detuned plasmon-like (PC-like) frequency; g pl and g pc are the coupling constants of the plasmon and PC mode to the TLS, respectively.To simplify the notation, we have set g pl,pc em = g em as the plasmon-PC mode coupling constant.The Lindblad dissipator is then with and where Γ pl(pc) is the diagonalized decay rate of the PC-like (plasmon-like) QNM, defined implicitly from Eq. (32).In contrast, a phenomenological dissipative JC model assuming [α µ , α † η ] = δ µη for µ, η = 1, 2 is represented by with gpl(pc) = −i ω pl(pc) /(2 0 )d a • fpl(pc) (r a ) using the untransformed QNM fields of the hybrid and ∆pl(pc) = ω pl(pc) − ω L .The Lindblad dissipator reads and we note that ω µ − iγ µ (µ = {pl, pc}) are the real and imaginary part of the original, individual hybrid QNM eigenfrequencies, respectively.We emphasize that the Lindblad dissipator L SE ρ associated to the vacuum spontaneous emission rate γ SE of the TLS is the same in both models.
The corresponding master equations are given by or for the QNM-JC model and phenomenological dissipative JC model, respectively.Next, to specify the differences between both master equations, namely Eqs.(60,61), we compare the occurring coupling parameters in the two different master equations for the same hybrid structure, which are summarized in Tab.I. We see that the overall behaviour of the dominating plasmon-mode related quantities are very similar before (phenomenological dissipative JC model) and after symmetrization and diagonalization (QNM-JC model).However, some of the PC-mode related parameters change drastically, as described below: (i) Dissipation.-Theeffective width Γ pc of the symmetrized PC-QNM is around one order of magnitude smaller compared to the original PC width γ pc .This is a consequence of the structure of the decay matrix χ (−) (cf. the text surrounding Eq. ( 25)), since its elements are linear combination of the original complex eigenfre-quencies ωµ .Indeed, for two modes, the exact decay eigenvalues Γ pl/pc for the plasmon mode and the PC mode are and, using the properties of the trace and the determinant, Here, R is defined through which is a dissipation-induced correction factor to the initial phenomenological damping.
The analytic form of the eigenvalues in Eq. ( 63) has an interesting implication: if γ pl and γ pc are very different from each other, e.g., γ pl γ pc , then the plasmonrelated eigenvalue Γ pl is only slightly shifted compared to γ pl , since the correction by the term involving γ pc is very small.In contrast, Γ pc is mainly influenced by the correction term corresponding to γ pl .On the other hand, when γ pc ∼ γ pl ≡ γ, then we find, as a first estimate (for a very small difference below 10%), Thus there is a symmetric splitting of γ, which depends on the detuning of both modes and the mode overlap |S pl,pc |.Strictly speaking, Eq. ( 65) is exact for degenerate QNM imaginary parts γ pl = γ pc , which is a technical interesting case, since it appears, e.g., in Fabry-Prot cavities.
In the example of Fig. 1, the former case γ pl γ pc applies, which explains why the effective PC-mode decay rate Γ pc is significantly shifted from γ pc (cf.Tab.I).It should be noted that these rate changes are also present in the full eigenvalues of the Liouvillians, which will be important for the response of the hybrid to an external optical field, as we will discuss later.
(ii) QNM-TLS coupling.-Another interesting observation is that the PC-TLS coupling constant g pc in the QNM JC-model (Eq.( 55)) is also nearly one order magnitude lower compared to the phenomenological dissipative JC parameter gpc (Eq.( 58)), using the original PC-mode eigenfunction fµ ; this is because g pc is formed by a linear combination of gpc and gpl , which also deviate by one order of magnitude to each other.Interestingly, since the effective PC decay rate changes by a similar amount, both parameter sets indicate the intermediate coupling regime of the photons and emitter, i.e., |g pc |/(2Γ pc ) ≈ |g pc |/(2γ pc ) ≈ 4 in the QNM-JC model as well as in the phenomenological dissipative JC model (cf.Eq. ( 55) and ( 58)).
To summarize this subsection, there are two key changes that occur by comparing the two master equations, Eq. ( 60), (61).We first recognize modifications of the mode and coupling parameters (mainly PC mode) due to symmetrization and diagonalization of the QNM annihilation and creation operators, and second, there are additional contributions in the Hamiltonian due to the presence of the photon-photon interaction part with the coupling constant g em (Eq.( 55)).The impact of these key changes, which results from a proper treatment of the QNM quantization (with real losses), will be shown below by explicitly comparing the phenomenological dissipative JC and QNM master equation simulations.First, in subsection III B, we inspect the weak light-electron coupling regime, where we apply the bad cavity limit as in Ref. 40 (but for a completely different hybrid structure), and adiabatically eliminate the cavity modes from the master equations, Eqs. ( 60), (61).In this limit, we compare the cavity-enhanced spontaneous emission of the quantum emitter obtained from the QNM-JC model, the phenomenological dissipative JCmodel and an independent semi-classical solution.Subsequently, in subsection III C and III D, we analyse the multiphoton regime in the intermediate to strong lightelectron coupling regime, where the bad cavity approximations are not valid anymore.Thus, we will use the full master equations, Eqs. ( 60), (61).The numerical results of the master equations were calculated using the library Quantum Toolbox in Python 70 (QuTiP) and we note again, that all parameters, which enter the master equations are summarized in Tab.I.

B. Weak light-exciton coupling regime: Purcell factors and radiative β factors
In a first step, we compare the QNM-JC model and the phenomenological dissipative JC model in the weak coupling limit, i.e., we reduce the QNM-emitter constants by choosing a small dipole moment of the TLS and leaving the mode parameters unchanged.In the weak coupling limit, we adiabatically eliminate both modes from the master equations, Eq. (60-61) (as in Ref. 40) to obtain the master equations in the bad cavity limit 71 , and for the QNM-JC model and phenomenological dissipative JC model, respectively.Here, Γ and Γ are the cavity enhanced spontaneous emission rates of the TLS, defined values are linearly interpolated.Also note that the maximum Purcell factor of the 2 nm gap hybrid is one order of magnitude higher compared to the 5 nm gap hybrid.The value of β rad class for the 2 nm gap has additional peaks for frequencies towards the low-Q mode, which likely come from non-modal quasi-static coupling and a constant background term with higher order modes, whose influence reduces with decreasing gap size.For all simulations, the classical dipole is at the gap center of the dimer and is z-polarized.
where ∆ µa = ω µ − ω a is the QNM-TLS detuning.We stress again that both models use the same original QNM parameter (ω µ , γ µ , gµ ).From the above equations, we clearly see, that in the limit S → 1, i.e., when there is vanishing radiative and non-radiative QNM overlap (cf.Eq. ( 19)), Γ and Γ coincide.However, we should note here, that this is really only the case, when γ µ → 0, which in a sense contradicts with the assumptions of a phenomenological dissipative JC model for finite loss.
Comparing the QNM-JC model and the phenomenological dissipative JC model with respect to the bad cavity limit master equations (Eqs.( 67),( 69)), we see that the differences can be summarized as additional offdiagonal terms in the QNM-JC cavity-enhanced spontaneous emission rate Γ of the TLS, Eq. ( 70).In contrast, in the full master equations (Eqs.( 60),( 61)), the differences of both models is not only present in the coupling constants, but also off-diagonal coupling between the different mode operator appear in the QNM-JC model.
For a demonstration of the influence of the QNM coupling terms in Eq. ( 70), we calculate the Purcell factor F P = Γ/γ SE and FP = Γ/γ SE (cf.Eq. ( 31)) as a function of the TLS frequency assuming the approximative bad cavity limit 40 , for the QNM master equation and the two-mode phenomenological dissipative JC master equation.To estimate the quality of the Purcell factor results and the underlying models in the bad cavity limit, we compare these results with a (independent) semi-classical Maxwell simulation.
However, before comparing the two different quantum models and semi-classical model for the specific hybrid in Fig. 1, we first discuss the choice for the design of the hybrid structure using results obtained solely from the full Maxwell simulations.To do so, we compare the hybrid design from Fig. 1 with a gap of ellipsoidal dimer of 2 nm and 5 nm with respect to the Purcell factor and radiative β-factor in Fig. 2.There are two main differences between both cases: First, as shown in Figure 2 (first and second row), the maximum Purcell factor of the hybrid structure with 2 nm gap is roughly one order of magnitude higher compared to the 5 nm gap hybrid.Second, the β-factor (cf.Fig. 2 (third and fourth row)) shows a constantly increasing contribution with additional resonances for frequencies towards the low-Q hybrid resonance for the 2 nm gap case, which are caused by effects beyond the two hybrid resonance description.Therefore, we have chosen an extreme case for the hybrid design for the quantum simulations, in that the gap between the two metallic ellipsoids is only 2 nm (cf.Fig 1).In such a small gap regime, one expects additional effects beyond the main two modes, and also the numerical calculations of the total beta factor are more difficult (cf.App.B).We choose such an extreme case, since it is more interesting for cavity-QED and emerging experiments 20,72 , since it allows one to reach the intermediate to strong coupling regime for realistic dipole strengths of the TLS, because of the large Purcell factor.
After having discussed the choice for the hybrid design in Fig. 1 tion.As shown in Fig. 3 the QNM master equation result is in very good agreement with the full but semiclassical Maxwell simulations and reproduces the pronounced Fano effect near the original PC-like mode frequency ω pc .Notably, the peak of the Purcell enhancement of the PC-like cavity mode is located at the eigenvalue (pc) em of the electromagnetic part of the Hamiltonian H em , in agreement with the derivation for the hybrid structure in Ref. 40.While the Purcell factor obtained from the phenomenological dissipative JC-model is in good agreement with the full Maxwell solution near the plasmonic-like mode frequency ω pl , it fails in the frequency regime where interference occurs.This is caused by the missing inter-mode coupling terms in the phenomenological dissipative JC model 40 (Eq.( 71)), which is present in the QNM-JC model.It should be noted, that S 22 was decreased from 1 to 0.77 in the phenomenological model to at least match the height of the phenomenological dissipative JC result at the PC peak.This leads to a slightly modified PC-TLS coupling constant for the phenomenological dissipative JC model of g2 → √ 0.77g 2 , but all other parameters are identical.72)) and the phenomenological dissipative JC model (green dashed line with β ≈ 0.14).Note, that for the JC-model, a phenomenological introduced constant is shown, which mimics the behaviour of the radiative output at TLS frequencies around ωpc with respect to the classical β rad class solution from Fig. 2. The inset shows a zoom-in around the PC-like mode frequency ωpc = 1.6052 (eV).In the QNM-JC model, a peak close to the PC-like mode with frequency ωpc appears, which is located at the position of the highly non-Lorentz interference effect (dip in Fp) from Fig. 3.
Next, we discuss discuss the difference between the QNM-JC model and the phenomenological dissipative JC model in terms of (modal) β-factors.For the QNM-JC model, we define the β-factor as β rad QNM = Γ rad /Γ with As one can see from Fig. 4, β rad QNM reproduces the overall shape of the full Maxwell solution from Fig. 2 (third row) through the presence of S rad µη and S nrad µη .In contrast, in a phenomenological dissipative JC model, it is not clear at all how to separate radiative and non-radiative contribution, and this model is not even able to predict a frequency-dependent modal β-factor: While adding a constant β-factor phenomenologically can yield appropriate results in the single-mode limit, this is non-trivial for more then one mode, as is shown by the (green dashed) line in Fig. 4, since off-diagonal effects between the modes can alter the output behaviour.

C. Intermediate light-exciton coupling: quantized system properties
Next, we turn to the case of the intermediate QNMemitter coupling regime; this regime is realized by the parameters shown in Tab.I (Γ pc/pl ∼ g pc/pl and γ pc/pl ∼ gpc/pl ), where we use the full master equation (Eq.(60-61)) beyond the bad cavity limit (Γ pc/pl g pc,pl and γ pc/pl gpc,pl ).For the following calculations, we choose the TLS frequency resonant to the PC-like eigenfrequency of both photon Hamiltonians (H QNM em and HJC em ), i.e., ω a = (pc) em .In the phenomenological dissipative JC model, the PC-like eigenfrequency of HJC em is simply the PC mode frequency itself, i.e.
(pc) em = ω pc .In contrast, in the QNM-JC model, the PC-related eigenfrequency of H QNM em is slightly red shifted compared to ω pc (cf.Fig 3)  and which we chose as the TLS frequency in the QNM-JC model.

Eigenenergies and eigenstates
To provide the essential background for the interpretation of our full master equation results, we first discuss the eigenenergies and eigenstates of the coupled TLS-QNM system with respect to the QNM-JC model and the phenomenological, dissipative JC model without an external pump, i.e., we set Ω = 0 in the Hamiltonians, Eq. ( 55) and (58).Note that the respective Hamiltonians without external pump are denoted as H QNM sys , H JC sys .To obtain a formal inside into the eigenstates of H QNM sys >@ >@ >@ >@ >@ >@ >@ >@ . Absolute square of the basis coefficients c i N kpc,k pl ,l for the one-and two-excitation manifold (N = 1, 2) with the TLS frequency ωa = Since the total number operator commutes with the system Hamiltonians without external pumping, i.e.
[ N , H QNM sys ] = 0, there exists a common eigenbasis of N and H QNM sys .However, since the degeneracy of the eigenvalues of H QNM sys is generally different compared to N , the eigenbasis formed by the eigenspaces from Eq. ( 73) is not necessarily a eigenbasis of H QNM sys .A common eigenbasis with eigenstates |ϕ N,i N , where i N = 1, . . .2N + 1, can be defined as a linear combination of basis elements |φ N,i N of the eigenspaces from Eq. ( 73) with respect to i N .Without loss of generality, we choose |k pl , k pc , l; j N δ N,kpc+k pl +l with j N = 1, . . .2N + 1 as the basis set of {|φ N } to construct the states for the N -th rung in the QNM-JC ladder (and formally equal for the phenomenological dissipative JC model).The states |ϕ N,i N are solutions to the eigenvalue problem where E N,i N are the (real) eigenenergies of the system Hamiltonian H QNM sys .Note that we order i N for a specific manifold N , such that E N,i N < E N,i N +1 for all i N .Although |ϕ N,i N are not eigenstates of the full Liouvillians (Eq.( 60)), they still reflect the effect of coupling between the subsystems and, in contast to the full eigenstates, constitute a orthonormal basis, which will later be used for the definition of a projection operator.
Since the parameters of the hybrid system (cf.Tab.I) indicate a regime well below the strong/ultrastrong 55,56 (cf.subsection II A) coupling regime (g pc/pl ω pl,pc ), each |ϕ N,i N is approximately dominated by a single bare state |n pl , n pc , l , as shown for N = 1 and N = 2 in Fig. 5.It is therefore instructive to rename the eigenstates corresponding to the one-excitation manifold (N = 1) as |ϕ Note, that in both subfigures, there is a sudden jump on the energy axis from the first to second rung and that the scaling is different in the upper and lower half.
The energies E N,i N are not eigenenergies of the full Liouvillians (from Eq. ( 60),( 61)), since the additional dissipative terms add imaginary parts to E N,i N and shifts the real parts due to coupling of the subsystems in the coherent part.However, the complex eigenenergies N,i N of the non-hermitian Hamiltonian  60)), excluding the transition energies 73 .This holds also true in the case of the phenomenological dissipative JC model, where, obviously, H QNM sys → HJC sys and Γ pl(pc) → γ pl(pc) in Eq. ( 76).The complex eigenenergies N,i N for the first and second rung are depicted in Fig. 6.We note, that the real part of these energies are very similar in the QNM-JC model and phenomenological dissipative JC model, although the full energies for the latter are not shown.However, the imaginary part of the eigenenergies, corresponding to the states dominated by the PC contributions, are very different (cf.Fig. 6, b), as is the case for the bare rates Γ pc and γ pc (cf.Tab.I).We em-phasize that this is a consequence of the symmetrization of the QNM operators in the QNM-JC model (necessary to construct Fock states), which yields symmetrized mode parameters as linear combination of the input QNM parameters, used by the phenomenological dissipative JC model.Furthermore, the increase of the imaginary part in higher rungs (N > 1) is also different in both models: While in the QNM-JC model, we calculate Im( pc−pc ) ≈ 3Im( pc ), for the phenomenological dissipative JC model, we find Im( pc−pc ) ≈ 2Im( pc ), which can lead to major differences between both models with respect to the response to an external laser in the higher rungs of the energy ladder.
We note that, in the following, we will also adopt the notation introduced above for the eigenstates of H QNM sys (H JC sys ) for N = 1, 2, to the complex eigenenergies, e.g., 1,1 = a .

Steady-state probabilities and occupation numbers
Having discussed the eigenenergies and eigenstates of the QNM-JC model and the phenomenological dissipative JC model without pump, we now apply an external optical driving on the system, and simulate the full master equations (Eq.(60-61)).In this situation, via external pumping, few photon effects can be studied.We choose a Rabi frequency of Ω L = 0.025|g pl | ∼ |g pc |, which is in a excitation regime, where effects from the two-excitation manifold of the JC ladder are visible (cf.App.D for discussion on excitation regimes) and we note, that gpl is the TLS-plasmon coupling constant in the original QNM basis using fpl .
Evaluating Eqs. ( 60) and ( 61) numerically, we now analyse the (total) probability to find the system in a 0, 1 or 2-exitation (photon) state as well as the occupation numbers of the two QNMs and the TLS for t → ∞ (steady-state regime) as a function of laser frequency ω L , so that we access the intrinsic quantum anharmonicities of the higher rungs of the JC ladder.The quantities of interest are the occupation numbers n pl ≡ A † pl A pl , n pc ≡ A † pc A pc and n a ≡ σ + σ − for the plasmon-like mode, PC-like mode and TLS upper level |e , respectively, as well as the probabilities on the i N -th eigenstate |ϕ N,i N of H QNM sys (or H JC sys for the phenomenological dissipative JC model).For N = 1, 2, we also adopt the notation introduced in the last subsection for the probabilities, e.g., P [1,1] = P a .
(i) One-excitation manifold.-First,we study the 1excitation manifold as well as the occupation numbers and choose a laser frequency regime around the PC mode frequency ω pc , where the most striking differences between the QNM-JC and the phenomenological dissipative JC model are visible.As explained in the last subsec-tion, the 0-excitation manifold only contains the trivial vacuum state |ϕ vac and the 1-excitation manifold yields the three eigenstates |ϕ a , |ϕ pc and |ϕ pl .
em is aligned to the eigenenergy of the full electromagnetic Hamiltonian Hem (ωpc in the phenomenological dissipative JC model) and the TLS is pumped with an external laser with Rabi frequency Ω ∼ gpc, where gpc ∼ Γpc is the TLS-PC coupling constant (in the intermediate photon-exciton coupling regime) and Γpc is the PC mode decay rate (cf.Tab.I).Additionally, we show Ppc for the QNM-JC model in the limit gem = 0 (dotted line), cf.Eq. ( 55).
The results are shown in Fig. 7, and we start with the analysis of P a and n a , corresponding to the probability and occupation of the TLS-like upper level, respectively.The peak heights of both quantities differ strongly in the two different models: Whereas in the phenomenological dissipative JC model, n a < 10 −2 , in the QNM-JC model n a ≈ 0.1, i.e. one order of magnitude larger.This is because of the stronger coupling (|g pc | ≈ 10|g pc |) between the TLS and the PC like-mode and the larger depopulation rates (γ pc ≈ 10Γ pc ) of the PC-like states in the phenomenological dissipative JC model.Quantitatively, this difference is similar for P a ; however, P a itself is much smaller than n a , which is consequence of |ϕ a being a linear combination of the bare states |1, 0, g , |0, 1, g and |0, 0, e , which leads to a reduction of the diagonal con-tribution |0, 0, e 0, 0, e| (cf.Fig. 5) in P a .
Next we look at P pl and n pl , corresponding to the probability and occupation number of the plasmon-like mode, respectively.The probability P pl has a negligible value, P pl < 10 −3 , for all values of ω L that we tried, which is the case for both models.Furthermore, the occupation number n pl of the plasmon-like mode is very small (< 10 −2 ) and also very similar in both models.The similarity for the plasmon-like mode is a consequence of the small differences in the plasmon decay rate (Γ pl ≈ γ pl ) and plasmon-TLS coupling constant (g pl ≈ gpl ) before and after symmetrization and diagonalization (see also Tab.I).The small values are a consequence of the high plasmon decay rates Γ pl , γ pl , which results in a fast depopulation of the plasmon-like states.
The most pronounced changes appear in P pc and n pc of the PC-like mode.First, the peak of P pc in the QNM-JC model is slightly detuned to higher frequencies compared to Ppc of the phenomenological dissipative JC model.This is a consequence of the small shift of the PC-like frequency (ω pc ) after symmetrization and diagonalization (Ω pc ).Second, the spectral width of the probability and occupation number dynamics is much broader in the case of the phenomenological dissipative JC model.This is again a consequence of the width γ pc of the PC-related resonance of the phenomenological dissipative JC model being about one order of magnitude broader compared to the effective PC width Γ pc of the QNM-JC model (γ pc ≈ 10Γ pc ).Therefore, the laser can effectively excite a much broader range of the first rung in the phenomenological dissipative JC model ladder.Third, the peak height of P pc and n pc between both models is completely different; P pc is about one order of magnitude larger than Ppc , which means that the system has a 10 times larger probability to be in the state |ϕ pc in the QNM-JC model.Of course, this is because of the fast depopulation of the first rung in the phenomenological dissipative JC model due to the large decay rates.Due to the same reason, there is also a pronounced difference in the peak height of n pc .
There are two further interesting observations that are connected to n pc and ñpc .First, in the phenomenological dissipative JC model, ñpc is nearly identical to Ppc close to ω pc , as shown in Fig. 7 (bottom, red and green dashed curve).To explain this, we recall that and in the phenomenological dissipative JC model, we observe that n pc ≈ |0, 1, g 0, 1, g| ≈ P pc .In contrast, in the QNM JC-model, n pc is different to P pc with respect to height and peak position, indicating that higher rung (N > 1) probabilities are also important here, as we will show below.Second, an additional indicator for processes on higher rungs in the QNM-JC model is a start of a spectral hole burning process at the peak of P pc , which comes from higher photon probabilities with a smaller laser excitation width. is aligned to the eigenenergy of the full electromagentic Hamiltonian Hem (ωpc in the phenomenological dissipative JC model) and the TLS is pumped with an external laser with Rabi frequency Ω ∼ gpc, where gpc ∼ Γpc is the TLS-PC coupling constant (in the intermediate photonexciton coupling regime) and Γpc is the PC mode decay rate (cf.Tab.I).We also show Ppc−pc for the QNM-JC model in the limit gem = 0 (dotted line), cf.Eq. ( 55), and the occupation number npc normalized to the same height as 0.1Ppc−pc.Note, that all other contributions corresponding to the twoexcitation manifold are negligible and not shown here.
(ii) Two-excitation manifold.-Next,we discuss the probabilities as a function of laser frequency around the PC mode frequency, connected to the 2-excitation manifold of the QNM-JC and phenomenological dissipative JC ladder, which consists of the five eigenstates: |ϕ a−pc , |ϕ pc−pc , |ϕ a−pl , |ϕ pc−pl and |ϕ pl−pl .The dominant probabilities with peak heights P > 0.001 are plotted in Fig. 8 (P a−pc , P pc−pc , P a−pl ).All other contributions (P pl−pc , P pl−pl ) are not shown.We notice, that quantitative differences of P a−pc , P a−pl between both models behave similar to the differences in P a from the 1-excitation manifold.However, in both models, these higher-rung contributions are negligible compared to the 1-excitation manifold contributions.Thus, below we concentrate on P pc−pc , which covers the most striking and interesting differences.
We observe that P pc−pc is roughly two order of magnitude higher compared to Ppc−pc , which constitutes a much larger difference compared to the case of the oneexcitation manifold.This leads to a different ratio of the maximum values R 1−2 = max (P 1,i ) /max (P 2,i ) between the one-and two-excitation regime: In the QNM-JC model, we obtain R 1−2 ≈ 2, while in the phenomenological dissipative JC model, we obtain R 1−2 ≈ 10 2 .To help explain this, we show the two photon resonance In the QNM-JC model, ω L at the peak of P PC−PC is not in the range of the PC-like eigenenergy Re( pc )±Im( pc ).However, the two photon resonance (2ω L ) is located in the range of the eigenenergy Re( pc−pc )±Im( pc−pc ) corresponding to PC-like eigenenergy on the 2nd rung.This means, that there is a direct population of the 2nd rung in the QNM-JC model via a virtual state with smaller energy compared to Re( pc ) − Im( pc ), leading to the relative high 2-photon probability.This is possible due to the anharmonicity of the JC ladder, i.e., the transitions Re( N +1,i ) − Re( N,i ) between different manifolds N and N + 1 depends on N due to the emitter-photon coupling.Of course, this anharmonicity is present in both models, the QNM-JC model and the phenomenological dissipative JC model, but takes a different value in the QNM-JC model as a direct consequence of H em not commuting with the photon number operators.However, in the phenomenological model, the widths of the eigenenergies in the one-excitation manifold are much larger than the difference D pc ≡ |Re( pc−pc ) − 2Re( pc )| of the transitions, and thus, the first rung is already majorly excited at ω L > ω L and the probability of an indirect population of the second rung via a virtual state is very small (cf.Fig. 6, b).We note, that this effect in the QNM-JC model is stable against a variation of ω a with respect to the PC-mode frequency.Additionally, it is worth to note, that the peak position of the occupation number n pc is nearly directly located at the peak of P pc−pc (cf.Fig. 8).
We briefly summarize the analysis in this subsection: For moderate pumping with respect to the PC-mode, Ω ∼ g pc , the QNM-JC model exhibits a relatively high probability to be in a 1 or 2 excitation state for the PClike mode (P pc ∼ 0.4 and P pc−pc ∼ 0.2), while the vacuum state |ϕ vac is surpressed, when the laser is tuned in the regime of the PC-like cavity mode frequency.In contrast, in the phenomenological dissipative JC model, the overall behavior of the resonance structure is similar, but qualitative and quantitative differences are present: The 2-photon probability in the phenomenological dissipative JC model is negligible over the inspected laser frequency regime, caused by the larger decay rate γ pc of the PC-like mode.In particular, the maximum peak of the 2-photon probability P pc−pc close to the PC-like mode frequency, is about two orders of magnitude smaller compared to the peak of the QNM-JC model (cf.Fig. 8).This significant difference between both models is caused by the presence of the inter-mode coupling in the Hamiltonian, Eq. ( 55), (the case g em → 0 is indicated by the dotted line in Fig. 7 and 8 (bottom)) and the shift of the PC decay rate (cf.Eq. ( 63)), but not caused by the TLS-QNM coupling renormalization: While in the phenomenological dissipative JC model the coupling constant gpc of the TLS to the PC-like cavity mode is about one order of magnitude larger (cf.Tab.I), the decay rate γ pc of the PC-like mode is also one order of magnitude larger, leading to the same ratios of both quantities, i.e. |g µ |/(2Γ µ ) ∼ |g µ |/(2γ µ ), as discussed earlier.
Therefore, both models are in the intermediate lightmatter coupling regime, and the difference in the response of the system to the external pump is mainly influenced by the different PC decay rates and the intermode coupling Hamiltonian with coupling g em .Thus, while the system described on the basis of the phenomenological dissipative JC model is (mainly) in the single-photon regime, the system described by the QNM-JC model has significant multiphoton properties.This shows, that the phenomenological introduction of a two-mode parameter set is not trivial in any way, since the input parameter drastically change because of dissipation itself in the course of deriving the QNM master equation, since decay rates and coupling strength are not independent from each other due to symmetrization (S) and diagonalization (U (−) ).Therefore phenomenological models can clearly miss important features in the multiphoton regime.
It is important to stress that our QNM quantization model can be used to predict new regimes in dissipative quantum optics, since all parameters are calculated on a solid foundation through the QNM eigenfunctions and eigenvalues without any phenomenological approaches.

D. Intermediate light-electron coupling: Output properties
Next, we analyse the impact of the off-diagonal QNM coupling on experimental observables, represented by the derived output electric fields, cp.Section II E. We focus on correlation functions of the cavity output field, measurable in specific detector setups.For this situation, a detector (e.g., a lens which collects light over a wide angle) is modelled as an intensity measurement device on a far field surface S.

Output far field intensity
The first quantity of interest is the mean value of the output intensity where S is the detector surface and E (+) out,i is given as in Eq. ( 52) (or Eq. ( 53)).Carrying out the surface integrals by using the definition of the dissipation matrix S rad µη , and by recalling that the input state associated to, e.g., A in µ , is chosen as the vacuum state, leads to the form (cf. App.C) where µ represent the annihilation (creation) operator of the symmetrized and diagonalized QNM µ = pl, pc.

Furthermore,
Clearly the output coupling matrix L µη is connected to the radiative part of the QNM decay through the dissipation-induced radiative coupling matrix S rad µη .In contrast, in a phenomenological dissipative JC output model, one assumes typically uncoupled output radiation, where each mode is coupled out to the surrounding environment independently 48,74 .In the following, we will compare the above formulas using the derived output electric field with the output intensity in a phenomeno- for the individual output intensities in the far field of the QNMs µ = pc, pl and in the same units as I out .Note, that the factor β rad (ω a ) was added phenomenologiaclly as the radiative β factor of the hybrid system at the TLS frequency (cf.Fig. 2, (circles) and Fig. 4, (dashed line)).This is another ambiguity of the phenomenological model, since there is no clear separation between radiative and non-radiative decay processes.This is because in this case, S µη is approximated as a Kronecker-delta from the beginning and therefore, these β factors have to be added phenomenologically.In fact, taking the beta factors in a hybrid structure at certain frequencies is a highly non-trivial choice, since β rad (ω) changes drastically as a function of frequency close to ω pc , as shown in Fig. 2 (bottom).This is usually a sign for non-Markovian output characteristics.This frequency dependent beta factor is captured (at least in the bad cavity limit, where a comparison to a full Maxwell solution is possible) by the specific form of S rad µη and S nrad µη in the QNM-JC model.To underline the differences, we show the different output intensities in the steady state as function of the laser frequency ω L in Fig 9.
The peak height of I out in the full QNM-JC model at the PC frequency is roughly one order of magnitude larger than in a phenomenological treatment using the formulas in Eq. (83).Obviously, the output coupling is drastically increased in the full QNM JC-model, when the laser is tuned to the PC frequency regime and results from the increase of the beta factor in the regime close to the PC-like eigenfrequency.This has a major impact on modelling of hybrid structures for nonlinear cavity-QED experiments, since the phenomenological model (even if the system master equation is fitted appropriately) highly underestimates the output coupling.

Second-order quantum correlation functions
Next, we turn to the stationary normalized secondorder correlation functions g (2) .For a detector as a sphere S (as explained above), this correlation function reads in general where again E (+) out,i is given as in Eq. ( 52) (or Eq. ( 53)).Carrying out the surface integrals similar to I out leads to the form (cf. App.C) which is in line with the assumptions used to obtain Ĩout from Eq. ( 83).However, we note that in field of quantum plasmonics and quantum optics, the second-order correlation functions are often assumed to be for the individual modes µ = pc, pl.We emphasize, that the quantities in Eq. (87) are not observables, but can in general be useful to characterize the emitter system properties, and are often used as if they were observables 50,51 .
In Fig. 10, we show the correlation functions g(2) pl , g pc and the output correlation functions g out , g out as functions of the laser frequency ω L .For off-resonant pumping, i.e., from ω L = 1.55 eV to ω L = 1.59 eV (sector I of Fig. 10), both the system and the output correlation functions show anti-bunched character g (2) (0) < 1 of the emitted light.In particular, the output correlation functions g (2)   out , g (2)  out and the plasmon-like correlation function g (2)   pl are nearly identical to each other.This is expected, since the hybrid is dominated in the frequency regime by the plasmonic part of the two fundamental QNMs.However, close to the PC-like resonance of the hybrid, for ω L in [1.59, 1.604] eV (sector IIb in Fig. 10), the light statistics is differently predicted as bunched light g (2) (0) > 1 by the PC-system correlations, and not as anti-bunched as shown by the emitted light through the full output correlation g (2)   out .Although the phenomenological output function g (2)  out shows also anti-bunched light, it still differs slightly from g (2) out .The plasmon-like correlation function shows also anti-bunched character in this laser regime, even one order of magnitude below g (2)   out .Another interesting laser frequency regime is ω L ∈ [1.604, 1.608] eV (sector IIb in Fig. 10); here the phenomenological output function g(2) out , Eq. ( 86), predicts anti-bunched light while the full output correlation is clearly above one, showing bunched light statistics.Furthermore, both, g(2) pl and g(2) pc show anti-bunched character.Thus, only the result from QNM-JC model predicts bunched light in this small but important laser regime.Between ω L = 1.608eV and ω L = 1.8eV (sector III in Fig. 10), the light statistics is again erroneously predicted by the PC-like correlation as bunched light.In addition, the phenomenological output function g (2)  out converges much faster to the plasmon-like correlation function g (2)  pl here.In a frequency regime towards the plasmon-like frequency, the plasmon system correlations become identical to both output correlations.Therefore, the light statistics has changed considerably especially in the regime close to the PC-like mode frequency ω pc .
To summarize this subsection, the PC-like correlation function g (2)  pc is not reliable to model the second-order output correlations in the 2-mode hybrid system over the whole laser excitation frequency regime under consideration, although it seems to be the dominant part of the TLS response to the electromagnetic field at the chosen TLS frequency (cf.Fig. 3).On the other hand, at a laser frequency regime sufficiently enough away from the PCmode frequency ω pc , the plasmon-like system correlation function g (2)  pc is a very good approximation to the output correlation functions g(2) out , g out , which is expected, since the main coupling regime of the QNMs is near the PC-mode resonance.Therefore, we can conclude, that in dissipative resonator structures with at least two fundamental (overlapping) modes, the system correlation functions do not reflect the actual quantum properties of the emitted light (at least in the overlap regime), since the output is formed by a linear combination of the coupled modes, which depends on the dissipation and the radiative and non-radiative properties of the system.Furthermore, even if a β factor is added to the formulas of a phenomenological dissipative JC model (Eq.(86), it also cannot reproduce the light statistics correctly, since there is an additional non-diagonal coupling in the output quantities, induced by the off-diagonal mode coupling through the symmetrization.

IV. CONCLUSIONS
We have presented a detailed extension of the QNM quantum model from Ref. 40 by including external pumping and by deriving an generalized input-output theory for multiple QNMs.Explicit expressions for correlation functions outside of the nanostructure were provided.Furthermore, pronounced differences in the correlation functions inside the system and at the outside detector were found.We analysed the cavity-QED behaviour of a metal-dielectric hybrid resonator coupled to a TLS in the intermediate coupling regime.We compared our full microscopic QNM model including mode coupling with a phenomenological dissipative JC model.Significant qualitative and quantitative differences, induced by the intermode coupling, were found between both models.
In the nonlinear pumping regime, we also studied the response of the hybrid cavity system to an external laser (in a excitation regime, that allows to study the twoexcitation manifold of the JC ladder), and we found that the phenomenological dissipative JC model is mainly in the single-photon regime, while the driven QNM-JC model has multiphoton character.Thus the quantum dissipation effects, that appear in the quantized QNM model, induce higher-order correlations.This shows that the inter-mode coupling, coming from the microscopic QNM model, is crucial to include in the quantum master equations.
coupling part for TLS frequencies near the PC-mode frequency, Ω ∼ g pc ∼ Γ pc indicates the intermediate coupling regime.However, we emphasize, that the hybrid system described by the phenomenological dissipative JC model has a larger TLS-PC coupling compared to the QNM-JC model, i.e., gpc ≈ 5g pc .Therefore the excitation regimes are slightly different.Thus, another possibility to compare both models would be to choose a Rabi frequency, that scales with the different TLS-PC coupling constants, i.e, choosing a different Rabi frequency for the QNM-JC model and the phenomenological dissipative JC model.To discuss the different ways of comparison, we show the peaks of the eigenstate probabilities P pc , P pc−pc (discussed in subsection III C) as function of Ω L in Fig. 11.In Fig. 11 (top), we show results for the treatment of the Rabi frequency, which is used in the main part, i.e. we choose the same external laser strength for both models.In this case, the same hybrid response completely different to an external laser in both models.In contrast, in Fig. 11 (bottom), a treatment, where the Rabi frequency scales with the different TLS-PC coupling constants is shown: Here, the probabilities are quiet similar, but a significant quantitative difference

FIG. 1 .
FIG. 1.(a) Top view on the metal-dielectric hybrid structure, consisting of a metallic dimer on top of a photonic crystal beam.(b) Side view on the metallic dimer, supporting a quantum dipole (z-polarized) in the middle of the gap.(c) Schematic of the quantum subsystems: The photonic crystallike mode (violet) is coupled to the plasmonic-like cavity mode (yellow) and the TLS (red) is coupled to both symmetrized QNMs.d) Excitation and input-output scheme.The TLS is driven by a cw-pump field with strength ΩL and decays with a rate γSE, while the QNM subsystem is characterized by the decay rates Γ pl and Γpc, naturally entering the model through the input fields A in pl and A in pc , respectively (see text).

FIG. 2 .
FIG.2.Classical Purcell factors (first and second row) and classical radiative β-factors (third and fourth row) from an embedded dipole emitter in the metal-dielectric hybrid structure from Fig.1with gap size of 2 nm and 5 mn, obtained from full Maxwell simulations over a frequency regime covering both hybrid resonances (left) with a zoom-in close to the high-Q resonance (right).Note that for the 5 nm (2nm) gap case, the numerical calculations were done for 81 (93) non-equidistant frequency points in the interval ω ∈ [1.55, 1.95] eV, and that the corresponding F num p

FIG. 3 .
FIG.3.Purcell factor of a z-polarized dipole emitter in the metal-dielectric hybrid structure from Fig.1, shown over a broad frequency range from full dipole simulations (dots), the two-mode QNM-JC model (solid magenta line) and the phenomenological dissipative JC model (dashed green line).Two peaks appear around the hybridized QNM frequencies ωpl = 1.6999 − 0.0479i (eV) and ωpc = 1.6052 − 0.0007i (eV), originating from the metallic ellipsoidal dimer and the photoniccrystal beam, respectively.The highly non-Lorentz interference effect (dip in Fp) is located near one of the eigenfrequencies, (pc) em , of the full electromagnetic Hamiltonian Hem of the QNM-JC model and is fully recaptured by the inter-mode coupling terms.

FIG. 4 .
FIG.4.Quantum modal β-factors of the metal-dielectric hybrid structure from Fig.1, shown over a broad frequency range for the two-mode quantum QNM model (solid magenta line, Eq. (72)) and the phenomenological dissipative JC model (green dashed line with β ≈ 0.14).Note, that for the JC-model, a phenomenological introduced constant is shown, which mimics the behaviour of the radiative output at TLS frequencies around ωpc with respect to the classical β rad class solution from Fig.2.The inset shows a zoom-in around the PC-like mode frequency ωpc = 1.6052 (eV).In the QNM-JC model, a peak close to the PC-like mode with frequency ωpc appears, which is located at the position of the highly non-Lorentz interference effect (dip in Fp) from Fig.3.

( 2 )
em aligned with the eigenvalue of the electromagnetic Hamiltonian Hem close to the PC frequency.The solid (QNM-JC model) and dashed (phenomenological dissipative JC model) bars reflect the contribution of the (bare) eigenstates |k pl , kpc, l of the uncoupled photon-exciton system to the eigenstates |ϕN,i N of the actual coupled photonexciton system.The eigenstates of the coupled system are similar to the respective bare states, since the coupling constants indicate the intermediate coupling regime (cf.Tab.I).(H JC sys ), we recall the bare state basis, in which the density operator is expanded, i.e. |k pl , k pc , l for the QNM-JC model and | kpl , kpc , l for the phenomenological dissipative JC model.Here, |k pl (| kpl ) represents the number state of the plasmon-like mode, |k pc (| kpc ) is the number state of the PC-like mode and l = g, e denotes the state of the TLS in the QNM-JC (phenomenological dissipative JC) model.These bare states are eigenstates of the photon number operators npl = A † pl A pl ( npl = α † pl αpl ), npc = A † pc A pc ( npc = α † pc αpc ) and the TLS number operator na = σ + σ − (corresponding to the occupation of the upper level |e ) with non-degenerate eigenvalues, respectively.However, they are not eigenstates of the Hamiltonians H QNM sys (H JC sys ), since the subsystems are coupled.The total number operator, e.g. for the QNM-JC model, N = npl + npc + na constitutes of the 2N + 1-dimensional eigenspace

1 , 1 =FIG. 6 .
FIG.6.(a) First (3-fold) and second (5-fold) rung of the QNM-JC energy ladder with ωa aligned to the eigenfrequency of the electromagnetic Hamiltonian Hem close to the PC mode frequency: The solid lines show the real part of the complex eigenenergies ( N,i N ), while the grey area covers the imaginary part and is bounded by Re( N,i N ) ± 0.2Im( N,i N ).(b) Eigenenergies 1,2 and 2,2: The solid (lightblue dashed) line shows the real part of the eigenenergies in the QNM-JC model (phenomenological dissipative JC model), while the grey (lightblue dashed) area reflects the imaginary part and is bounded by Re( N,i N ) ± Im( N,i N ).Direct 2-photon transitions via an external laser with frequency ωL(ω L ) are sketched.Note, that in both subfigures, there is a sudden jump on the energy axis from the first to second rung and that the scaling is different in the upper and lower half.
) yield a subset of the eigenenergies of the full Liouvillian L QNM = H QNM sys + L QNM diss (cf.Eq. (

FIG. 7 .
FIG.7.Occupation numbers na, npc, n pl and probabilities Pa, Ppc, P pl of the TLS (corresponding to upper level |e ), plasmon and PC-like mode in the steady state, obtained with the QNM-JC model (solid) and phenomenological dissipative JC model (dashed) as functions of laser frequency ωL around ωpc.The TLS frequency ωa = (2ω L ) close to the peak of the 2-photon probability P PC−PC of the QNM-JC model in the eigenenergy diagram from Fig. 6 (b), showing the eigenstates |ϕ pc and |ϕ PC−PC :

FIG. 9 .
FIG. 9. Output intensity function I out measured on a sphere in the far field as function of laser frequency ωL for the QNM-JC model (solid) and the phenomenological dissipative JC model (dashed) in the steady state with the TLS frequency ωa = (pc) em aligned to the eigenenergy of the full electromagentic Hamiltonian Hem.The Rabi frequency of the external pump is ΩL ∼ |gpc|, where the TLS-PC coupling is in the intermediate couling regime, i.e., |gpc| ∼ Γpc (cf.Tab.I).The lower plot shows a zoom-in close to the PC-mode frequency ωpc.The quantities are normalized to max(I out ) from the QNM-JC model.