Single Quantum Emitter Dicke Enhancement

Coupling $N$ identical emitters to the same field mode is well-established method to enhance light matter interaction. However, the resulting $\sqrt{N}$ boost of the coupling strength comes at the cost of a"linearized"(effectively semi-classical) dynamics. Here, we instead demonstrate a new approach for enhancing the coupling constant of a \textit{single} quantum emitter, while retaining the nonlinear character of the light-matter interaction. We consider a single quantum emitter with $N$ nearly degenerate transitions that are collectively coupled to the same field mode. We show that in such conditions an effective Jaynes-Cummings model emerges, with a boosted coupling constant of order $\sqrt{N}$. The validity and consequences of our general conclusions are analytically demonstrated for the instructive case $N=2$. We further observe that our system can closely match the spectral line shapes and photon autocorrelation functions typical of Jaynes-Cummings physics, hence proving that quantum optical nonlinearities are retained. Our findings match up very well with recent broadband plasmonic nanoresonator strong-coupling experiments and will therefore facilitate the control and detection of single-photon nonlinearities at ambient conditions.


I. INTRODUCTION
Strongly coupled light-matter systems, described by the Jaynes-Cummings (JC) model [1], are of fundamental interest in testing the quantised nature of coupled light and matter degrees of freedom in the context of cavity QED [2]. Importantly, JC physics requires that a single two-level emitter interacts with a single-mode electromagnetic field. Such systems can be relevant for quantum information processing, for example in interfacing flying and static qubits [3] or for quantum-state engineering [2]. Yet strong coupling of light and a single emitter is not easy to achieve in practice, due to the large mismatch in the spatial extension of free-space photons and typical emitters (i.e. artificial or real atoms) that feature electronic transitions in the optical domain [4,5].
A possible route to overcome this difficulty is to engineer an "effective" emitter with large dipole moment. This is most commonly implemented by combining N 1 identical emitters coupled to the same field, resulting in the well-known "Dicke enhancement" of light-matter interaction: the effective coupling strength is boosted by a factor of order √ N, compared to the single-emitter case [6][7][8][9]. However, the relevant theory in such systems is the many-emitter limit of the Tavis-Cummings (TC) model [10], whose low-energy dynamics is effectively linear, i.e., indistinguishable from that of coupled harmonic oscillators [11].
The key achievement of our work is to bypass this seemingly inevitable trade-off between coupling strength and (quantum) nonlinearity: we open an innovative route to a 'Dicke-like' enhancement of the coupling constant which preserves the precious quantum optical nonlinearities of the JC model. * tommaso.tufarelli@nottingham.ac.uk † hecht@physik.uni-wuerzburg.de In this paper we indeed show how the exploitation of a multilevel emitter, with N nearly degenerate excited states coupled to the same field mode, increases the effective lightmatter coupling constant by a factor of oredr √ N. While this is the same scaling found in the Dicke [12] and TC models [10] with N emitters, the crucial difference is that in our case the quantum non-linearity of the interaction is preserved: we demonstrate this based on emission spectra and second-order photon autocorrelation functions. More precisely, our system can closely approximate the behaviour of a JC model, in principle, for arbitrary values of N. These findings, in their simplest form, are schematised in Fig 1. Additionally, we find that our predictions are robust against typical sources of decoherence (dephasing) that are common at elevated temperatures. Thus, while our results are fully general and independent of the physical implementation of the model, we anticipate that our approach can be particularly useful in the context of room temperature strong coupling. We thus find it useful to provide some context on the current state-of-the-art of this research field, and to highlight how our research can help it move forward.
Recently, room-temperature strong coupling with many emitters [13][14][15][16] and even single quantum emitters [17][18][19][20][21] have been experimentally demonstrated. These experiments rely on the remarkable coupling strengths achievable in plasmonic nanoresonators, thanks to the deeply sub-wavelength mode volumes of the latter. In order to relax the mode volume requirements (which are currently pushing the limits of nanofabrication techniques) and obtain more robust experimental realizations, it would be of great interest to find emitter systems with increased dipole moments. This would result in larger coupling constants while still allowing to harness the single-emitter quantum nonlinearity of the JC model (JCM) [22,23]. So far, only the "first rung" of the JC ladder, i.e. the single-excitation regime of the model, has been experimentally detected at ambient conditions. Efforts are still ongoing Standard Jaynes-Cummings model, in which a single two-level emitter couples resonantly to a single mode field, with strength g. The resulting anharmonic spectrum is a cornerstone of modern quantum optics, but challenging to probe experimentally: the required regime g κ is difficult to achieve, where κ is the cavity decay rate. Center: Multilevel model, the focus of this paper. A quasi-degenerate multilevel emitter, with N closely spaced excited states, interacts with a single mode field. When all the transitions are approximately resonant with the field, an effective coupling constant g eff , enhanced by a factor ∝ √ N, is established. In principle, the resulting low energy spectrum can closely approximate the JC spectrum for any value of N. Right: Tavis-Cummings and Dicke Models, in which N 1 nearly-identical two level systems interact with the same field. A boosted effective coupling constant g eff ∝ √ N is established also in this case, but the spectrum becomes approximately harmonic as N is increased.
to detect higher rungs of the ladder to shed light on (and harness) the quantum non-linearity of the model.
In order to realise the JCM experimentally, great efforts have been made over the years in the fabrication of artificial atoms -such as semiconductor quantum dots [24] -that are close to ideal two-level emitters. Indeed, so far the multi-level nature of such emitters has mainly been treated as a modelling nuisance, i.e. as a source of errors with respect to an idealised two-level scenario. In this context, Madsen et al. [25] and our previous paper [19], have recently mentioned the possible effects of a more complicated energy level structure of the emitter. To the best of our knowledge, however, our previous work [19] was the first to explicitly suggest that an emitter with multiple (> 2) energy levels can even bring about some advantages: in the right conditions it may be used to simulate a JC system with enhanced coupling constant. With reports of quantum dots featuring as much as 64 quasi-degenerate excited states [26,27], a general quantum-mechanical analysis of such a system (multilevel model for brevity) is in order.
The findings presented here will thus be of importance for the selection and the design of suitable emitters for nanoscale strong-coupling experiments at ambient conditions and for the interpretation of such experiments. Importantly, our approach lives independently of the contentious issues regarding the calculation of mode volumes and coupling constants in quantum (and classical) nanophotonics, see for example [28][29][30]. These approaches will of course determine the numerical values of coupling constants in different experimental situations, but they do not affect the nature of our quantum optical model and the general conclusions we can draw from it.
The paper is organised as follows. In Section II we introduce a Hamiltonian that models an emitter with N excited levels coupled to the same ground state via a single mode field. We will show that our system can be understood as a JC model which is weakly interacting with N − 1 "dark states" of the emitter. In Section III we generalise our model to an open system via a master equation. This takes into account incoherent processes such as photon loss, dephasing and emitter pumping, all of which may be required for describing realistic experimental scenarios. In Section IV we support our general arguments with detailed analytical results for the special case N = 2: we find that our system's emission spectrum features approximately the same resonances as a JC model, plus additional features (due to the dark states) that we fully character-ize. Section V presents numerical studies for N = 3, comparing common experimental observables in our model with the same quantities calculated in reference JC and TC/Dicke models. In Section VI we draw our conclusions. For completeness, Appendix A contains some useful reminders and formulas pertaining JC, TC and coupled-oscillator models, since those are referenced throughout the paper.

