Unconventional saturation effects at intermediate drive in a lossy cavity coupled to few emitters

Recent technological advancements have enabled strong light-matter interaction in highly dissipative cavity-emitter systems. However, in these systems, which are well described by the Tavis-Cummings model, the considerable loss rates render the realization of many desirable nonlinear effects, such as saturation and photon blockade, problematic. Here we present another effect occurring within the Tavis-Cummings model: a nonlinear response of the cavity for resonant external driving of intermediate strength, which makes use of large cavity dissipation rates. In this regime, $(N+1)$-photon processes dominate when the cavity couples to $N$ emitters. We explore and characterize this effect in detail, and provide a picture of how the effect occurs due to destructive interference between the emitter ensemble and the external drive. We find that a central condition for the observed effect is large cooperativity, i.e., the product of the cavity and emitter decay rates is much smaller than the collective cavity-emitter interaction strength squared. Importantly, this condition does not require strong coupling. We also find an analytical expression for the critical drive strength at which the effect appears. Our results have potential for quantum state engineering, e.g., photon filtering, and could be used for the characterization of cavity-emitter systems where the number of emitters is unknown. In particular, our results open the way for investigations of unique quantum-optics applications in a variety of platforms that neither require high-quality cavities nor strong coupling.


I. INTRODUCTION
At the heart of quantum optics lies the interaction of light with matter at the level of individual quanta.As a result of the light-matter interaction between a single or an ensemble of two-level emitters and a resonant single-mode cavity, the emitters introduce nonlinearity to the otherwise linear cavity spectrum.This nonlinearity results in a splitting of eigenenergies known as the Jaynes-and Tavis-Cummings ladders [1,2].Under weak cavity and emitter excitation, one effect of these ladders is vacuum Rabi splitting in the spectrum of the system.Three other well-known quantum-optical effects also arise from this nonlinearity: saturation [3], photon blockade [4], and unconventional photon blockade [5][6][7].These effects are all of great interest for quantum control of light fields with important applications such as singlephoton switches [8][9][10][11][12] and transistors [9,11,13,14] and the generation of specific quantum states [15][16][17][18][19][20][21][22].In this paper, we demonstrate yet another effect, reminiscent of the saturation effect, which shows potential for applications in, e.g., quantum state engineering or the characterization of the number of quantum emitters in the cavity.
The saturation effect occurs when an emitter or nonlinear medium in a cavity cannot absorb more photons and thus has become saturated.In the spectrum, this is revealed as a merging of the vacuum Rabi doublet * therese.karmstrand@chalmers.se into a single Lorentzian peak at the cavity resonance when increasing the intracavity field [3].Ideally, a single photon incident on the system is needed to saturate a single emitter in the cavity.In this case, singlephoton saturation could implement, e.g., a single-photon transistor [13] or a single-photon sensor [23,24].Reference [25] demonstrates experimental progress approaching the single-photon limit.However, the saturation effect is typically associated with a very strong drive.The need for a strong drive is due to the generally low probability for photon-emitter interaction [26] and, in the many-emitters case, to the fact that the entire medium must be saturated [27].This is problematic for applications, especially if the systems exhibit large dissipation rates.
In photon blockade [4], on the other hand, the anharmonicity in the spectrum blocks the absorption of a subsequent photon.The effect occurs for a resonant drive on one of the polariton transitions.A characteristic of photon blockade is nonclassical photon-counting statistics, which can be probed via the normalized second-order correlation function g (2) in the weak drive regime [28][29][30].Two signatures of nonclassical light are photon antibunching [g (2) (τ ) > g (2) (0)] and sub-Poissonian photon statistics [g (2) (0) < 1] [31,32].Thus, photon blockade could be exploited for the generation of nonclassical photon states, e.g., a single-photon source.The single-photon blockade has been extensively explored theoretically [4,[33][34][35][36][37] as well as demonstrated experimentally [17,38,39].Stimulated by the potential for quantum state engineering including more than one photon [26], there have recently also been several works on the multi-photon blockade [40][41][42][43][44]. Additionally, a breakdown of the photon blockade has been observed for strong external driving [37,45,46], and has been studied with mean field theory in the limit of large quantum emitter numbers [47].One basic condition for both single-and multi-photon blockade is that the decay rates of the system should be much smaller than the cavity-emitter interaction strength.These conditions require high-quality cavities as well as small emitter dephasing.For that reason, demonstration of photon blockade in dissipative systems remains difficult.
An alternative approach to the generation of nonclassical states of light, exploiting the anharmonic Jaynesor Tavis-Cummings spectrum, is through the so-called unconventional photon blockade [5][6][7].In contrast to traditional photon blockade, the unconventional photon blockade effect relies on the interference between two transition pathways (see, e.g., Refs.[22,48,49]) when the drive is tuned in between the two polariton transitions.Thus, being an interference effect, the overlap due to broader transition linewidths can be exploited.Similar to photon blockade, unconventional photon blockade displays nonclassical photon statistics in g (2) measurements in the weak drive regime.Originally, unconventional photon blockade was found for coupled Kerr resonators [5,50].Since then, it has been predicted [22,48,51] and demonstrated experimentally [42,52] with dissipative cavity-emitter systems described by the driven Jaynes-and Tavis-Cummings Hamiltonians.It has also been predicted for large ensembles of emitters provided large enough individual cavity-emitter interaction strength [53,54].Nevertheless, demonstrating unconventional photon blockade remains difficult, due to fast oscillations of g (2) that exceed the resolution of state-of-the-art detectors and the requirement of fine-tuning of intrinsic system parameters [7].
In this work, we demonstrate a different approach to harness the nonlinearity introduced by one or a few twolevel emitters interacting with a dissipative cavity.Our scheme is simple, employing a continuous-wave (CW) coherent drive, requiring only detection of the steadystate cavity population.We base our analysis on numerical solutions of the corresponding Lindblad master equation.We use the rotating-wave approximation for the drive and coupling terms, but otherwise no further approximations that would limit us to the weak drive regime [48,53,54].Therefore, we can explore the intermediate drive regime, where we find a saturation-like effect on the cavity population, due to destructive interference between two excitation pathways.The cavity can be excited either directly by the drive or by the excited emitters.The interference between these two transition pathways has similarities with the interference that gives rise to unconventional photon blockade.Therefore, we name the effect observed here unconventional saturation.
The unconventional saturation effect is revealed in the cavity response to resonant driving of intermediate strength and arises due to the intermittent saturation of the destructive interference, leading to direct cavity excitation.Already visible in the weak-excitation regime, well before traditional saturation, the effect leads to a strong nonlinear dependence of the intracavity field on the drive strength.Moreover, it is not limited to strong cavity-emitter coupling.Instead, we find that the basic requirements for observing unconventional saturation are: few quantum emitters, large cooperativity C ≡ 4g 2 col /γ c γ e , and intermediate drive strengths.The second condition, large C, is naturally found in many lossy cavities where the cavity decay rate γ c is large compared to the emitter decay rate γ e , such that they fulfill γ c γ e 4g 2 col .Here, g col is the collective cavity-emitter interaction strength.In comparison to unconventional photon blockade, which often involves systems with small C ≈ 0.5−2 and weak driving, the unconventional saturation effect becomes notable for C 10 with intermediate drive strengths and grows more prominent for higher C.
The signature of unconventional saturation is the dominance of (N + 1)-photon processes in scattering from an N -emitter-cavity system.Somewhat hand-wavingly, the emitter ensemble can be seen as a saturable mirror, which only can reflect states with up to N photons.We identify the origin of this effect as the same type of quantum interference that explains unconventional photon blockade.Nevertheless, the fact that unconventional saturation can be detected in steady-state scattering could facilitate a more straightforward experimental demonstration than the more elaborate photon-correlation measurement g (2) needed for verifying unconventional photon blockade.Moreover, as opposed to vacuum Rabi splitting and photon-blockade techniques, our approach unambiguously differentiates between different numbers of emitters with the same collective interaction strength.This property makes it a promising scheme for characterization of cavity-emitter systems where the number of emitters is unknown, e.g., counting of NV-centers in diamond [55,56], localized emitters in hBN [57,58], or molecules in a Fabry-Pérot cavity [59], as well as for the verification of fundamental differences between singleand few-quantum-emitter systems.Other possible applications are within technologies such as quantum imaging [60], quantum metrology [61] and more, which rely on the generation of nonclassical light fields.
Besides being a -to our knowledge -novel quantumoptical effect, we also see the potential for the use of unconventional saturation for progressive quantum state engineering that could find a natural place in hybrid quantum systems, similar to the setup proposed in Ref. [13].The generation of specific quantum states of light with dissipative systems has already been proposed for other setups, including single- [62,63] and multiphoton [21,22,64] generation.We believe that our work offers a foundation for further explorations of hitherto unknown effects that could complement and improve existing schemes.Our results already suggest a form of photon filtering that can be achieved using a setup (CW drive and scattering) that is simpler compared to many other schemes.
A potential platform for demonstrating unconventional saturation is hybrid light-matter systems using, e.g., broad-linewidth surface plasmons and narrow-linewidth localized two-level emitters.Light-matter hybridization is a growing research field that utilizes hybridized states of light and matter such as surface plasmons as the carrier of the photonic component.Part of the attraction of such setups is the sub-wavelength confinement of the light mode that can greatly enhance the interaction with optical emitters [65][66][67][68].Even ultrastrong coupling [69,70] has been demonstrated [71] in these systems.The potential of hybrid systems for quantum technology has already been demonstrated in Ref. [13], proposing a singlephoton optical transistor.Lastly, the fact that strong coupling between cavities and emitter ensembles can be observed at room-temperature in these dissipative hybrid systems, motivates the search for observation and possible application of quantum-optical phenomena beyond cryogenic temperatures [72][73][74][75][76].
This article is organized as follows.Our theoretical framework is presented in Section II, including the driven Tavis-Cummings model (Section II A), the master equation used for numerical calculations (Section II B), scattering from the cavity (Section II C), and an analog classical coupled-oscillator model used for analytical calculations in the weak-drive regime (Section II D).In Section III, we present the results from our explorations of the driven Tavis-Cummings model, which show a saturation-like response of the cavity population in the intermediate drive regime.First, we study the spectrum in Section III A and note that there is a sweet spot for quantum effects on resonance, which differentiates between different system sizes.Thereafter, the system response is examined for resonant driving in the weak-to intermediate-drive regimes in Section III B. The observed nonlinear response is analyzed in Section III C in terms of an effective drive acting on the cavity.In Section III D, we build on this description to derive an analytical expression for the critical drive strength required for entering the nonlinear regime.We also find a figure of merit for the observed effect in Section III E and show that it can be explained by quantum interference effects in Section III G. Finally, the conclusions from our investigations are presented in Section IV.We give additional details for some calculations in appendices: Appendix A shows the mapping between the classical and quantum models used in Section II D, Appendix B reviews the quantum theory for a propagating laser beam, and Appendix E contains plots further demonstrating the effect of coupling strengths and cooperativity on the unconventional saturation.

II. THEORETICAL FRAMEWORK A. Coherently driven Tavis-Cummings model
Figure 1 shows a schematic illustration of the driven dissipative Tavis-Cummings system considered in this work.The Tavis-Cummings model describes the dynamics of an ensemble of N identical quantum emitters interacting with a common single-mode cavity field [2].No interaction between the individual emitters is included, which is motivated in circumstances where the cavityemitter interaction is the dominant interaction governing the dynamics.Including a coherent drive Ē cos(ω d t), with spatial amplitude Ē and frequency ω d , the driven Tavis-Cummings Hamiltonian can be written within the rotating-wave approximation (RWA) as Here, â and â † are annihilation and creation operators, respectively, for the cavity mode, ω c is the cavity frequency, σ−i and σ+i are the Pauli lowering and raising operators, respectively, for the ith quantum emitter, ω e is the transition frequency of the emitters, Ω d is the strength of the cavity drive, and g is the strength of the coupling between the cavity mode and a single quantum emitter.
The cavity drive strength is given by Ω d = (ā • Ē)/ , where the spatially dependent parameter ā is cavityspecific.Thus, the exact form of Ω d is determined by the explicit drive and cavity configuration.No external driving of the emitters is considered.This assumption is natural for most experimental setups where the emitters are located inside the cavity, but works as well for open cavities such as plasmonic nano-cavities that typically have much larger transition dipole moments than most quantum emitters.Furthermore, spatial variations of the cavity-emitter dipole interaction is neglected.Thus, we take g = μe • Ēc / for all emitters, with the transition dipole moment μe interacting with the cavity field with amplitude Ēc .This approximation is sufficient for many situations involving only a few localized quantum emitters, and in situations where the emitters are small compared to the cavity.With equal interaction rates g, the structure of the interaction term in the Tavis-Cummings Hamiltonian [Eq.(1)] leads to the collective interaction strength g col = √ N g between the cavity and the collective bright mode of the emitter ensemble.

B. Master equation
In this work, the scattering from the system under weak to intermediate driving is investigated.To solve for the cavity-emitter state including dissipation, an openquantum-system approach is employed, using the master equation Here the operator • acting on the density matrix ρ is the standard Lindblad superoperator for dissipation associated with the operator ô [77].With this master-equation approach, it is also possible to treat the case of strong driving, for which the traditional saturation effect would be found.
The first term in Eq. ( 2) describes coherent evolution with the Tavis-Cummings Hamiltonian.The second term describes radiative and non-radiative dissipation of the cavity mode, making the total dissipation rate γ c = γ r c + γ nr c .In the third term, the individual dissipation rates γ e for the emitters are assumed to be equal.The form of Eq. (2) neglects the contribution of thermal photons to the system dynamics and is therefore valid for low temperatures or high-frequency quantum systems with ω c , ω e k B T , such that thermal fluctuations do not particularly affect the dynamics.In experimental realisations, this condition is naturally met, e.g, for optical frequencies at room temperature.
A more compact way of writing Eq. ( 2) is in terms of the Liouvillian superoperator: Then, the task of finding the steady state is reduced to the eigenvalue problem with a Hermitian density operator ρss satisfying the normalisation condition C. Probing the cavity For applications in quantum photonics, the scattering from the system is of great interest.In cavity-emitter systems where the cavity interacts much more strongly with the environment, the collection of emitted photons from the emitters may be neglected.This complies with the condition γ r c γ e , which is what is considered in this work.Moreover, in most experimental setups, the collection of emitted photons from the driven system can be located such that the incident laser field is filtered out.The collected scattering S from the system will therefore be proportional to the radiative cavity decay rate and the average cavity population:

D. Analogue classical coupled oscillator model in weak drive regime
For adequately weak drive, much of the phenomenology associated with coupled cavity-emitter systems can be described by a classical coupled-oscillator (CO) model [78][79][80].Here a CO model will be used for comparison when analyzing the quantum effects that arise beyond the weak drive regime.
The CO model considered involves N + 1 mechanically coupled masses on springs.The corresponding classical coupling constant and drive strength are −2g √ m c m e ω c ω e and Ω d √ 2m c ω c , respectively.For simplicity, the cavity and emitter masses, m c and m e , are set to 1.The mapping of the quantum parameters to the classical model can be found Appendix A. Letting index 0 denote the oscillator representing the cavity mode and index 1, ..., N the emitters, the equations of motions for the classical analog of N identical emitters coupled to a coherently driven cavity mode are ẍi + γ e ẋi + ω 2 e x i + 2g The set of coupled equations ( 7)-( 8) is easily solved by making the ansatz z i = C i e iω d t for all i = 0, .., N and noting that x i = Re{z i } and cos (ω d t) = Re e iω d t .The solutions for the amplitudes are Equations ( 9) and ( 10) can be used to calculate the classical oscillator energies E c/e = 1 2 ω 2 c/e |C 0/i | 2 , which can be compared with the average energies E qm c = ω c â † â and E qm e,i = ω e σ+i σ−i in the cavity mode and emitter ensemble, respectively, calculated using the Tavis-Cummings model.Since the average energies in both models must be the same, the classical analogue to the populations is given by Equations ( 11) and ( 12) will be useful for comparing the classical and quantum results in this article.Note that Eq. ( 12) represents the total ensemble average population.

III. UNCONVENTIONAL SATURATION EFFECT AT RESONANT DRIVING
A. Scattering spectrum Large loss rates generally limit experimental investigations to weak excitation, â † â 1.In this regime, strong coupling with the emitter ensemble will lead to vacuum Rabi splitting in the spectrum.This effect can be seen in Fig. 2(a) for the steady-state cavity population â † â ss = Tr â † âρ ss under continuous driving within the weak-excitation regime, with emitters on resonance with the cavity (ω c = ω e ).The steady state ρss is found from the master equation by solving Eq. ( 4) for N = 1, 2, 3, and 4 quantum emitters, and is compared with the classical solution n c given by Eq. (11).As can be seen, there are only small differences in the spectra between different N and the classical solution.
On the other hand, examining the same spectra on a logarithmic scale in Fig. 2(b), large deviations (several orders of magnitude) from the classical model can be seen when the drive is resonant with the cavity and the emitters.Despite having the same collective interaction strength g col , large differences can also be seen between the spectra for different numbers of emitters in the ensemble.The spectrum for N = 1 shows the largest deviation from the classical case; adding more emitters yields spectra approaching the classical response.Thus, we have found a sweet spot for quantitative quantum effects that differentiate between different emitterensemble sizes N in the weak-excitation regime.In fact, it turns out that a strongly N -dependent nonlinear response can be accessed for resonant driving in the steady state, as will be shown in the next section.
For this simulation and throughout the main text, the emitters are taken to be on resonance with the cavity mode, i.e., ω e = ω c .The other parameters used in Fig.  were γ c /ω c = 0.03, γ e /ω c = 0.0003, g col /ω c = 0.03, and Ω d /g col = 0.25.

B. Mean cavity response for increasing drive rate
Encouraged by the visible quantum effects on resonance in the spectrum, we here further explore the optical response of the Tavis-Cummings model for resonant driving.The spectra in Fig. 2 are calculated with a drive strength that is often considered to be in the weak-drive regime.Nonetheless, it is perhaps more instructive to discuss in terms of an intermediate drive regime, Ω d < g col , prompted by the large losses that retain the system response in the weak-excitation regime.
In Fig. 3, we show the evolution of the steady-state cavity population â † â ss and the steady-state total ensemble population σ+ σ−  of drive strengths from truly weak-drive conditions to strong drive (Ω d > g col ).The populations are calculated for N = 1 − 4 and are compared to the corresponding classical populations n c and n ens given by Eqs.(11) and (12).Intriguingly, as the drive strength increases, the results in Fig. 3(a) show the cavity response entering into a nonlinear regime for each system size in turn.This effect is visible as an N -dependent break from the linear, classical response.Contrary to what naively could be expected in the weak-excitation regime, this implies that the few-level nature of a quantum-emitter ensemble plays an important role for the system dynamics, similar to the traditional saturation effect, already well before the system enters the strong drive regime.
More specifically, Fig. 3(a) shows how the cavity response transition between two linear regimes at weak and strong driving, respectively.In the weak-drive regime, the cavity population follows the classical population n c calculated with the CO model.Then, depending on the number of emitters N , a nonlinear regime is passed before the cavity again has a linear response described by an uncoupled driven damped harmonic oscillator with population n o c = Ω 2 d /γ 2 c .Comparing the cavity population in Fig. 3(a) with the ensemble population in Fig. 3(b) shows that the nonlinear regime appears well before the emitter ensemble has saturated, i.e., when σ+ σ− ens ss N/2.It is not until strong drive conditions have been reached, Ω d /g col > 1, that the ensemble saturates and the cavity population approaches the response of an uncoupled driven harmonic oscillator.
A study of the slopes for the cavity population in the intermediate drive regime [Fig.3(a)] shows the first evidence of the emitter-ensemble origin of the strongly nonlinear behavior.In linear response, the cavity population is expected to have a linear dependence on the driving intensity I ∝ Ω 2 d .This dependence is precisely what Fig. 3(a) shows for a sufficiently weak drive, where all systems have a slope of two in the log-log plot.If a multiphoton process of order n is dominant, on the other hand, the cavity population would instead be proportional to the nth power of the driving intensity: â † â ss ∝ I n ∝ Ω 2n d .Simple log-fits of the new slopes in Fig. 3(a) give the inclines ∼ 4, 6, 8, and 10 which would correspond to 2-, 3-, 4-, and 5-photon processes for the cases of N = 1, 2, 3 and 4 respectively.This indicates the dominance of (N + 1)-photon processes facilitated by the ensemble of N emitters.

C. Multiphoton processes due to unconventional saturation effects
To gain a better understanding of the observed dynamics, we study the system for weak external driving.By applying appropriate weak-drive approximations, the master equation in Eq. ( 2) simplifies to the same type of coupled equations of motion for the density-matrix elements [81] as in the coupled-oscillator model described in Section II D. Therefore, more insight into the Tavis-Cummings dynamics in the weak-drive regime can be gained by observing the simple analytical solutions to the classical equations of motion presented in Eqs.(7) and (8).
This section includes two parts.In Section III C 1, we show that the observed cancellation of cavity population in the weak-drive regime can be understood as a destructive interference between the external drive and the effective driving from the emitter ensemble due to the cavity-emitter coupling.Specifically, we formulate an effective drive on the cavity inspired by the classical equations of motion, which explain the cavity response in the weak-drive regime well.In Section III C 2, we discuss how the breakdown of the destructive interference leads to the dominance of (N + 1)-photon processes in the cavity response at intermediate drive strengths.We also provide a simple phenomenological picture of the cavity response.In this picture, we neglect the ensemble and describe the observed (N + 1)-photon process by multiphoton absorption events followed by cavity decay through single-photon processes.

Destructive interference
In the coupled-oscillator model, an effective drive on the cavity can be defined by combining the external drive term with the coupling to the emitters, i.e. rearranging the terms in Eq. ( 7): The solution for x i is the real part of the ansatz z i = C i e iω d t with the coefficient C i given by Eq. ( 10).Inserting the solution for x i on resonance (ω d = ω c = ω e ≡ ω 0 ) gives the effective drive Equation (14) shows that the external drive and the ensemble oscillator will interfere destructively.For γ c γ e g 2 col , the effective amplitude is much smaller than Ω d .This result explains the deep dip observed at resonance in the spectrum presented in Fig. 2(b).Moreover, it elucidates the suppression of the coupled cavity population n c compared to the uncoupled cavity population n o c shown in Fig. 3(a).A similar effective drive can be found in the Heisenberg picture for coupled quantum oscillators described by the position and momentum operators {x c , xe,i , pc , pe,i }.The effective quantum drive has the same form as Eq. ( 13), but with the classical position variable x i replaced by the quantum operator xe,i .Since the eigenenergy spectrum of the emitter ensemble resembles a harmonic oscillator up to the same order of excitations as the number of emitters N , the validity of this coupled-oscillator picture is motivated.As such, the emitter ensemble behaves like a harmonic oscillator for weak excitation where higherorder terms are negligible.By employing this coupledoscillator picture in the weak-excitation regime, a quantum analog to the classical effective drive described above can be formalized utilizing the properties of coherent states.
A coherently driven damped harmonic oscillator will also be in a coherent state.Hence, we can make a coherent-state approximation of the emitter ensemble to order N in the weak excitation regime.In terms of Fock states, this coherent state can be written as The complex amplitude α ens with |α ens | 2 ∝ σ+ σ− ens ss is defined by the emitter ensemble.
For this work, a rectangular 'time-bin' temporal mode with duration T is a sufficient description for the mode of the ensemble state in the weak-drive regime.This mode choice gives a simple expression for the complex amplitude: The time duration T is a characteristic timescale set by the system.For the considered Tavis-Cummings system, the occupation of the ensemble is related to the cavity field through the collective coupling g col .Therefore, the natural choice of T for the ensemble state is To formulate a quantum analog to the classical effective drive, we have to relate the approximate coherent state for the ensemble to the classical external drive on the cavity.The relation can be found by considering an idealized laser for the external drive.The state of an idealized laser beam propagating through free space can be represented as a continuous-mode coherent state.This continuous-mode coherent state can be partitioned into an infinite set of discrete-mode coherent states with amplitude α d .The partitioning into discrete temporal modes can be performed with a large freedom of choice under the condition that the characteristic mode timescale T d and oscillator frequency ω d obey ω d T d 1.This freedom of mode choice is discussed in detail in Appendix B, where we also give the partitioning into rectangular time bins as a specific example.Hence, in accordance with Appendix B, we can also choose rectangular temporal modes for the external drive with a time duration T d , which can be chosen arbitrarily as long as the condition ω d T d 1 is fulfilled.To compare the effects of the two sources of driving, the same choice of time duration must be made for the ensemble and the external drive.Thus, we take T d = T , which gives the coherent state amplitude with T = 1/g col as given by Eq. ( 18) above.
Taken together with the coherent state approximation for the ensemble and the mode-matched partitioning of the laser beam, Eq. ( 13) suggests that the cavity can be seen as driven by an effective coherent state in the linear regime.In terms of Fock states, this effective drive state can be written down as for small |α d | 2 and |α ens | 2 .For sufficiently weak drive, only lower-order Fock states (n ≤ N ) contribute notably to the scattering dynamics.In this regime, Eq. ( 15) gives a classically derived analytical expression for the effective drive Ω eff , which agrees well with the numerical calculations using the master equation.By combining this picture of classical destructive interference with the idea of two modematched coherent states, it can be seen that Eq. ( 21) describes an effective, coherent drive on the cavity with amplitude where T = 1/g col is the characteristic timescale identified for the ensemble given in Eq. ( 18).The validity of this choice of T is confirmed by the good agreement between our analytical predictions using Eq. ( 18) and the exact numerical calculations in the weak-drive regime presented in Fig. 3.

Breakdown of destructive interference at intermediate drive strengths
Under the weak drive conditions discussed so far, the classical and quantum models give the same result for the cavity field due to external driving.Nevertheless, there are distinct differences between the classical and quantum models that become clear in the effective-drive picture.In the classical picture, the effective drive in Eq. ( 14) represents two harmonic fields acting on the cavity with opposite phase, which therefore cancels.Since both fields are harmonic, the relative amplitude Ω eff , written in Eq. ( 15), will not change when the drive strength increases.Thus, the classical cavity population n c maintains a linear dependence on the external drive in the strong-drive regime (Ω d > g col ).
The emitter ensemble, on the other hand, only resembles a harmonic oscillator up to photon number N .This truncation of the harmonic spectrum is reflected in the effective coherent drive in Eq. (21), where the summation only goes to N .Hence, the destructive interference between the external drive and the ensemble breaks down at order (N + 1).However, the breakdown of the destructive interference is only visible when the (N + 1)th state in the Fock-state expansion in Eq. ( 19) has become significant.In Fig. 3(a), the cavity response reveals this effect as (N + 1)-photon processes for intermediate drive strengths.For even stronger drive, the emitter ensemble will saturate, and the uncoupled cavity response n o c will be approached.
The (N + 1)-photon processes in the intermediate drive regime can be seen as a result of an unconventional saturation effect where the emitter ensemble intermittently saturates on the cavity-emitter interaction timescale identified above.This unconventional saturation is not visible in the spectrum since the emitter ensemble is still weakly populated and has not saturated in the traditional sense, i.e., can not absorb more energy.Instead, the unconventional saturation effect can be described as the destructive interference between different excitation pathways which occurs on the characteristic timescale T given in Eq. (18).
In a simplified picture, the unconventional saturation effect can be understood as a sequence of (N + 1)-photon pulses driving the cavity due to the intermittent saturation of the destructive interference.In Fig. 4(a), we provide a naive sketch illustrating one cycle of this (N + 1)photon process.The portrayed dynamics contain three distinct parts: (i) cancellation, (ii) (N + 1)−photon absorption, and (iii) exponential decay.First, the cavity is mainly in the ground state with an average population nc weak due to the emitter-drive interference, which cancels the cavity population.However, when (N + 1) photons arrive from the drive, the destructive interference breaks down.On the timescale of the cavity-emitter interaction, the emitters intermittently saturate, which leads to direct absorption of the (N + 1)-photon state in the cavity.Following the absorption event is exponential decay, where the (N + 1)-photons leak out of the cavity photon-by-photon.
The dynamics illustrated in Fig. 4(a) can be modeled with a simple phenomenological master equation for the probabilities P n (t) of occupying the nth Fock state in the cavity.Despite its simplicity, considering only cavity processes, this master equation qualitatively captures the unconventional saturation effect for drive strengths approaching the collective coupling g col .The details of this approach are shown in Appendix C.Here we present the main results.
Most importantly, the phenomenological master equation presented in Appendix C provides analytical solutions to the time-dependent probabilities P n (t).These solutions allow us to make an analytical prediction of the steady-state cavity population by calculating the timeaveraged contribution from having a stream of (N + 1)photon pulses driving the cavity due to the unconventional saturation effect.The steady-state cavity population in this naive picture can be found as where nc weak is the suppressed cavity population due to the coupling to the ensemble and P T N +1 is the probability of having (N +1) photons in the external drive during the time T given by Eq. ( 18).Here we have also introduced the notation nc ss = â † â ss for the steady-state cavity population.
An expression for nc weak can be found by employing the coupled oscillator model explained above.In Fig. 3(a), we have already seen that the steady-state cavity population is well described by the classical result n c for a sufficiently weak drive.Therefore, Eqs. ( 9) and (11) give an expression for nc weak which accurately predicts the cavity population in the weak-drive regime.On resonance, this expression is The probability P T N +1 is given by the Poisson distribution for the external drive with discrete-mode amplitude α d = Ω d T , Figure 4(b) shows the steady-state cavity population for N = 1−4 emitters obtained with Eq. ( 23) (red dashed curves).The analytical results are compared with the exact numerical calculations using the master equation in Eq. (2) (blue solid curves).The comparison shows that the simplified dynamics presented in Fig. 4(a), leading to the analytical prediction for nc ss given in Eq. ( 23), capture the unconventional saturation effect qualitatively.It can be seen that Eq. ( 23) accurately captures both the weak-drive behavior and the dominance of (N +1)-photon processes in the intermediate drive regime.The simple picture of (N + 1)-photon pulses arising from the breakdown of destructive interference at order (N + 1) can thus give a qualitative intuition for the unconventional saturation effect.
The analytical prediction overestimates the photon number.However, this is not surprising since the phenomenological model employed to derive Eq. ( 23) entirely neglects all effects from the coupling to the ensemble beyond the cancellation effect.For example, the possibility of excitation transfer to the ensemble during the decay process is completely overlooked.Nevertheless, the simplified model described in this section is a good tool that can be used to gain an intuition about the (N + 1)-photon processes associated with the unconventional saturation effect.
Equation (23) can also be extended to account for higher-order photon pulses (n > N + 1).In that case, the last term becomes a sum of the contributions.See Appendix C for details.Since the external drive is coherent, the probabilities P T n for the higher-order photon states (n ≥ N + 1) follow a Poisson distribution.Hence, the (N +1)-photon pulses contribute the most to the cavity response at intermediate drive strengths (Ω d < g col ).Including higher-order photon pulses in this simplified picture will, therefore, not qualitatively change the cavity response in this regime.In Appendix C, we show a calculation including the contribution of photon pulses up to order N + 5, which confirms our argument above.

D. Critical drive strength
The effective-drive picture, and the breakdown of the destructive interference discussed above, can also be used to write down a condition for entering into the nonlinear regime.Under the assumption that the emitter ensemble behaves as a driven harmonic oscillator up to order N , we would expect the cavity response to enter the nonlinear regime when the missing (N + 1)th term in the coherentstate approximation |α ens becomes comparable to the cavity population.From the coupled-oscillator perspective, this condition is easy to understand.That is, up to N excitations, the system behaves classically, and the ensemble can interfere destructively to cancel out excitation of the cavity.On the other hand, when the ensemble fails to interfere destructively due to its few-level spectrum, the cavity population becomes comparable to the missing (N +1)th term in |α ens .Formally, this condition for the critical drive can be written down as nc ss = (N + 1)P αens (N + 1).( 26) is the Poisson probability distribution for finding n excitations in the coherent state when |α| 2  1 and the factor (N +1) comes from having (N +1) excitations with probability P α (n).
Since we are approaching the nonlinear regime from the weak-drive regime, we can take nc ss = nc weak and use the expression for nc weak given in Eq. (24).We can also find an expression for the effective ensemble drive amplitude α ens = Ω ens T using Ω eff in Eq. ( 15): The condition in Eq. ( 26) thus becomes which gives the expression for the critical drive strength Ω cr that indicates the onset of the unconventional saturation regime.
As can be seen in Fig. 5, the calculation of the critical drive with Eq. ( 29) predicts very well the onset of the nonlinear regime.Figures 5(a)-(d) show the steadystate cavity response and the calculated Ω cr for a wide range of emitter decay rates γ e , whereas Figs.5(e)-(h) show the same for several different coupling strengths g col .In all panels, the analytically calculated Ω cr (N ),  collective interaction strengths g col .The remaining parameters were held fixed, using the same values as in Fig. 3.The red stars mark the derived analytical expression for the critical drive Ωcr.As can be seen, the analytical results predict exceedingly well the onset of the nonlinear regime.
marked with red stars, lie very close to the beginning of the nonlinear regime.Thus, we have not only found a new intriguing regime for performing quantum nonlinearity measurements, but we can also, with high accuracy, predict its onset for a wide range of parameters.

E. Figure of merit
So far, we have discussed unconventional saturation as an effect of the destructive interference (or the competition of) two distinct excitation pathways: cavity-drive and cavity-emitter.That interesting quantum effects can arise in dissipative Tavis-Cummings type systems due to quantum interference has already been shown with the so-called unconventional photon blockade effect [7,22,42,48,[52][53][54].The unconventional photon blockade, however, is observed for weak excitation and strong-coupling conditions.The unconventional saturation effect observed here, on the other hand, is instead present for intermediate drive strengths and appears for resonant driving in a parameter regime where unconventional photon blockade is absent.
It turns out that a good figure of merit for unconven-tional saturation is the cooperativity Why Eq. ( 30) is a good figure of merit can be seen by studying the induced transparency (reduced cavity population) on resonance due to the effective drive in the classical case.By taking the ratio between the classically derived coupled-and uncoupled-cavity populations, it can be seen that the cavity response will be suppressed with a factor depending on C: To arrive at this relation, we have used the expression in Eq. ( 24) for n c , found by solving the classical coupledoscillator equations of motion.A similar calculation for an uncoupled cavity, driven by the same external drive on resonance gives as already mentioned above.In Fig. 3(a), it can be seen that Eq. ( 31) governs the region in which the unconventional saturation effect can be observed.For small cooperativities, i.e., C ∼ 1 or smaller, the suppression of the cavity response due to the interaction with the emitter ensemble is too small for observing unconventional saturation.However, for C 10 the unconventional saturation effect starts to become clearly visible, and (as would be expected) it grows more distinct for increasing C. In Appendix E, additional simulations that show how the cavity response changes with the cooperativity can be found.

F. Suppression of the cavity response
Equation (31) in Section III E shows a classically derived expression for the suppression of the coupled-cavity response in the linear weak-drive regime.This result, together with the effective drive |α eff found in Eq. ( 21), underlines the expectation of a transition in the cavity response from a coupled coherent state to an uncoupled coherent state as the drive is increased.And indeed, this is what we see in Fig. 3(a).In Section III C, we also identified the timescale T = 1/g col , which defines the characteristic time for unconventional saturation.This timescale provides a simple relationship between the cavity decay time T c = 1/γ c and the effective (suppressed) drive α eff , which is easily found by rewriting the expression for n c in terms of the effective drive amplitude Ω eff in Eq. ( 15) and employing the relations for α ens in Eq. (22): The arguments above explain the two asymptotical behaviours observed and demonstrate the competition of timescales causing the unconventional saturation effect.In the following section, we show that the populations of the reduced density matrix for the cavity lie between the Poissonnian distributions for the two asymptotical coherent states.Moreover, we discuss how the destructive interference with the N -emitter ensemble appears in the populations and how they evolve as the drive increases.

Suppression of populations
As already noted, the nonlinear cavity response will appear in the transition between the response of a coupled cavity and that of an uncoupled cavity.In the two linear regimes, at weak and strong drive, respectively, the cavity will be described by the two coherent states with amplitudes α c and α o c .To find the amplitudes, we can use the well-known property |α| 2 = n for coherent states and use the classical results derived above.Thus, we find for the coupled coherent state and for the uncoupled coherent state.
To gain information on the state inside the cavity, we examine the populations ρ n ≡ n|ρ c |n of the reduced density matrix ρc ≡ Tr ens {ρ ss }.In Fig. 6, the cavity populations for N = 1, 2, 3, and 4 emitters (symbols) are plotted against the Poisson distributions for the coherent states with amplitudes α c (dashed line) and α o c (solid line) for three different drive strengths as indicated (dashed vertical lines) in Fig. 6(b).
Because of the destructive interference, we expect the populations to approach the coupled Poisson distribution, in the weak-drive regime.However, due to the truncation of the ensemble spectra at N excitations, the higher-order photon-state 'tail' will be pulled towards the uncoupled Poisson distribution.This effect can be viewed as a sequential flow of excitation from the ρ n≥N +1 populations, that are not affected by the cancellation, to the ρ n<N +1 populations via decay processes.
When the system is not in the unconventional saturation regime, we expect the distribution of the populations to be 'Poisson-like' (ρ 0 ρ 1 ρ 2 ...) except for the truncation-induced cross-over explained above.On the other hand, when the system is in the unconventional saturation regime, the breakdown of the destructive interference at order N + 1 will facilitate direct (N + 1)photon excitation of the cavity; see Fig. 4(a).Following the absorption is cavity decay, which is again a singlephoton process.Thus, in the unconventional saturation regime, we expect the ρ N +1 population to be comparable to the lower-order populations ρ n<N +1 , that is, (N + 1)ρ N +1 ≈ N ρ N ≈ ... ≈ ρ 1 .This picture agrees well with the breakdown of destructive interference discussed in Section III C 2. In Appendix C, we additionally show that the simple phenomenological master equation discussed in Section III C 2 also predicts the populations qualitatively.

G. Relation to exciton-induced transparency
The underlying destructive interference effect, giving rise to the observed suppression of the cavity population on resonance, can in the weak drive regime be understood as a classical analogue to electromagnetically induce transparency (EIT) [80].Indeed, the coupled set of equations presented in Eqs. ( 7) and ( 8) is the same as those used for modeling classical EIT.EIT is more commonly discussed in systems with more moderate lightmatter coupling than illustrated in Fig. 2. In Fig. 7, the cavity population spectrum is plotted for N = 1, 2, and 3 emitters together with the classical analogue for increasing interaction strength g col .At the top, in Fig. 7(a) and (b), the typical EIT regime with a rather moderate interaction strength g col < γ c /2 is shown.As is clearly visible throughout the log-log plots in the right panel, the Ndependent quantum interference effect at resonance persists for a wide range of interaction strengths.
For a system with one emitter, a semiclassical analog to classical EIT, treating the emitter quantum mechanically and the cavity classically, has already been referred to as 'exciton-induced transparency' (ExIT).ExIT has been discussed within the context of plasmon-exciton coupling [51,82,83] and demonstrated with a plasmonic nano-cavity and quantum-dot system [84].The saturation of the ExIT effect has been discussed in Ref. [85], although the exact power dependence was not analyzed.The regime with dominant two-photon processes for this single emitter system was thus not identified.However, the systems exhibiting ExIT should indeed be suitable for demonstrating unconventional saturation.

H. Relation to dressed-state polarisation
In Fig. 3(a) two critical points can be identified.The first point, at weak drive, marks the onset of the unconventional saturation effect, and it depends strongly on the emitter number, as discussed in Section III D. The second point, at strong drive, marks the conventional saturation of the emitter ensemble, as can be evidenced by also studying the steady-state ensemble population in Fig. 3(b).This point is related to the collective coupling strength.Thus, it is primarily the same for all system sizes with the same collective coupling.
In a resonantly driven Jaynes-Cummings model (N = 1), the critical point at strong drive has additionally been identified as a first-order phase transition for a strongly coupled system [37,45,46].Due to the strong cavityemitter coupling, the Jaynes-Cummings system can be seen as an atom-cavity 'molecule' that becomes dressed by the resonant external drive.Beyond the critical drive strength, which in our parameters translates to Ω d /g col = 1, the system exhibits phase bimodality in the steady state.This effect is known as dressed-state polarisation and has been analysed in detail in Refs.[45,46].The dressed-state polarisation effect has also been connected to the breakdown of the photon-blockade effect [37], and it has been studied with mean-field theory in the limit N → ∞ [47].
Unlike the unconventional saturation effect, the dressed-state polarisation effect is a strong-coupling effect in the strong drive regime, that relies on climbing the dressed Jaynes-Cummings ladder.The unconventional saturation effect, on the other hand, is present in, but not limited to, the strong-coupling regime.Instead, as discussed in Section III E, the relevant figure of merit for unconventional saturation is large cooperativity.By virtue of the single-emitter decay rate in the denominator in Eq. ( 30), a large cooperativity can be obtained even though the cavity loss is large (i.e., weak coupling) if the emitter decay is small.Thus, unconventional saturation can be present in systems which are prohibited from displaying dressed-state polarisation in the strong drive regime due to the large cavity dissipation.Moreover, in contrast to dressed-state polarisation, the unconventional saturation effect appear at intermediate drive strengths, well before the critical point, Ω d /g col = 1, marking the phase transition to the dressed-state polarisation phase in the resonantly driven Janyes-Cummings model.

I. Limitations of the model
With the aim of formulating the simplest model that predicts the main physical phenomena arising with small emitter ensembles coupled to a cavity in an environment, the master equation employed in this work is derived under two important conditions: (i) the cavity and emitters interact with the environment independently of the cavity-emitter coupling, and (ii) each emitter interacts with local environments.This approach follows a criterion of simplicity and is appealing since it keeps the discussion of the physics simple.Below, we discuss the limitations of the conditions (i) and (ii).Moreover, we also discuss how small changes to these conditions do not change the effect presented in this work.
In the following, we refer to the eigenstates of the coupled cavity-emitter system as dressed states, and we refer to the eigenstates of the cavity and emitter subsystems as bare states.
, cavity, respectively.As can be seen for all system sizes N = 1-4, the emitter ensemble behaves as a true coherent state in the coupled system only up to first order.Then, for increasing drive strength, each system saturates in turn, visible as the higher-order photon states approach the uncoupled-cavity distribution.The red circles illustrate when ρN+1 have become comparable to ρ1, i.e., ρ1 ∼ (N + 1)ρN+1 , which is a signature of the unconventional saturation regime.As is clearly seen throughout the log-log plots, the quantum interference effect is present for a wide range of interaction strengths below and at the strong-coupling limit g col ≥ γc/2.

The rotating-wave approximation
Condition (i) states that the dissipative part of the master equation is derived without taking into account the cavity-emitter interaction.The dissipators [D â and D σ− in Eq. ( 2)] obtained through this method will induce transitions between the bare cavity and emitter states.The master equation obtained in this fashion can accurately describe many cavity-QED experiments in regimes where the RWA can safely be applied.
Generally, the RWA is justified for systems with g col /ω c < 0.1, which is the conventional limit for the ultra-strong coupling (USC) regime [69].In this article, we demonstrate an unconventional saturation effect with a set of parameters giving g col /ω c = 0.03, which is clearly below the limit of USC.Moreover, as discussed in Section III E, the effect is also present for lower g col provided that the cooperativity C is large enough.Thus, the RWA should, in general, be justified in the regimes where the unconventional saturation effect typically appears.
However, the limit g col /ω c < 0.1 is merely a historical convention, not a strict boundary for when the RWA can safely be used.Therefore, when a high degree of accuracy is demanded, e.g., in the case of low photon numbers, the validity of the RWA may have to be reevaluated.When the RWA is not applied, the system Hamiltonian is no longer excitation-number conserving.The mixing of the bare ground state with higher-order states leads to virtual excitations in the dressed ground state (see, e.g., Ref. [69] for a review on USC discussing virtual excitations in the ground state).For cavity-emitter systems that are not ultra-strongly coupled, the effects from the mixing of states containing different numbers of excitations are typically negligible, but the effect demonstrated in this manuscript appears at very low excitation numbers.The change of the ground state when the RWA is not applied could therefore be important.
The retraction of the RWA requires proper adjustments to the master equation in Eq. ( 2), which would otherwise give unphysical excitation out of the coupled ground state even without driving at zero temperature.In Appendix D, we provide a master equation that induces transitions between the dressed states and accurately brings the system to the dressed ground state, including the counter-rotating terms in the Hamiltonian.Employing this master equation, we find that the counter-rotating terms do not qualitatively affect the unconventional saturation effect.For a comparison of the results obtained with the different master equations, see Fig. 10 in Appendix D.

Local emitter environments
The second condition (ii) neglects any effects of collective interaction between the emitter ensemble and the environment.For example, in optical bistability, which also appears within the drive-dissipative Tavis-Cummings model, the individual dissipation approach is the textbook procedure [86,87].Optical bistability was also recently studied with the individual-emitter-decay master equation in a low-photon-density regime [88] similar to ours.Thus, individual emitter decay seems to be an adequate method to model the central physical phenomena in the driven-dissipative Tavis-Cummings model.The connections between optical bistability and unconventional saturation in this regime also remain an interesting question for future study.
Nevertheless, the work of Dicke [89] highlighted the importance of collective dissipation for atoms contained in a small volume compared to the wavelength, which interacts with a single mode in free space.In that setting, the collectiveness of the interaction gives rise to the wellknown superradiance effect [90].Since then, the concept of collective decay has been transferred to cavity-QED systems, where quantum emitters are confined to a small volume inside the cavity.
In this work, all emitters interact uniformly with the cavity mode.This approximation is strictly only valid for emitter ensembles that are localized in a small volume compared to the mode wavelength, but serves as a good first approximation for many other experimental configurations.Because of small emitter confinement, the collectiveness of the interaction with the environment could be non-negligible for accurately describing the system dynamics.However, with the current set of parameters that show the unconventional saturation effect, the interaction between the emitters and the environment is weak and much smaller than the cavity dissipation (γ e γ c ).Therefore, the main decay channel for the emitters is through the cavity.Under these circumstances, where emitter emission is a rare event mainly shielded by the cavity, it is not unreasonable to assume that the remaining dissipation from the emitters will appear through local decay channels.Another situation where individual emitter dissipation would be a good approximation is for localized emitters, e.g., on a substrate.
Another important aspect is the treatment of only a few emitters.In the few-emitter regime, the difference between collective and individual dissipation for the emitters is minimal.So, a collective-dissipation picture does not qualitatively change the unconventional saturation effect; see Fig. 10 in Appendix D. The exact limits to the validity of individual versus collective dissipation is an interesting question that would require further theoretical work.

IV. CONCLUSIONS
In this work, the stationary response from a coherently driven cavity coupled to an ensemble of N quantum emitters, described by the Tavis-Cummings model, has been studied.The steady-state density matrix was calculated numerically using a master-equation approach without making the frequently applied weak-drive approximation.
Additionally, a classical coupled-oscillator model was applied to give analytical insight into the dynamics in the linear regimes.
For resonant drive frequency and intermediate drive strength, our results show strongly N -dependent nonlinear scattering.Specifically, we see the dominance of (N + 1)-photon processes in the nonlinear regime of the cavity response when it couples to an ensemble of size N .In contrast to observing Rabi splitting in the spectrum, this effect clearly differentiates between different ensemble sizes N with the same collective interaction strength g col .
Exploiting analytical results from a classical coupledoscillator model, and properties of coherent states, we found that the origin of this effect could be explained by the destructive interference between the ensemble and the coherent drive up to the order N .Thus, the ensemble behaves as a saturable mirror that can only reflect photon states up to order N .This unconventional saturation effect occurs due to a competition of interaction rates and arises for weak ensemble population, well before traditional saturation.We also derived an analytical expression for the critical drive Ω cr that to good accuracy predicts the onset of the nonlinear regime.Moreover, we find that a basic condition for the observed unconventional saturation effect is large cooperativity C.This condition can be met without the requirement of strong coupling, if the decay rates of the emitters are not too large.
The observed effect implies a simple continuous-wave method that could characterize dissipative cavity-emitter systems where the number of quantum emitters is unknown.The N -dependent interference effect and the resulting (N + 1)-photon processes in the cavity response could also be exploited for photon filtering.Thus, our results show great promise for the use of dissipative cavityfew-emitter systems for quantum state engineering.For this, further theoretical work that investigates the specific output state for different input states would be of interest.
Lastly, we note that a first-order expansion of the Holstein-Primakoff transformation (HPT) [91] for weakly excited two-level emitters would give the same coupledoscillator model as employed in this manuscript.Thus, the investigation into higher-order expansions of the HPT could be an interesting continuation of the present work.
The classical model involves N + 1 coupled oscillators when N emitters interact with the cavity mode.If all emitters are assumed identical, we may take ω i = ω, m i = m, k i = k for all emitters with index i = 1, 2, ..., N .Letting index 0 denote the cavity oscillator, the corresponding classical Hamiltonian is written In the classical model, no rotating-wave approximation (RWA) is applied.To compare with the quantum model, we therefore start with the quantum Rabi Hamiltonian with N identical quantum emitters, which simplifies to the Tavis-Cummings Hamiltonian when applying the RWA.In a weak-drive approximation, we may replace σ−i → bi , σ+i → b † i with bi and b † i being the annihilation and creation operators, respectively, of the ith quantum oscillator.Then, introducing quantum position and momentum operators {x c , pc , xe,i , pe,i }, the annihilation operators can be written Using these relations, we can write the quantum Hamiltonian in the xp-representation:

Coherent drive terms
Finding the classical analog for the strength of an external drive is done in the same manner by comparing the interaction Hamiltonians.Considering a semiclassical model where the quantized cavity mode is driven via dipole interaction with a classical coherent drive Ē(t) cos(ω d t), the quantum and classical interaction Hamiltonians are In Eq. (A8), the dipole interaction is rewritten in terms of the transition dipole moment for the cavity oscillator, which in the single-excitation manifold is μc = q c 2mcωc êk , where êk is a unit vector along the polarization direction.Then, defining the drive amplitude Ω d ≡ μc Ē(t) , we can write Appendix B: Quantum theory for the propagating laser beam Inside the laser cavity, the light field is well defined and can readily be quantized.In typical experiments, on the other hand, a beam of light has to travel through free space over distances that are much longer than the characteristic length scales of the studied system.The quantum theory for such a propagating light beam traveling in a straight line in free space can be found in, e.g., Chapter 6 of Ref. [29] and will be presented briefly below.
Consider a single propagating laser beam under circumstances where transversal effects are irrelevant to the experiment.Then, the quantization geometry can be taken as a finite cross-sectional area A (defined by the experiment) perpendicular to the propagation axis and a quantization axis of infinite length parallel to the propagation axis.This geometry corresponds to a onedimensional continuous-mode variable that can be taken as the frequency ω k with a mode spacing ∆ω = 2πc/L that goes to zero as the quantization length L tends to infinity, L → ∞.In this limit, the conversion from sum to integral is and the discrete Kronecker delta is related to a continuous Dirac delta-distribution as It follows that the continuous-mode annihilation and creation operators are related to the discrete operators as which fulfil the continuous-mode commutation relation, Under the assumption of a narrow bandwidth laser, i.e., the excitation bandwidth is much smaller than its central frequency, the lower integration bound can be extended from 0 to −∞ to cover the entire frequency axis in the integrals above.Thus, the corresponding timedomain operators are obtained as the Fourier transform of â † (ω) and â(ω): which have the commutation relation, The quantized continuous-mode electromagnetic field operator can then be written as where h.c.denotes Hermitian conjugate.The state of the laser field inside the laser cavity can be taken as a coherent state [92].The output field from the laser will then be a one-dimensional continuous-mode coherent state due to the lack of confinement along the propagation axis.Such a state can be represented using the Fock-space basis kets {|n } and is created from the vacuum state |0 with the continuous-mode annihilation and creation operators according to |{α(t)} = e dωα(ω)â † (ω)−α * (ω)â(ω) |0 . (B10) Here α(ω) is the continuous-mode spectral amplitude.The corresponding time-domain state with wavepacket amplitude α(t) is found via Fourier transform as The coherent-state mode functions satisfy the normalization condition where n is the number operator n = dωâ † (ω)â(ω) = dtâ † (t)â(t ).(B13) 1. Idealized continuous-wave laser In the following, a continuous-wave single-mode laser in a coherent state as described above will be considered.In typical optical experiments, the linewidth of the laser mode is much narrower than the other components of the observed quantum system.Therefore the spectral amplitude can be taken as: Here α is the coherent state amplitude, ϕ the phase, and ω d the center frequency.The corresponding wavepacket amplitude α(t) is obtained via the Fourier transform of α(ω) and is thus a propagating plane wave, For an ideal stationary beam in a coherent state, the mean photon flux, f (t) = â † (t)â(t) , will be timeindependent: As is expected for a stationary beam, with a constant photon flux for all times, the mean photon number n defined in Eq. (B13) is infinite, and the spectral amplitude cannot be normalized.These facts make calculations using the quantum representation of the stationary beam problematic.

Partitioning infinite temporal modes
The infinite mean photon number is problematic for calculations as it makes the photons in the external drive field ill-defined.A solution to this problem is to define a complete set of discrete, orthonormal basis functions {Φ i (t)}, which partition the continuous-mode coherent laser beam into an infinite tensor product state of discrete-mode coherent states [93].
If the basis functions {Φ i (t)} satisfy the orthogonality and completeness relations they form a non-continuous basis set with which a discrete set of annihilation operators may be created according to ĉi = dtΦ * i (t)â(t).(B19) Equation (B18) gives the inverse relation Naturally, an eigenstate of â(t) with eigenvalue α(t) is also an eigenstate of ĉi with eigenvalue It follows that a continuous-mode coherent state can be equivalently expressed as an infinite tensor product of discrete-mode coherent states: The result in Eq. (B22) is an important property for mode matching the continuous-mode coherent state to a discrete-mode.The freedom in choosing the set of basis functions {Φ i } is large.This facilitates mode matching of the discrete-mode coherent states |{α i } with a large variety of mode functions.One of the simplest examples is the partitioning into rectangular time-bins with duration T , which are described by the set of functions {Ψ m (t)} defined as [94] Ψ Above, the label z 0 denotes an arbitrarily chosen reference point along the propagation axis.The set of functions in Eq. (B23) can be extended to form a complete set that satisfies Eqs.(B17) and (B18).The corresponding eigenvalue for each of these discrete-mode coherent states can be obtained from Eq. (B21) as The duration T can be chosen arbitrarily as long as it is much larger than 1/ω d .Thus, the continuous wave laser described by the travelling plane wave α(t) can be expressed as a sequence of M → ∞ copies of the discretemode coherent states |α 0 defined by the functions in Eq. (B23).The benefit of going through all the trouble of reaching this representation is that we now have a well-defined wavepacket amplitude α 0 for each partitioned piece of the laser beam.

Appendix C: Phenomenological master equation
In this appendix, we write down a phenomenological master equation for the probabilities P n (t) of occupying the nth Fock state in the cavity.This master equation provides an analytical approach to gaining an intuition for the (N + 1)-photon processes associated with the unconventional saturation effect.The solution to this master equation was used to obtain the analytical results presented in Fig. 4(b) in the main text.
The setting for the phenomenological master equation is the effective-drive picture described in Section III C. In summary, the cavity population becomes suppressed in the weak-excitation regime due to the coupling to the emitter ensemble.This suppression can be understood by an effective drive, where the external drive and the effective driving from the ensemble interfere destructively and thus cancel the cavity population.However, the ensemble can only cancel photon numbers up to order N .Hence, the destructive interference breaks down at order N +1.The higher-order photon states (n > N ) in the external drive can be seen as multiphoton pulses driving the cavity.The effects of such pulses on the cavity dynamics can be studied with a simple phenomenological master equation, which only considers direct cavity absorption and exponential decay.Even though this model neglects all effects from the coupling to the emitters beyond the cancellation of population, it qualitatively captures the unconventional saturation effect, as shown in Fig. 4(b).
Consider a cavity in the unconventional saturation regime, which is approximately in the ground state because of the destructive interference.Additionally, consider a single pulse with k photons, which is directly absorbed into the cavity due to the intermittent saturation of the destructive interference.Neglecting further effects from the coupling to the emitters, the k photons will subsequently leak out of the cavity through single-photon processes with the rate γ c .In this scenario, the dynamics for the probabilities P n to occupy the nth Fock state can be described by the master equation on vector form ∂t P (t) = Γ P (t), (C1) where P = (P 0 , P 1 , ..., P k ) T is a column vector with the Fock states {|n , n ∈ [0, k]} and is a matrix describing the in-and out-flow of probability to occupy each state.The general solution to Eq. (C1) is P (t) = e Γt P (0), (C3) which for P (0) = P k = 1 gives Equation (C4) describes the dynamics of the exponential decay of k photons from the kth Fock state to the ground state.

Prediction of the total cavity population in the steady-state
The time-averaged contribution from a stream of independent k-photon pulses to the total cavity population in the steady state can be calculated according to The first term is the steady-state cavity population in the weak-drive regime, given by Eq. ( 24).The probability P T k of having a k-photon pulse during the time T is given by the Poisson distribution for the external drive with amplitude The factor T ss /T gives the number of time bins with duration T in the time T ss , over which we average.Taken together, the term (T ss /T )P T k gives the fraction of kphoton pulses that contributes to the time average.Using the explicit forms for the time-dependent probabilities P n (t) given in Eq. (C4), the sum in Eq. (C5) can be evaluated as k n=0 nP n (t) = ke −γct .Additionally, we can assume that the time it takes to reach the steady state is much longer than the time scale for the decay dynamics, T ss 1/γ c .Thus, we can take the limit T ss → ∞, which yields the expression When N emitters couple to the cavity, the destructive interference between the emitters and the external drive breaks down at order N + 1.Hence, the lowest-order photon pulses contributing to the cavity population in the steady state will be N + 1. Taking k = N + 1 in Eq. (C7) gives the expression for the cavity population given in Eq. ( 23) in the main text.
Equation (C7) can easily be generalized to include the contribution from independent k-photon pulses up to or- 2. Prediction of the cavity populations in the steady state Using the analytical solutions in Eq. (C4), we can also find an expression for the probability of occupying the nth Fock state in the steady state.Including only the contribution from (N + 1)-photon pulses, the expression is Performing the integral with P n (t) given by the solution to the master equation given in Eq. (C4), and taking the limit T ss → ∞, gives The probability P T N +1 is given by Eq. (C6) with k = N +1 and P n weak is given by the Poisson distribution for a coherent population with amplitude |α weak | 2 = nc weak , The second term in Eq. (C10) is proportional to Ω 2(N +1) d .Thus, this term explains the cross-over from the coupled distribution to the uncoupled distribution [P αc (n) and P α o c (n) in Fig. 6, respectively] at n = N + 1.In the weak-drive regime, before the system has entered the unconventional saturation regime, the first term P n weak will be much larger than the contribution from the (N + 1)-photon pulses for n < N + 1.Therefore, the P n ss will still have a Poisson-like distribution, but will deviate from the coupled distribution because of the contribution from the (N +1)-photon pulses given by the second term in Eq. (C10).On the other hand, when the system is in the unconventional saturation regime, the probability P T N +1 of having (N + 1) photons in the drive has become large enough for the second term in Eq. (C10) to dominate.Hence giving P 1 ss ≈ P 2 ss ... ≈ P N +1 ss .Thus, Eq. (C10) explains the observation of the relationship (N + 1)ρ N +1 ≈ N ρ N ≈ ... ≈ ρ 1 between the populations in the unconventional saturation regime.In Fig. 9, we compare the P n ss given in Eq. (C10) (right panel) with the cavity populations ρ n ≡ n|ρ c |n from Fig. 6[(b),(c),(d)] (left panel).This comparison shows that the simple phenomenological master equation described in this appendix gives a good analytical intuition to the cavity response to external driving observed in the main text.The analytical results also confirm the expectations on the populations ρ n discussed in Section III F. regime cannot be neglected.Instead, we have performed a similar derivation as in Ref. [96,97] to obtain the following master equation: Here, the counter-rotating terms â † σ+i and âσ −i are included in the Hamiltonian, and The positive-frequency operators X+ c and X+ ei are defined as where |j and |k are eigenstates to the undriven Tavis-Cummings Hamiltonian before making the RWA.This master equation brings the system to the correct dressed ground state, including counter-rotating terms, and is suitable for systems with overlapping transitions.The derivation of a collective dissipator, including counter-rotating terms, is analogous to the derivation of Eq. (D4).The result can be directly written down by replacing σ+i and σ−i with the collective pseudo-spin operators S + ≡ N i=1 σ+i and S − ≡ N i=1 σ−i .The corresponding collective dissipation rate is γ col = γ e /N .
In Fig. 10, we compare the cavity response to external driving obtained with the different master equations presented above.Solid, blue, and grey curves show the cavity response to external driving for N = 1 − 4 quantum emitters obtained with Eq. (D1), which is suitable for the Tavis-Cummings model and individual emitter decay as presented in the main text.Red diamonds show the cavity response calculated with counter-rotating terms in the cavity-emitter interaction and individual decay of the emitters.Lastly, dotted red curves show the cavity response calculated with counter-rotating terms and collective decay of the emitters.As can be seen, the presented changes to the master equation do not affect the unconventional saturation effect qualitatively.Thus, the RWA captures all the essential physics we demonstrate in our manuscript.
Note that the quantity plotted on the y-axis is X− X+ ss and not directly â † â ss .The operators X+ = j,k>j j|â † + â|k |j k| and X− = ( X+ ) † correspond to the positive and negative frequency component of the photon-like field inside the cavity, respectively.With a correct treatment of input-output theory, the field emitted from the cavity is proportional to X+ [69,98].Thus, X+ is the relevant field operator for experiments.When the RWA is applied, X+ = â when the cavity and emitters are on resonance.On the other hand, when the RWA is not applied, X+ = X+ c , with X+ c given in Eq. (D3).Thus, the notation X− X+ ss allows us to compare the results obtained with the different master equations.Following this note, Fig. 10 further tells us that the unconventional saturation effect appears in the excitation of the photon-like field (counted with the operator X− X+ ) that leaks out of the cavity and not directly in the population of the cavity Fock states (counted with the operator â † â).Therefore, there is no distinction between applying or not applying the RWA for experiments.
To obtain the cavity response with counter-rotating terms in the Hamiltonian (diamonds and dotted curves in Fig. 10), we have slightly changed the drive frequency ω d .Including the counter-rotating terms in the Hamiltonian, we find a minor shift of the transparency dip observed in Fig. 2(b).This shift slightly changes the drive frequency at which the effect is most pronounced.With the current set of parameters, we find that this shift is less than 0.5 % for N = 1 and less than 0.3 % for N = 4. Since this shift is so small, not correcting the drive frequency will still show the effect and we do not find any qualitative changes in the unconventional saturation effect.

Appendix E: Cooperativity
The cooperativity is defined as In Fig. 11, the panels with break-point predictions from Fig. 5 in the main text are shown together with the corresponding cooperativities.As can be seen, the observed saturation effect and our analytical expression for the critical drive are both robust to a wide range of cooperativities, 25 ≤ C ≤ 1600.It is also evident from Fig. 11, that a large cooperativity facilitates the observation of the saturation effect as it pushes the emergence of the nonlinear effect to lower drive strengths.
To find out how robust the observed effect is at low cooperativity, we performed a few simulations with one and two emitters in the cavity for different C. The results from these investigations are presented below and are structured as follows.First, a set of plots with a semilossy cavity, i.e., γ c ≈ 0.03 ω c , is presented in Fig. 12.This cavity loss rate is the same as was used to produce the plots in Fig. 11.After that, a second set of plots, where γ c ≈ 0.17 ω c , is presented in Fig. 13.This loss rate corresponds, e.g., to a localised surface plasmon mode in a metal nanoparticle.Moreover, both Fig. 12 and 13 have two columns.The first column shows the response with low-loss emitters, γ e = 0.01 γ c , and the second column shows intermediate-loss emitters with γ e = 0.1 γ c .According to Eq. (E1), a more lossy emitter ensemble can be compensated to have the same cooperativity as a less lossy one (γ c fixed) by increasing the cavity-emitter interaction strength.Therefore, the two columns could also be regarded as corresponding to very weak coupling in the left column, and weak coupling in the right column.

FIG. 2 .
FIG. 2. Master-equation calculations of the cavity populations in the steady state for coupled cavity-N -emitter systems with the same collective interaction strength g col .The spectra are compared with the analogue, classical calculation.(a) The spectra for N = 1 − 4 quantum emitters show only minor differences between different N and the classical CO model when plotted on a linear scale.(b) The spectra viewed on a logarithmic scale, on the contrary, show considerable differences of several orders of magnitude for resonant driving.

N 1
σ+i σ−i ss as a function

FIG. 3 .
FIG. 3. Log-log plots of the average steady-state populations as a function of the normalized drive strength Ω d /g col , calculated with the master equation for N = 1 − 4. (a) The cavity population shows an N -dependent transition through a nonlinear response regime for intermediate drive strengths between two linear asymptotes nc (dashed black line) and n o c (dotted black line).The dashed vertical line marks Ω d /g col = 0.25, which was used for calculating the spectra in Fig. 2. (b) The total ensemble population for the same range of drive strengths as in panel (a).This population saturates at N/2 when the drive is strong.

FIG. 5 .
FIG.5.Calculations of the steady-state cavity population for different [(a),(b),(c),(d)] emitter decay rates γe and [(e),(f),(g),(h)] collective interaction strengths g col .The remaining parameters were held fixed, using the same values as in Fig.3.The red stars mark the derived analytical expression for the critical drive Ωcr.As can be seen, the analytical results predict exceedingly well the onset of the nonlinear regime.

FIG. 6 .
FIG. 6. Visualization of the evolution of [(b),(c),(d)] the multiphoton populations ρn as the drive strength Ω d is increased; the total population is shown in panel (a).The dashed and solid lines mark the Poisson distributions P (n) calculated with the coherent-state amplitudes corresponding to a coupled, |αc| 2 = Ω 2 eff /γ 2 c , and uncoupled, |α oc | 2 = Ω 2 d /γ 2 c, cavity, respectively.As can be seen for all system sizes N = 1-4, the emitter ensemble behaves as a true coherent state in the coupled system only up to first order.Then, for increasing drive strength, each system saturates in turn, visible as the higher-order photon states approach the uncoupled-cavity distribution.The red circles illustrate when ρN+1 have become comparable to ρ1, i.e., ρ1 ∼ (N + 1)ρN+1, which is a signature of the unconventional saturation regime.

FIG. 7 .
FIG.7.Spectra of the average cavity population â † â for N = 1, 2, and 3 emitters plotted against the classical analogue nc on a [(a),(b),(c)] linear and [(d),(e),(f)] logarithmic scale.As is clearly seen throughout the log-log plots, the quantum interference effect is present for a wide range of interaction strengths below and at the strong-coupling limit g col ≥ γc/2.

FIG. 8 .
FIG. 8. Comparison of the steady-state cavity population calculated withEq.(C7), including the contribution of multiphoton pulses up to order M = N + 5, and the cavity population obtained with the Lindblad master equation from the main text.
for N = 1 − 4, calculated with Eq. (C8) and M = N +5, are plotted as dashed red curves in Fig. 8. Solid blue curves show the results calculated with the Lindblad master equation from the main text.A comparison of Fig. 8 with Fig. 4(b) in the main text shows that including a few more orders of multiphoton absorption in the simple model described in this section qualitatively captures the unconventional saturation effect for drive strengths close to the traditional saturation point Ω d ∼ g col .The extension does not, however, change the behaviour in the intermediate drive regime (Ω d > g col ).

FIG. 10 .
FIG.10.Visualization of the effects of counter-rotating terms (diamonds and red dotted curves) and collective emitter dissipation (red dotted curves) on the unconventional saturation effect.Solid curves reproduce the results of Fig.3(a) from the main text.

FIG. 11 .
FIG.11.The critical drive panels from the main text with the corresponding cooperativities C. In panel (a -d), g col and γc have been held fixed while γe has been varied from 0.5 -5% of γc going left to right.In panel (e -h), γc and γe have been held fixed while g col has been varied from 2 -0.25 γc going left to right.

FIG. 12 .
FIG.12.Steady-state cavity populations as a function of drive strength for cavity loss rate γc ≈ 0.03 ωc and varying emitter loss rate and collective coupling.

FIG. 13 .
FIG.13.Steady-state cavity populations as a function of drive strength for cavity loss rate γc ≈ 0.17 ωc and varying emitter loss rate and collective coupling.