II. MODEL
Our multilevel model features a quantum emitter with electronic ground state |G and a collection of excited states |e k , also called sublevels, where k = 1, 2, . . . , N. Each internal transition of the emitter, |G ↔ |e k , is coupled to a singlemode field (cavity for brevity) with strength g k ∈ R. The typical physical system we have in mind is a multilevel quantum dot embedded in a plasmonic nanoresonator, where the cavity resonance is broadband, in the sense that it overlaps with all the emitter transitions, and where strong field gradients may allow for light-matter couplings beyond the dipole approximation (hence the possibility of coupling a large number of levels). Yet, it is important to note that our quantum optical formalism applies to any type of emitter or resonator that satisfy our modelling assumptions. As usual the cavity mode is described via the annihilation operatorâ, with [â,â † ] = 1.
Taking = 1 for convenience, the cavity resonance has average frequency (hence energy) ω 0 , the energy of |G can be set to zero without loss of generality, while each state |e k has energy ω 0 + ∆ k . Here ∆ k is the detuning between the k-th emitter transition (|e k → |G ) and the cavity. For convenience we assume ∆ 1 ≤ ∆ 2 ≤ ... ≤ ∆ N , we indicate the full range of energy detunings as ε ≡ ∆ N − ∆ 1 and the average detuning as ∆ = 1 N k ∆ k . we can write the total Hamiltonian as where H 0 , V are, respectively, the free and interaction Hamiltonians, and we have neglected counter-rotating terms in the light-matter interaction. The physics of Hamiltonian (1) can in general be quite complicated, despite H being numberconserving -i.e. H commutes with the number operator N tot =â †â + k |e k e k |. However, here we are interested in the special case of quasi-degenerate emitters that are nearresonant with the cavity, i.e. we assume the parameter regime |∆ k | g k . This approximation holds e.g. for colloidal quantum dots coupled to plasmonic nanoresonators. Here, ∆ k are on the order of a few meV, while the coupling strengths, g k , are on the order of 100meV [19]. In this scenario a number of qualitative observations can be made. In the bare emitter basis the two levels |e 1,2 are not directly coupled to each other, while ε plays the role of a relative detuning between them. Both |G ↔ |e 1,2 transitions couple to the cavity field, and we are assuming for simplicity that the associated coupling constants are the same: In the radiation basis, the bright and dark states |B , |D have the same frequency, but are coupled to each other with strength ε/2. The transition |G ↔ |B is coupled to the field with enhanced strength 2g, while the dark state |D does not directly interact with the field.

A. Bare emitter basis versus radiation basis
To aid understanding, the main ideas behind this subsection are illustrated in Fig. 2 for the simplest case, N = 2. We begin our study of the general N-sublevel model by considering the case of exact degeneracy of the excited states, ∆ k = ∆ ∀k, or equivalently ε = 0 in terms of the full range of energy detunings. The resulting mathematics follow in a straightforward manner from the well-known JCM. This simple starting point will set the stage to investigate small deviations from an idealized scenario, i.e. cases with small but nonzero ε. For ε = 0, we find that any linear combination of the states |e k is itself a valid emitter excited state, or more precisely, an eigenstate of H 0 with eigenvalue ω 0 + ∆. Among all the states that can be constructed in this way, of particular significance is the coherent superposition where the denominator ensures normalization. We shall refer to |B as the bright state, or superradiant state (mimicking the terminology of the Dicke model). State |B can be interpreted as the only emitter excited state that couples directly to the cavity field. Indeed it is easy to check that the light-matter interaction term in Eq. (1) may be rewritten in the simple form where we have defined an effective coupling constant Hence, through an appropriate change of basis for the emitter, the N degenerate transitions |G ↔ |e k can be replaced by a single transition |G ↔ |B , characterized by an enhanced coupling to the cavity field. Eq. (6) indeed shows that the effective coupling scales like g eff ∼ √ N with the number of sublevels. This boosting of the effective light-matter coupling is analogous to what happens in the Dicke [12] and TC [10] models, where N degenerate emitters are coupled to the same field. A crucial difference from the TC model is that Eq. (5) retains the structure of a JC interaction Hamiltonian, so that in the model studied here there is in principle no loss of anharmonicity as the number of sublevels N is increased.
To fully specify the mentioned change of basis, we can define N −1 further linear combinations of the excited states |e k , say {|D 1 , |D 2 , . . . , |D N−1 }, which together with {|G , |B } form a complete orthonormal set for the emitter internal levels. When ε = 0 each |D k is also an eigenstate of H 0 with eigenvalue ω 0 + ∆, due to the N-fold degeneracy. By construction, these additional basis states do not directly interact with the cavity field, i.e.
and may be called subradiant states or dark states. In practice, explicit expressions for the dark states may be found by standard orthogonalization procedures (details not shown). The above discussion suggests that the set {|G , |B , |D 1 , |D 2 , . . . , |D N−1 }, which we will call radiation basis, may be a more meaningful conceptual tool to study light-matter interaction in our system, as compared to the bare emitter basis {|G , |e 1 , |e 2 , . . . , |e N }.
For the purpose of modelling the optical properties of the system, it is now tempting to simply neglect all the dark states and retain only the two emitter levels |B , |G . This is justified if there are no additional processes, coherent or incoherent, that are able to populate the dark states and/or couple them to the two levels |B , |G . When all these conditions are satisfied, the system is exactly described by a JCM with coupling constant g eff .
We are now in a position to study small deviations from the ideal case. As anticipated, the most obvious deviation from pure JC physics occurs when the excited states |e k are not exactly resonant with each other (i.e. ε 0). In the radiation basis, the (small) relative detunings between the levels |e k will translate into (weak) couplings between the bright and dark states, and between the dark states themselves. In other words, |B and |D k cease to be eigenstates of H 0 when ε 0. Note however that Eqs. (5) and (7) remain valid in general, so that the dark states may interact with the field only indirectly, i.e. through the mediation of state |B . These additional interactions will bring about deviations from the well understood JC physics described above. While a general calculation of these effects would be cumbersome and beyond the scopes of this paper, the case N = 2 sublevels is simple enough to be treated analytically and already contains all the essential ingredients in our problem -see Section IV below and Fig. 2.
In some cases, the exact equivalence to a JC system may be lost even in the case of perfect degeneracy between the excited FIG. 3. Schematic of the investigated setup. The system may be excited through an incoherent pump of intensity Λ applied to the emitter, or coherently driving of the (plasmonic) cavity with field amplitude Ω and frequency ω L . The bright state of the emitter is coupled to the cavity with the effective strength g eff , and the emitter sublevels are spread in frequency over a range ε. κ is the loss rate of the cavity and γ the decay rate of all the emitter sublevels, which are also subject to dephasing at rate γ d (not shown). All the observables discussed in this work can be measured by collecting the cavity output field.
sublevels. This can happen if incoherent processes couple the pair {|G , |B } to the dark states. We shall see in Section IV, for example, that emitter dephasing provides a common mechanism for interconverting bright and dark states. Yet, our analytical and numerical results suggest that the multilevel model can provide an excellent approximation to JC physics even when these imperfections are taken into account.

III. MASTER EQUATION
In order to predict relevant observables for future experiments, we now introduce an open system generalisation of our model, via the well established Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) master equation [31]. As anticipated above, the main experimental scenario we have in mind is a quantum dot interacting with a plasmonic nanoresonator, yet the ideas presented here may be applicable to a wider variety of setups. We will include cavity photon loss as the primary channel for energy dissipation, while the decay rate of the emitter into free-space modes will be assumed negligible (however, in simulations it will remain finite for reasons of numerical stability). Such an assumption is reasonable due to the large β-factors achievable in plasmonic resonators. On the other hand, we shall assume that the emitter excited states can be broadened well beyond their radiative linewidth by additional dephasing processes -typically thought of as a drawback of room temperature operation. We consider two ways in which energy can be supplied to the coupled system: incoherent pumping of the emitter, or a coherent driving field applied to the cavity. While the first is a common method to feed excitations into quantum dots, the second is arguably the standard avenue to test nonlinearities in a quantum optical sys-tem. For a visual summary, the main ingredients of our open system model are sketched in Fig. 3. In formulas, the following master equation will be used to model the emitter-cavity system, where H is the Hamiltonian (1), and is an interaction term describing coherent driving of the cavity by a classical laser at frequency ω L (Ω quantifies the driving field strength). The L's appearing on the right hand side of Eq. (8) are the jump operators, used to model incoherent processes. They are collected in the set J, which features the following elements: • Cavity decay (photon loss) at rate κ: L = √ κâ • Radiative decay of the emitter γ: L = √ γ |G e k |, • Emitter pumping of (total) rate Λ: L = Λ N |e k G|, where k = 1, 2, . . . , N. We assumed that the decay and dephasing rates are the same for all sublevels, and that the pump power Λ is split symmetrically over the N sublevels. We also consider equal couplings N for all k, while the detunings ∆ k are equally spaced in the interval [∆−ε/2, ∆+ε/2]. All the assumed symmetries between the sublevels are of course an idealisation. Many generalisations of this basic model can be included in principle. However, we do not expect that these can bring about significant qualitative departures from the physics described here. Hence, in this work we confine our analysis to the above setting for simplicity and computational efficiency.

IV. ANALYTICAL RESULTS FOR 2 SUB-LEVELS
We here present the analytical diagonalization of Hamiltonian (1) for the case N = 2 with equal couplings g 1 = g 2 = g, as well as the associated implications for the system's open dynamics. This instance of the model, illustrated in Fig. 2, is particularly transparent and fully captures the essence of our study. We verified that analytical solutions are also available for the cases N = 3, 4, using similar methods, but these are significantly more complex without adding fundamentally new insights. For N = 2 the detunings are given by ∆ 1 = ∆ − ε/2 and ∆ 2 = ∆ + ε/2, respectively. Furthermore there are only one bright and one dark state given by |B = |e 1 +|e 2 √ 2 , |D = |e 1 −|e 2 √ 2 . In the radiation basis representation, the Hamiltonian may be recast as where in this case g eff = √ 2g. As for the general case treated in the previous section, a JC interaction describes the coupling between the two levels {|G , |B } and the cavity field. However we can also see that the relative detuning of the sublevels, ε, translates into a small coupling between |B and |D , so that even the dark state |D can indirectly influence the optical behaviour of the system. The ground state of H is easily spotted as |G, 0 , with eingenvalue E = 0. Exploting excitation number conservation, i.e. [H,N tot ] = 0, wherê N tot =â †â + |B B| + |D D|, one can show that the remaining eigenvalues and eigenvectors of H can be found by diagonalising the 3 × 3 matrices (n = 1, 2, . . . , ∞) which describe the Hamiltonian restricted to the subspace The behaviour of the eigenvalues of H 1 is shown in Fig 4, as a function of the detuning ∆. Even at the non-negligible ratio ε/g eff = 0.25 used in the figure, the commonalities with a JCM spectrum are evident. Note also that the eigenvalues of H n for n ≥ 2 will follow the same qualitative structure, as a generic H n may be obtained from H 1 by a simple rescaling g eff → g eff √ n, ω 0 → nω 0 . While a general analytical diagonalization of Eq. (13) is cumbersome, for the resonant case ∆ = 0 one can obtain simple analytical expressions. For this case, the eigenvalues of H n read E n,± = nω 0 ± g eff n + ε 2 The corresponding normalised eigenvectors |E n,± , |E n,D are The above expressions suggest that in the regime ε g eff it is natural to identify {|E n,± , E n,± } as a weakly perturbed JC eigensystem (the latter is indeed recovered exactly for ε → 0). For brevity, we refer to those as the 'JC model-type' eigenvalues and eigenvectors. The eigenstates |E n,D , which instead feature the same energy levels of an empty cavity, may be understood as dressed quasi-dark states of the emitter, accompanied by approximately n − 1 cavity photons. They indeed result from a hybridization of states |G, n and |D, n − 1 , with FIG. 4. Eigenvalues of H 1 as a function of the average detuning ∆, for ε = 0.25g eff (blue continuous lines). The anticrossing behaviour of eigenvalues E 1,± around ∆ = 0 approximates well (the maximum relative error in the plotted range being ∼ 1.3%) the reference JCM with coupling g eff (red circles). The energy E 1,D of the quasi-dark state is instead approximately linear in the detuning and sits in between the JCM-like energies. An elementary perturbation theory calculation indeed yields E n,D ∆ 1 − ε 2 /4g 2 eff n , which agrees well with the exact value in the considered regime (relative error below 0.2%). Qualitatively similar results are obtained for higher excitation numbers n ≥ 2, and in fact it is easy to see that the relative weight of ε and ∆ decreases with n.
the second state being dominant in the superposition. It is interesting to notice that the two states involved in such superposition are only coupled indirectly in the Hamiltonian, through the mediation of state |B, n − 1 . Exploiting the analytical results obtained here, it is possible to calculate the Hamiltonian dynamics of a generic initial state via standard methods. In fact, while our work focuses on the regime ε g eff , we note that the above diagonalization procedure is valid for a general choice of the model parameters. The closed expressions obtained so far are particularly helpful in analysing the open system behaviour of our multilevel model, as well as its relationship to the JC model. Since in this work we focus on the regime of low pumping (or low coherent driving) and negligible emitter decay, we can gain valuable intuition about the properties of the master equation (8) by setting γ = Λ = Ω = 0. First, we observe that the introduction of a photon loss rate κ results in a finite lifetime for the excited states of H, such that eigenstates withN tot = n will tend to decay towards states withN tot = n − 1, until the ground state |G, 0 is finally reached. For a given eigenstate |E j , an estimate of the decay rate induced by photon loss can be obtained via the perturbation theory formula which at second order in ε yields (n ≥ 1) We may immediately notice that the lowest quasi-dark state, |E 1D , is long lived, while all the remaining excited states have approximately the same lifetime of either the JCM excited states or the Fock states of an empty and lossy cavity. Next, emitter dephasing will induce further broadening of the Hamiltonian energy levels, which can be estimated as Some straightforward algebra yields where again we have expanded our results at second order in ε. In this case we can see that the quasi-dark states are the most affected, since they feature a larger overlap with the bare excited states of the emitter. Note also that the δE j 's are not associated with a decay towards states of lower excitation number (more below). The exact broadening of the energy levels, induced by the combined effect of photon loss and emitter dephasing, can be obtained as −2Im(Ẽ j ), wherẽ E j are the eigenvalues of the non-Hermitian operator The resulting expressions are unwieldy, however we checked via a numerical optimization that −2Im(Ẽ j ) Γ j + δE j , with a relative error below 3.1% in the range 0 ≤ ε ≤ 0.25g eff , 0 ≤ κ ≤ 2g eff , 0 ≤ γ d ≤ 2g eff , which is more than enough for our scopes.
In order to gain qualitative understanding of photon emission processes in our system, we next look at the (squared) matrix elements of the annihilation operatorâ, between an initial energy eigenstate |E i and a final one |E f . The quantum jump approach to GKSL master equations [31] provides a useful interpretation for these quantities: A i→ f is proportional to the probability that the transition |E i → |E f occurs upon loss of a cavity photon [in the master equation (8), such transitions are triggered by the term κâρâ † ]. When the transition occurs, a photon is emitted into extra-cavity modes with (approximately) frequency ω f i = E f − E i and bandwidth The A i→ f terms involving only JCMtype eigenstates do not present any particular surprises: they approximate the well-known JCM results for ε g eff , and converge to them for ε → 0. We thus omit such cases for brevity and focus only on the matrix elements involving quasi dark states, which yield the following expressions: A nD→n−1,D = n − 4g 2 eff n 4g 2 eff n + ε 2 , From Eqs. (27), (28) and (29) we notice that, once the system is in one of the quasi-dark states |E nD , the loss of a cavity photon can only produce the quasi dark state |E n−1,D for n ≥ 2, or |G, 0 for n = 1. This corresponds to the emission of a photon at the cavity frequency ω 0 , with the associated matrix element of the order ∼ n − 1. Assuming for the moment γ d = 0 (no dephasing), this implies that once the system has evolved into a quasi-dark state, it essentially reproduces the open dynamics of an empty cavity until state |1, D is reached. From there the system will decay slowly towards the lowest state |G, 0 , resulting again in the emission of a photon at the cavity frequency ω 0 , but with a much narrower bandwidth Γ 1D κε 2 /4g 2 eff . On the other hand, from Eq. (30) we can also see a small probability that a JCM-like state decays into a quasi-dark state by loss of a cavity photon. This process results in the emission of photons at average frequency ω 0 ±g eff √ n: these exactly overlap with the JCM spectral resonances for n = 1, while they provide novel emission peaks for n ≥ 2. Putting back finite dephasing γ d 0 into the picture, we see that, in addition to introducing further broadenings of the levels, dephasing events are capable of interconverting states |B and |D . More in detail, we find | f |L|i | 2 = γ d /2 for L = √ γ d |e k e k | with k = 1, 2 and i, f ∈ {D, B}. Referring again to the quantum jump approach [31] this will swap the states |B and |D with probability 1/2, if a dephasing event occurs. In summary, we have provided a detailed characterization of the processes occurring in our open system for the emblematic case N = 2, ∆ = Λ = γ = Ω = 0. While the multilevel model can approximate all features of a JCM for ε g eff , it displays additional phenomena due to the indirect interaction of dark states and cavity field. In the next section we explore these issues quantitatively via numerical simulations.

V. CLIMBING THE JC LADDER
In the following numerical studies, we will investigate how closely the physics of our system -as relevant for typical quantum optics experiments -resembles a JCM with enhanced coupling constant g eff ∝ √ N. For the sake of comparison, we shall also display the behaviour of a TC model with N atoms and the same effective coupling constant g eff -see Appendix A for further details. Simulations are done for N = 3: this is one step up in complexity compared to the previous section, but simple enough that the resulting calculations can be handled by a standard computer.
Our analysis focuses on two observables that are ubiquitous in both theoretical and experimental investigations of the JC model (proper definitions will be given below): • The steady state spectrum S a (ω) of the cavity output field, when the system is excited via incoherent pumping of the emitter; • The photon autocorrelation function g (2) (0), when the emitter is excited via a coherent driving field with variable detuning.
The first observable provides valuable information about the energy level structure of a quantum system, in particular highlighting its optically allowed transitions. Instead, the g (2) (0) may be seen as a witness of optical non-linearity: we will use it in our model to show that the √ N enhancement of the coupling constant, g eff ∝ √ N, does not come at the cost of "linearising" the light-matter interaction (something that instead happens in the TC scenario as N is increased, as we shall see below). To support this last point, we will indeed compare our autocorrelation functions with the g (2) (0) obtained for a reference "coupled oscillator model" -i.e. the quintessential linear system in our context (see Appendix A for details). Note that we are focusing only on observables that can be measured via the cavity output field. The latter has indeed the same definition in all models considered here, facilitating direct comparisons between them. On the other hand, emitter observables do not have an unambiguous correspondence in the different models, as they are defined on Hilbert spaces of different dimensions. In detail, according to standard input-output theory [32], the portion of the cavity output field that reaches our detector, saŷ E out (scalar for simplicity) may be expressed as a function of the cavity annihilation operatorâ as follows: where E c is a constant depending on the details of the cavityplus-detector setup, while the operatorÊ vac captures all "vacuum noise" contributions from the detector's environment. Thanks to the assumption of a vacuum environment, which is justified for optical frequencies at room temperature,Ê vac does not contribute to any normally ordered observable constructed fromÊ out -and these are precisely the types of observables employed here.
In the interest of brevity, we restrict ourselves to an experimentally desirable situation in which ∆ = 0 and the hierarchy ε κ g eff holds. In other words we consider a strongly coupled system where the average transition frequency of the emitter is at resonance with the cavity, while the range of sublevel detunings is within the bare cavity bandwidth. This scenario is indeed appropriate for the case of quantum dots embedded in plasmonic nanoresonators, as anticipated in Section II. We shall also take the driving strengths (Ω for the coherent and Λ for the incoherent case) as the smallest parameters in the problem. This ensures that the driving will not induce significant modifications to the intrinsic resonances of the coupled system, nor to the conditions for achieving strong coupling [33,34]. In this setting we focus on two representative parameter regimes: (i) g eff = κ, or 'single-excitation regime', (ii) g eff = 20κ, or 'bi-excitation regime'.
Loosely speaking, case (i) gives access to the 'first rung' of the energy ladder, thus the resulting Physics should be similar to that of a pair of coupled oscillators in all models. Case (ii) instead allows us to explore the 'second rung' of the ladder, where the anharmonic nature of the investigated models becomes more prominent. In both situations we shall analyse the role of emitter dephasing by considering three scenarios: no dephasing, low dephasing (i.e. one order of magnitude smaller than the cavity decay rate κ) and high dephasing (i.e. dephasing and cavity decay rates of the same order). This is particularly relevant for quantum-dot implementations of our models, especially at room temperature, in which case emitter dephasing rates can be comparable to the cavity decay rate. Finally it is important note that we shall adopt a logarithmic scale in our plots, since the examined quantities display features of very different magnitudes. This implies that the weaker features displayed in our plots (e.g. the second-rung spectral resonances) may be challenging to detect in experiments.

A. Steady state spectra
As anticipated, our first set of numerical examples investigates the optical emission properties of our multilevel model under incoherent pumping of the emitter, as relevant in many quantum dot experiments also at ambient conditions. Specifically, we set Ω = 0 (no coherent driving) and Λ 0. In this setting we are interested in the steady state spectrum of the cavity output field, i.e., where Eq. (31) has been used, the two-point correlation function â † (τ)â(0) ∞ is calculated via the quantum regression theorem [31], and the expectation value is calculated on ρ ∞ , the steady state of the master equation. The latter is found by solving the linear system Lρ ∞ = 0, Tr[ρ ∞ ] = 1, where the superoperator L is implicitly defined by rewriting the master equation (8) asρ = Lρ. We perform such calculations numerically with Python, by truncating the Hilbert space of the cavity to a sufficiently high dimension to obtain numerical convergence. In addition to the cavity spectrum, we shall also be interested in the quantity i.e. the steady-state probability that the emitter has settled in a dark state (or a mixture thereof). The resulting spectra, in logarithmic scale, are plotted in Fig. 5 together with those of the reference TC and JC models (again, see Appendix A for details). For a fairer visual comparison between the three models we have also rescaled the spectrum of JC and TC models by a factor 1 − p dark . We indeed recall from Section IV A that the multilevel model cannot display JC physics once it gets "stuck" into a dark state.
In the single excitation regime, we find that all three models display the expected resonances at ω ω 0 ± g eff , but the multilevel model presents additional sharp features at ω ω 0 . In the case of no dephasing, these resonances account for a significant portion of the total emitted energy, and correspondingly p dark is close to one. We argue that the ω ω 0 features can be ascribed to the slow decay of the long-lived quasi-dark states of the single-excitation sector: referring to our discussion and notation in Section IV A, this involves transitions of the form |E 1D → |G, 0 (where |E 1D can be any of the two quasi-dark states occurring for N = 3). As dephasing is increased, going left to right in Fig. 5, we find that the three models show increasing agreement, and at the same time p dark is significantly reduced. Referring again to our analytical discussion in Section IV A, we recall that dephasing is indeed able to swap bright and dark states. This opens a new decay channel for the quasi-dark states |E 1D : as dephasing is increased, we may reach a parameter regime where the quasidark states |E 1D are more likely to decay by a process of the form |E 1D → |E 1± → |G, 0 , where the first step is due to dephasing, the second to photon emission. Hence, in such conditions JC-like features can become more prominent in the measured spectrum. At high dephasing (γ d = κ), the dark-state features are effectively invisible in the emission spectrum. The bi-excitation scenario is where the multilevel model truly shines, compared to the reference TC model: we can see in the bottom of Fig. 5 that multilevel and JC models are indeed in agreement on the location of both single-excitation (ω ω 0 ± g eff , gray dot-dashed lines) and bi-excitation [ω ω 0 ± ( √ 2 ± 1)g eff ], gray dotted lines) resonances, while the bi-excitation peaks of the TC model are visibly shifted. Also in this case the multilevel model presents an additional resonance (or perhaps a collection of closely spaced resonances) at ω ω 0 . Referring again to the notation of Section IV A, these can be ascribed to two processes: the photons emitted in the |E 1D → |G, 0 transition, which is associated with a small decay rate, and those emitted in the "second rung" process |E 2D → |E 1D , which we recall is very similar to photon emission by an empty cavity, and hence has a singnificantly higher decay rate κ. We again notice that the ω ω 0 resonance can be weakened in the presence of dephasing, however it remains prominent even in the "high dephasing" scenario. Our interpretation for this is the following: while dephasing of order κ can suppress the |E 1D → |G, 0 process, as discussed above, the |E 2D → |E 1D decay channel remains a probable process since it is also associated with a decay rate of order κ. It is worth pointing out that if one goes beyond our symmetric pump hypothesis (populating all the emitter excited states with equal probability -see Section III) it may of course be possible to reduce the value of p dark and hence the prominence of dark-state features in the cavity emission spectrum.
To summarize, the results of this subsection tell us that the multilevel model is a viable system to test the characteristic ∝ √ ng eff splitting of energy levels that is still the holy grail of experimental JC research. In this light, the main strength of the multilevel model is that the effective coupling g eff may be significantly boosted by considering emitters with many closely spaced sublevels.

B. Photon autocorrelation function
As a second step, we will look at the steady state behaviour of the zero-delay photon autocorrelation function, under weak coherent driving of the cavity. In all three models this quantity is defined as follows: where again we made use of Eq. (31) in order to obtain an expression that only features system operators. In detail, in this subsection we consider the master equation (8), setting Λ = 0, Ω 0, and vary the driving field frequency ω L . While the resulting master equation is explicitly time dependent, it is still possible to define a steady state by considering an interaction picture with respect to the Hamiltonian H L = ω Lâ †â +ω L n k=1 |e k e k |. This is a standard procedure in quantum optics and it amounts to replacing V drive → V (int) drive = Ω â +â † and H → H (int) = H − H L in the master equa-Single-excitation regime -κ = g Bi-excitation regime -κ = 0.05g FIG. 6. steady-state photon autocorrelation g (2) (0) for N = 3 emitter sublevels and a sub-level energy spread fixed at ε = 0.05g eff . We set the incoherent pump to zero and introduce a coherent cavity drive, corresponding to model parameters Λ = 0, Ω = 0.01g eff . Dephasing increases from left to right as before. Solid blue lines = multi-level model; dashed red lines = reference JC model; dot-dashed green lines = reference TC model; dotted black lines = reference coupled oscillators model. For details on all the auxiliary models used here see Appendix A tion (8). In turn, this procedure implicitly defines a time independent superoperator L (int) which does admit a steady state ρ (int) ∞ . Such a steady state is again found numerically with the same process described in the previous section, and it can be easily shown that it provides the correct expectation values involved in Eq. (34). The obtained autocorrelation functions for multilevel, JC and TC models are displayed in Fig. 6 as functions of the driving field detuning. We note that in this case no rescaling of the TC and JC plots is required, since the g 2 (0) is normalised by definition.
We observe that, excluding the region ω ω 0 , the multilevel model does a good job at matching the JC model predictions, capturing both the position and magnitude of resonances and dips. These observations are true both in the single-and bi-excitation regimes. In contrast, the g (2) (0) of the TC model displays features that are significantly less pronounced than their JC counterparts (recall that we are using a logscale), as well as being shifted in frequency in the biexcitation regime. This is compatible with the intuition that the TC model is "more linear" than both the JC and multilevel ones. Finally, the coupled oscillator model is of course even farther off from the JC predictions. As expected, the latter indeed displays features that are several orders of magnitude smaller than those in all the other models, plus it satisfies g (2) (0) > 1 at all times. As anticipated, however, the multilevel model displays additional features in the region ω ω 0 : at low or no dephasing a double-peak structure can be clearly seen, likely due to processes involving the quasi-dark states (we recall that for N = 3 there are two quasi-dark states for each value of the excitation number n ≥ 1). As dephasing is increased the prominence of these peaks is reduced, presumably due to the introduction of additional decay channels for the quasi-dark states (see also previous Subsection). However, dephasing also results in the appearance of a central peak in the g (2) (0), which may be due to an increased likelihood that a pair of driving photons excite a quasi-dark state of the system, and are later emitted at the frequency of an empty cavity by a process of the form |E 2D → |E 1D → |G, 0 . While a more rigorous analysis of these effects would be certainly worthwhile to confirm our intuitions, we feel that it goes beyond the scopes of this work.
As anticipated, the g (2) (0) is one of the most commonly adopted witnesses of nonlinearity in emitter-cavity systems such as the ones we are investigating. The results of this subsection thus confirm that the multilevel model closely matches the nonlinear signatures of JC physics, and in this respect it clearly outperforms the reference TC model with the same values of N and g eff .

VI. CONCLUSIONS
We explored a novel route to realise the strong coupling regime of light-matter interaction, exploiting a multi-level emitter coupled to a single-mode electromagnetic field. When the N excited sublevels of the emitter are nearly degenerate, we showed both analytically and numerically that our proposed system can closely approximate many aspects of JC physics. The associated light-matter coupling constant scales as g eff ∝ √ N and, crucially, the quantum nonlinearities of our model are not suppressed with increasing N. Indeed, we are able to observe the characteristic energy level splitting of order ∼ √ n, where n is the number of combined (emitter + field) excitations, as well as the associated nonlinear signatures in emission spectra and photon autocorrelation functions. Moreover, we have been able to characterise the main effects associated with the additional levels of the emitter: the appearance of very interesting sharp resonances in the middle of the emission spectrum, associated with the slow decay of "quasi dark states". While it may be possible to observe such resonances in an experiment, our analysis shows that these are quickly washed out in the presence of dephasing. Observing darkstate signatures would thus require cryogenic temperatures.
Our findings are particularly relevant in solid-state and plasmonic implementations of Cavity QED, where the emitters of choice are often multi-level quantum dots, and may pave the way towards the observation and exploitation of light-matter strong coupling at room temperature.

VII. ACKNOWLEDGMENTS
We thank M. G. Genoni for the useful discussion on modelling dephasing. Support by the Science Foundation Ireland (SFI) under grant 18/RP/6236 is gratefully acknowledged. TT acknowledges support from the University of Nottingham via a Nottingham Research Fellowship.