Identifying unbound strong bunching and the breakdown of the Rotating Wave Approximation in the quantum Rabi model

We use a recently derived gauge-invariant formulation of the problem of a two-level system coupled to an optical cavity, to explore the transition between the weak, and the ultra-strong coupling regimes of light-matter interaction. We explore this transition using the intensity correlations $g^{(2)}(\tau)$ of the emitted light, and find strong, unbounded bunching of the emission from systems governed by the Rabi Hamiltonian. Surprisingly, this effect is observed not only in the ultra-strong coupling regime, but also for weakly coupled systems, where the Jaynes-Cummings Hamiltonian would predict the opposite, antibunched emission. This suggests that the higher-order correlations are a particularly sensitive probe of the divergence between the Jaynes-Cummings and Rabi Hamiltonians, and can serve as an indicator of the breakdown of the rotating wave approximation. Our findings indicate also that the boundary between the weakly, strongly, and ultra-strongly coupled dynamics, is much richer than currently accepted.


I. INTRODUCTION
Many systems studied in cavity quantum optics are variations of its fundamental workhorse: a two-level system (TLS) coupled to an optical cavity [1,2].The behavior of the cavity-TLS system (CTS) is dictated by the relationships between the coupling strength (g) and the characteristic resonant frequencies (ω σ and ω a ) and dissipation rates (γ and κ) of the TLS and the cavity, respectively.In particular, if g is smaller than the dominant losses of the system (g ≲ (κ + γ)/2, defining the weak coupling (WC) regime), any energy that enters the system is likely lost before the exchange of excitations between the cavity and the TLS can occur.On the other hand, in the strong coupling (SC) regime for which g exceeds the losses of the system (g ≳ (κ + γ)/2), the cavity and the TLS can coherently exchange excitations before the energy is dissipated, inducing a hybridization in the response of the CTS [3][4][5].
Large coupling strengths can lead the system into the so-called ultra-strong coupling (USC) regime, where g becomes comparable to the resonant frequencies of the cavity and TLS [6].The phenomenological limit for the onset of the USC is typically defined as g ≳ 0.1ω a , or g ≳ 0.1ω σ (USC systems are commonly studied in a resonant configuration where ω a = ω σ ).Importantly, systems operating in the USC regime exhibit new characteristics and non-trivial properties, such as the existence of a non-vacuum ground state [7,8], which would be largely missed if attempting to describe them using the typically approximations applied in the WC or SC regimes -most importantly, the Rotating Wave Approximation (RWA).
The emergence of new effects in the USC, and its boundary with WC and SC regimes, is typically quantified through the study of one-photon spectra S (1) (ω) of the emitted light, by tracing the intensities and frequencies of the spectral features, leading to the usual criterion for USC (g ≳ 0.1ω a ).Alternatively, we could choose to focus on the statistics, rather than spectrum, of the emitted light, by measuring its intensity correlations g (2) (τ ), through the Hanbury Brown and Twiss (HBT) interferometer [9][10][11][12] shown schematically in Fig. 1(a).Intensity correlations should constitute a more sensitive probe to the non-number-conserving dynamics generated by the Rabi Hamiltonian [13], and might delineate an far more complex boundary of the USC regime.
In this work, we analyze the emission of an incoherently-pumped CTS, operating in the WC, SC, and the USC regimes.We focus on studying the intensity correlations g (2) (τ = 0) obtained with two models: the Quantum Rabi model (QRM), and its widely used approximation -the Jaynes-Cummings model (JCM).We embrace the formulation of the QRM recently derived in a series of papers which have reconciled long-standing questions about ensuring gauge invariance [14][15][16][17][18][19], and proposed a complete description of the interaction between an USC system and the environment.This formulation of the gauge-invariant QRM offers an opportunity to carefully study two questions, which are of critical importance to the field of cavity electrodynamics: (i) how does the statistics of emission from a CTS change as we transition between different coupling regimes?(ii) how does the JCM break down in the USC regime?To illustrate these effects, in Fig. 1(b) we plot g (2) (0) calculated using the QRM (solid lines) and the JCM (dashed lines), and two different pumping rates (details of the models and excitation schemes are discussed in Section II).In weakly pumped strongly pumped Setup of the studied system.An incoherently driven TLS interacts with an optical cavity.The emission of the cavity is analyzed with a Hanbury Brown and Twiss (HBT) interferometer (further discussed in the text), giving the intensity correlations g (2) (τ ).(b) Intensity correlations g (2) (0) as a function of the coupling parameter η = g/ω0, calculated using the QRM (solid lines; see details in Section II A) and the JCM (dashed lines; Section II B).In the limit of strong incoherent pumping rate (Γ/γ = 10, orange lines), the results of the two models coincide in the WC regime (η ≲ 2.5 × 10 −2 ), and become very different for USC η ≳ 0.1; for the weak pumping rate (Γ/γ = 10 −6 , dark blue lines), correlations obtained with the JCM and the QRM differ in the USC and WC regimes, for any η ≳ 5 × 10 −3 .
contrast to the reported works on the emission from a thermally pumped CTS [20,21], we find that under incoherent illumination of the TLS, the emission of the CTS appears to exhibit a seemingly unbounded bunching g (2) (0) ≫ 1 [22].Furthermore, we find that the two models (JCM and QRM) can deviate significantly both in the USC (see the dashed and solid orange lines diverging in Fig. 1(b) for η = g/ω 0 ≳ 0.1), as well as in the WC regime (blue lines diverge for η ∼ 10 −2 ; for the parameters used in Fig. 1, the WC is defined as η ≲ 2.5 × 10 −2 , see Section II A for details), depending on the incoherent pumping rate.
The strong deviation between the JC and QR models, observed even in the WC regime in Fig. 1(b), is a surprising result with far-reaching consequences.The JCM is the default model for the majority of quantum-optical systems operating in the WC and SC regimes, including the circuit QED [23], exciton polaritons [24], or quantum plasmonics [25].It is therefore critical to understand how the Jaynes-Cummings model unravels when probed using higher-order correlations, by comparing it to more complete Quantum Rabi model, and carefully analyzing their relationship.
This work is structured as follows: in Section II we formally introduce the system under study, the Quantum Rabi and Jaynes-Cumming Hamiltonians, and describe the formulation of the excitation, emission, and dissipation of the system.In Section III we identify the key mechanism which gives rise to this strong bunching.Finally, in Section IV we probe the extent of this effect in the USC, SC, and WC regimes, and discuss how the intensity correlations can help us identify the breakdown of the RWA.

II. FRAMEWORK OF THE MODELS
In this Section we briefly review the formulation of the Quantum Rabi (Subsection II A) and Jaynes-Cummings (Subsection II B) Hamiltonians.We describe how to formally address the interaction of such systems with the environment, and how to access the spectra, and the intensity correlations of the emitted light.

A. Quantum Rabi model
The exact description of the interaction between the cavity and the TLS is given by the Quantum Rabi Hamiltonian, which we choose to write in the Coulomb gauge [14, 15][26] as Here ω a and ω σ are the resonant frequencies of the cavity and TLS, respectively.In this work, we consider a resonant system with ω 0 ≡ ω a = ω σ , and introduce the normalized coupling parameter η = g/ω 0 .Operator â denotes the cavity photon annihilation and σz = σσ † − σ † σ and σy = i(σ † − σ) are the z and y Pauli operators associated with the TLS.The Ĥ QR Hamiltonian does not conserve the number of excitations, but conserves its parity.
To characterize the statistical properties of the emission of the CTS, we study the intensity correlations of the emitted photons, g (2) (τ = 0), measured as where I 1 (t + τ ) and I 2 (t) are the photocurrents registered by the two detectors of the HBT interferometer, and τ is the time delay between the detection events (see Fig. 1(a)).For zero time delay (τ = 0) and sufficiently large t (so the system reaches the steady state), this quantity is related to the statistics of photons inside the cavity as [12]: where ⟨ Ô⟩ ss denotes the expectation value of operator Ô in the steady state (ss).The operator is the dressed operator of the cavity, which mediates the losses, and emission from the cavity [15,17].It is expressed in the basis of eigenstates |µ⟩ R because the emission from the cavity occurs due to the transitions between the eigenstates of the Hamiltonian (see Refs. 16 and 27 for an in-depth discussion of this formulation and the role of the secular approximation).Kets |ν⟩ R and |µ⟩ R in Eq. ( 4) are the eigenvectors of the QRM Hamiltonian, and ℏω ν > ℏω µ are their respective eigenvalues, plotted as a function of η in Fig. 2(b).The notation for QRM eigenstates used throughout the text, |ν⟩ R ≡ |n±⟩ R , is chosen to recall the JCM polaritons (|n±⟩ = (|n, g⟩ ± |n − 1, e⟩)/ √ 2), as the two match in the limit of vanishing coupling g (see Fig. 2(a) and (b)).Note that throughout the text we keep the labeling of the eigenstates even after the eigenvalues crossing points (see color code in the figure).
The dynamics of the system, described through the density operator ρ, is given by the master equation with the gauge-invariant formulation of the Hamiltonian [14,17] and dissipation terms [15,16,27]: with L denoting the Lindbladian superoperator, and the Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) terms expressed through the dressed operators xa (Eq.4) and xσ [15]: The GKSL terms are introduced to describe the interaction between the TLS, cavity, and the bath following [15,19,28], including modelling the incoherent pumping of the TLS (with rate Γ), the dissipation of the TLS (γ) and the cavity (κ).The dissipation is assumed additive, and we neglect any bath-mediated TLS-cavity interaction [7,16].We define the steady state of the system via the density operator ρss as ∂ t ρss = 0. We choose the pumping mechanism and parameters which correspond to recent experiments with plasmonic CTS systems, which have reportedly reached values of η close to 0.1 [29][30][31][32][33][34][35][36][37][38].In particular, we explore the emission from a CTS under incoherent pumping of the TLS, corresponding to driving the atomic systems with a far blue-detuned laser, followed by a spontaneous cascade to the excited state of the TLS; see Fig. 1(a).In the plasmonic CTS considered here, the cavity is characterized by resonant frequency ω 0 (set to ℏω 0 = 1 eV), and a low quality factor Q = ω 0 /κ, defined by the dissipation rate κ, and set to Q = 20.The emitter (such as a single molecule [29,31], or a quantum dot [34,37,38]), is modelled as a TLS with decay rate γ/ω 0 = 10 −3 .These parameters establish the upper limit for the WC regime η = g/ω 0 < κ/(2ω 0 ) = 0.025.For reference, we note that the characteristic cooperativities C = 4g 2 /(κγ) in our study are C = 0.08 (for η = 10 −3 ) and C = 800 (for η = 0.1).

B. Jaynes-Cummings model
The Jaynes-Cummings model (JCM) can be derived from the QRM by taking three approximations.First, we expand the interaction term in the QRM Hamiltonian as σz cos 2η(â + â † ) + σy sin 2η(â + â † ) and drop the terms nonlinear in η.Next, we introduce the rotating wave approximation (RWA) by removing the so-called non-number-conserving terms σâ+ σ † â † , to find the JCM Hamiltonian Figure 2(a) shows the eigenvalues of the JCM Hamiltonian which, for small η, converge with the eigenvalues of the QRM Hamiltonian (Fig. 2(b)).
The third approximation in the JCM regards the emission, dissipation, and absorption of the system, which, within the JCM, are mediated by the bare â and σ operators (instead of the dressed xa and xσ operators).The correlations arising in this model can thus be calculated as where the steady-state ρ(JC) ss is obtained from the standard master equation [42], as the solution to ∂ t ρ(JC) ss = 0.The JCM summarized in Eqs. ( 8) and (10) has been extensively used to describe the properties of weaklycoupled CTSs [12].In this work, it establishes a point of comparison to identify the features that can arise in the QRM.We use the JCM to calculate the intensity correlations g (2) JC (0), shown in dashed orange and blue lines in Fig. 1(b) for the strong and weak incoherent pumping, respectively.In the former case, the JCM correctly reproduces the results of the exact QRM below the USC threshold η ≲ 0.1, but fails for larger η, where the g (2) (0) obtained with the JCM saturates below g (2) (0) < 1 (see Appendix A for analytical equations), while, in the exact QRM, g (2) (0) keeps increasing with η.Furthermore, for the weak incoherent pumping (blue lines), we find significant differences between the JCM and QRM even in the WC regime (the JCM and QRM differ for η ≳ 5 × 10 −3 ).This unexpected breakdown of the JCM for small η is discussed in detail below.

III. ORIGIN OF THE BUNCHING IN ULTRASTRONGLY COUPLED SYSTEMS
In this Section, we demonstrate that the strong bunching identified in the USC regime in Fig. 1 exact values of g (2) (0) (solid blue line); the diagonal approximation truncated to the states 11) (solid orange line); the diagonal approximation considering only the to characteristics (decay pathways and population) of the single |3−⟩ R eigenstate.To justify the focus on that particular eigenstate of the QRM Hamiltonian, in Fig. 3 we plot the exact values of g (2) (0) (solid blue line) obtained for an intermediate pumping Γ/γ = 10 −3 together with approximated results.We start by approximating the steady-state density matrix as being diagonal in the basis of the |ν⟩ R eigenstates of the QRM Hamiltonian, The approximation is justified by noting that the system is driven with an incoherent pumping mechanism, andmuch like any weakly coupled, thermally populated system -should exhibit limited coherence.We explore the validity of this approximation in Appendix B, where we show the values of the steady state density matrix as a function of the coupling strength.We can thus approximate g (2) (0) (Eq.( 3)) as (see derivation in Appendix E) where |µ⟩ R and |ν⟩ R are the eigenstates of the QRM Hamiltonian.The intensity correlations calculated by truncating the double sum in the numerator up to |3−⟩ R are shown with the solid orange line in Fig. 3 -this approximation gives a very good agreement with the exact calculation for η ≳ 2.5 × 10 −2 ; as we have numerically verified, the deviation observed in the WC regime η ≲ 2.5 × 10 −2 is not due to the truncation of the basis, but rather due to the effect of the off-diagonal terms of ρss .
We can further approximate g (2) (0) by by limiting the double sum in Eq. ( 12) over ν and µ to the |ν⟩ = |3−⟩ R and |µ⟩ = |1−⟩ R term: This approximation, denoted in Fig. 3 with the green line, explores the role of the particular, correlated two-photon emission from |3−⟩ R to the |1−⟩ R state.Equation (13) shows that the singular role of the |3−⟩ R polariton in the onset of bunching can be ascribed to two effects: the presence of a two-photon emission pathway towards |1−⟩ R , and a comparatively large, non-thermal population of this state, expressed in our model as R 3− .
To analyze the latter effect in detail, we plot in Fig. 4 the populations R ν of all the relevant eigenstates of the QRM Hamiltonian as a function of the normalized coupling η for the intermediate pumping Γ/γ = 10 −3 .Initially, the eigenstates are populated according to their respective eigenvalues, with the lower-energy states being more populated.While this mechanism holds for η < 0.05, Fig. 4 shows that for larger coupling the population R 3− (orange solid line) rapidly grows at a rate proportional to η 4 , becoming larger than R 2+ (purple solid line), and eventually approaching R 1+ (red solid line) for η > 0.1.This increase of R 3− far exceeds that observed within the JCM (see the orange dashed line denoting the population of |3−⟩).
This large population of the |3−⟩ R polariton can be attributed to the new, direct excitation pathway introduced in the QRM (see Fig. 5(a)).In this model, the |3−⟩ R polariton can be directly driven from the ground state |0⟩ R , through the Γ/2D x † σ incoherent pumping term introduced in the master equation (5).We illustrate this effect in Fig. 5 The direct excitation pathway should be proportional to the incoherent pumping rate Γ, and compete with the sequential excitation pathway |0, g⟩ R → |1±⟩ R → |2±⟩ R → |3−⟩ R , with total rate ∝ Γ 3 .We can identify this competition in Fig. 5(c), by plotting the population R 3− as a function of Γ, for a range of coupling parameters η.The direct mechanism dominates the pumping for small Γ, where its linear dependence on the pumping rate makes it more efficient than the sequential pumping mechanism.Only when we increase Γ, does the latter process become more efficient and we recover R 3− ∝ Γ 3 [43].The transition from the linear to the cubic dependence on Γ, or from the direct to sequential excitation mechanisms, shifts towards larger Γ as we increase the coupling strength η.This is because the overall efficiency of the direct excitation, governed by the matrix element  (2) (0) < 1) and bunching (g (2) (0) > 1) regions, respectively.The color scale of g (2) (0) is linear from g (2) (0) = 0 to 2 and logarithmic from g (2) (0) = 2 to 10, where it saturates.
In Section III we traced the bunching in the USC regime seen in Fig. 1(b) and Fig. 3 to the direct excitation, and two-photon emission from the |3−⟩ R polariton.Here, we explore the extent of this new effect, and identify the threshold for the deviation between the JCM and the QRM.Figures 6(a) and (b) compare the values of g (2) (0) calculated with the QRM and the JCM, respec-tively, for a range of the normalized coupling strengths η and pumping rates Γ, and demonstrate that the deviation between JCM and QRM also depends strongly on the latter parameter (Γ), in a manner that can be explained using the formulation laid out in Section III.

A. Qualitative dependence of the bunching on the pumping
The map of intensity correlation shown in Fig. 6(a) indicates that the bunching observed in the QRM depends on the rate of incoherent pumping Γ for a wide range of coupling parameters η.We explore this effect in more detail in Fig. 7(a), where we plot vertical cross sections of Fig. 6(a) -the dependence of g (2) (0) on Γ -, for η = 0.1 to 1, and find that the intensity correlations follow This dependence results from the direct excitation mechanism of the polariton |3−⟩ R : In the USC regime, the contribution from that polariton dominates the numerator of g (2) (0) approximated as in Eq. ( 13), with the numerator proportional to R 3− ∝ Γ (as discussed in section III).Conversely, the denominator of Eq. ( 13) is dominated by the contribution from the |1±⟩ R eigenstates ∝ R 1± .Since both |1−⟩ R and |1+⟩ R are populated directly from the ground state, we find R 1± ∝ Γ.Thus we recover the g (2) (0) ∝ R 3− /(R 1± ) 2 ∝ 1/Γ dependence, and g (2) (0) is therefore unbounded.Equation ( 14) can be used as a marker for the breakdown of the JCM, since in that model, the direct pumping mechanism of high-order polaritons is absent, and consequently the intensity correlations are largely constant for small Γ (see Fig. 6(b)).This is further explored in the following subsection.

B. Breakdown of the JCM Hamiltonian in the WC regime
Figure 6(a) shows that in the QRM, the strong bunching region can also be identified for η in the WC and SC regimes.In particular, for very weak incoherent pumping (Γ/γ ≲ 10 −5 ), the strong bunching appears for couplings as small as η ≈ 2.5×10 −2 , highlighting a significant deviation between the predictions of the JCM (where g (2) < 1 for all η) and of the QRM.To demonstrate this effect more clearly, in Fig. 7(b) we plot horizontal cross sections of Fig. 6 -the dependence of g (2) (0) on η calculated for different pumping rates Γ using the QRM (solid lines) and in the JCM (dashed lines) models.
The difference between the QRM and JCM in the WC and SC regimes has the same origin as in the USC regime: the direct excitation of the |3−⟩ R eigenstate in the QRM.As we discuss in Section III, this direct mechanism can be dominant for any η, providing that the pumping is sufficiently small.Conversely, in the JCM, for small Γ (Γ ≪ γ), g (2) (0) is determined by the population of the polariton with two-photon terms |2±⟩ (we neglect the far smaller contribution from |3±⟩).Since this polariton is excited through a sequential process, its population is proportional to Γ 2 .Normalized by the square of the population of the |1±⟩ states (∝ Γ 2 ), the intensity correlation in the weak pumping and coupling limit of the JCM is independent of Γ, and does not support bunched emission pathways (we derive this result explicitly in Appendix A).
It is worthwhile to consider how these striking differences can emerge in the WC limit, and far below the USC threshold g ≪ 0.1ω 0 , where the QRM and JCM would be conventionally expected to match.To derive ĤJC , in Eq. ( 8) we (i) performed the expansion of the QRM Hamiltonian Ĥ in a power series of η, and dropped terms scaling with higher powers of η, and (ii) applied the RWA to remove the non-number-conserving terms.Figure 6 shows that QRM and JCM can diverge for arbitrarily small η, suggesting that the error in g (2) (0) introduced by the series truncation in (i) can be made arbitrarily small.Thus, we can trace the observed differences to the application of the RWA.This is, to our knowledge, the first proposal for observing the breakdown of the RWA in the WC limit.An experimental assessment of this effect would constitute a challenge, mostly due to the very short decoherence time of the molecules, further reduced by the highcooperativity coupling to the cavities, resulting in sub-ps timescales for the detection of bunching (see discussion in Appendix G).These requirements can be partially relaxed by considering atoms with smaller γ, and reduced cooperativity.Alternatively, a weaker pumping -while reducing the emission and coincidence rate -should boost the bunching effect (see Fig. 6), making its observation possible even after averaging over the longer detection window.

C. Emission spectra
To conclude our characterization of the QRM, we now study the dependence of the one-photon emission spectrum S (1) (ω) on the coupling strength.We show that the differences below the USC between the QRM and the JCM are much smaller for S (1) (ω) than for the intensity correlations g (2) (0), emphasizing the interest of using the intensity correlation to identify the breakdown of the JCM.
In Fig. 8, we plot the emission spectra S (1) (ω) for a range of coupling strengths η.We calculate S (1) (ω) within the QRM as [15,17,42] where xa and x † a are the dressed annihilation and creation operators of the cavity, respectively (see Eq. ( 4)).For a detailed study of the emission from the TLS (S (1) σ calculated from the two-time correlators between xσ operators) in the USC regime, we direct the reader to Ref. [19].Each spectral feature in the figure corresponds to a transition between the eigenstates of the system, which we identify by matching the frequency of the peak with the difference between the eigenfrequencies of the system (see Fig. 2(b)).We plot these spectra as functions of the rescaled frequency (ω − ω 0 )/η to compensate for the ηdependence of the Rabi splitting between the polaritons.
The smallest coupling strength considered in Fig. 8 corresponds to the low limit of the USC, η = 0.1.This spectrum displays two peaks of similar intensity which, as in the JCM, correspond to the emission of a single photon via the |1−⟩ R → |0⟩ R (lower frequency feature) and |1+⟩ R → |0⟩ R (higher frequency) transitions.Only when we increase the coupling strength to about η ≈ 0.3, do the spectra develop additional features: the strength of the spectral peak corresponding to the The visibility of each peak can be compared to the populations of the initial states participating in the emission processes.For instance, the two new peaks follow the same dependence on η as the populations R 3− and R 2− , respectively (Fig. 4), so that, for η > 0.05, the two populations begin to grow very rapidly with increasing η.Simultaneously, as we increase the coupling η, spectral features continuously shift, reflecting the changes to the spectrum of the QRM Hamiltonian (Fig. 2(b)).These additional peaks emerge only in the USC regime, for η ≳ 0.3.To check whether this threshold can be significantly modified, as we observed for the intensity correlations, in Fig. 9 we compare the spectra obtained from the JCM (dashed lines) and the QRM (solid lines) using weaker incoherent pumping Γ, and various couplings η = 10 −3 , 10 −2 , and 0.1 (note that η = 10 −3 and 10 −2 are below the minimum value considered in Fig. 8, and give no peak-splitting).The spectra point to small differences between the two models for the larger η for both considered pumping rates, but these differences would be likely difficult to identify in experimental settings.On the other hand, as shown in Figs. 6 and  the g (2) (0) obtained for Γ/γ = 10 −3 and Γ/γ = 10 −6 show qualitative differences for η ≳ 2.5 × 10 −2 due to the direct excitation pathway introduced in Section III.The lack of sensitivity of S (1) (ω) to the direct excitation pathway is due to the fact that the emission from the lower |1−⟩ R and |1+⟩ R states predominantly govern the emission spectra for η ≲ 0.3.The excitation and emission from these eigenstates are not impacted by the direct excitation mechanism discussed above, and do not flag the breakdown of the RWA.We thus conclude that the characterization of the correlations is a more powerful tool than measuring the one-photon emission spectra for the identification of phenomena caused by the non-numberconserving terms of the QRM Hamiltonian for coupling below the traditional USC threshold η > 0.1.

V. CONCLUSIONS
In this work, we analyze the statistics of the emission from a generic quantum system comprising an incoherently driven two-level emitter interacting with a cavity.We identify an emergence of unbounded bunching as the system approaches the USC regime.By expressing the dynamics of the system in the basis of the polaritonic eigenstates of the quantum Rabi Hamiltonian, we can attribute the bunching to the singular behavior of the individual eigenstate |3−⟩ R , which (i) decays through a correlated two-photon emission, and (ii) is very strongly populated by a new, direct excitation mechanism from the ground state.
Our work shows that intensity correlations g (2) (0) are a much more sensitive tool for observing the phenomena induced by the non-number-conserving terms in the QRM, than the one-photon emission spectra.Indeed, we find that the intensity correlations can identify a breakdown of the rotating wave approximation far below the conventional limit of the USC, with the exact limit determined by the rate of incoherent pumping.These findings calls for an experimental verification, and further theoretical studies, to verify the robustness of the identified excitation and emission mechanisms.Our model can be extended to account for the more complex decay dynamics and energy structure of the quantum emitter, involving dark excitonic states, or pure dephasing, as well as the interaction with a structured reservoir. with and vector b = (0, Γ, 0, 0, 0, 0, 0, 0) T .Here we truncated the set of equations by considering that the pumping rate of the system, Γ, is very small.Thus, in the steady-state, the TLS is mostly in the ground state.Hence, we approximate [45]: (i) ⟨σσ † ⟩ ≈ 1, (ii) ⟨â † â † ââσσ † ⟩ ≈ ⟨â † â † ââ⟩, and (iii) ⟨â † â † ââσ † σ⟩ ≈ 0. In the steady state d dt ⟨v⟩ ss = 0, and we can derive closed expressions for ⟨â † â⟩ ss and ⟨â † â † ââ⟩ ss .Instead of writing their full forms, we are interested only in the values of ⟨â † â⟩ ss and ⟨â † â † ââ⟩ ss in the two opposite limits of g ≪ Γ (vanishing coupling), and g ≫ κ (beyond the SC regime).We obtain an expression that is valid in both limits by retaining the terms dependent on the leading powers of the two free parameters Γ and g: where C = 4g 2 /(κγ) is the cooperativity, and The Γ 2 dependence of ⟨â † â † ââ⟩ ss reflects our intuition that the two-photon states in the JCM should be excited via a sequential excitation from the ground state.This intuition fails for the QRM Hamiltonian, as discussed in Section IV.The intensity correlations are found, in the two limits of interest, as Note that g (2) here is derived by neglecting the direct emission from the TLS, which is an invalid approximation for very small g in equation (A5) (unless the emission of the TLS is filtered-out).If direct emission of the TLS is included (A5) should be modified.The more clear situation is in the limit of g = 0.In this case, the emission of the system is only given by the TLS, which emits one photon at a time resulting in g (2) = 0. Appendix C: Effective thermal pumping While our work focuses on the experimentally-viable mechanism of incoherent driving of the TLS, it is worth comparing this framework to other driving mechanisms.In particular, we can consider the pumping of both the TLS and the cavity due to the coupling with two thermal baths, both at the same temperature, as discussed in Refs.20 and 21.In this case, the steady state of the system is given by the statistical mixture of the eigenstates of the QRM Hamiltonian, with parameters where E ν is the energy of the νth eigenstate (calculated with respect to the energy of the ground state), k B is the Boltzmann constant, and T is the temperature of the thermal bath.
Figure 11 shows the correlations obtained with the QRM considering the pumping by a thermal bath at different temperatures (from T = 750 K to T = 3000 K, for reference we also include a secondary vertical axis in the figure with k B T /(ℏω 0 ) values).For each temperature, we show the dependence of g (2) (0) on the coupling strength η.For η ≲ 0.1, Fig. 11 shows that regardless of the coupling strength between the cavity and the TLS, the emission of the system follows a thermal statistic, g (2) (0) = 2 [20].On the other hand, for η ≳ 0.1, the system emission becomes non-classical: For low temperatures T ≲ 1750 K and large couplings η ≳ 0.1, the system emission results in antibunching with g (2) (0) < 1.For high temperatures, T ≳ 1750 K and large couplings η ≳ 0.1, the emission of the system results in a strong bunching with g (2) (0) > 2. This behavior is very different to the one studied in the main text for incoherent illumination of the TLS.The intensity correlation under thermal pumping in Fig. 11 attains a maximum value of 4 (compared to ≈ 10 7 in Fig. 7) and it does not show any signature of the unbound increase for weaker pumping.
The general behavior we show in Fig. 11 is in quali-FIG.11. g (2) (0) of a thermally pumped CTS as a function of the coupling strength η, and of the temperature T of the thermal bath (for reference we give a secondary vertical axis with the values of kBT /(ℏω0)).The results of the calculation are obtained within the QRM.The colormap indicates in blue and red the antibunching (g (2) (0) < 1) and bunching (g (2) (0) > 1) regions, respectively.The colorscale of g (2) (0) is linear from g (2) (0) = 0 to 2 and logarithmic from g (2) (0) = 2 to 4, where it saturates.
tative agreement with the results reported in Ref. [20], where the authors also analyze the emission statistics of a thermally pumped CTS as a function of temperature and coupling strength.There are, however, some qualitative differences due to the description of the Hamiltonian and emission operators in the QRM [14][15][16] (Eqs.( 3)-(4)).For example, neglecting the Gauge correction of the QRM Hamiltonian result in an antibunched emission of the system for higher temperatures than in Fig. 11 (results not shown here).(0) FIG.
12. Intensity correlations calculated with different definitions of the intensity correlations.The blue line corresponds to the scenario where the detectors couple to the dressed operators of the photonic excitations in the cavity, xa (see Eq. ( 3)).The red and green line corresponds to the scenario where the detectors coupled to the first (Eq.(D2)) and second (Eq.(D4)) time derivatives of xa, respectively.All the calculations are done within the QRM for an intermediate pumping rate, Γ/γ = 10 −3 .

Appendix D: Other detection schemes
Several recent contributions (Refs.[46,47]) have discussed alternative formulations of the intensity correlations, as defined by the time derivatives of xa :  1)) and ω ν > ω µ their respective eigenvalues.This scenario where g (2) (0) depends on the ẋa operators describes an HBT interferometer with capacitive detectors that couple with the time derivative of the photons emitted by the cavity [46] (instead of coupling directly with the photons emitted in the cavity).
This coincidence rate can be also estimated by assuming that the two-photon emission is exclusively due to the relaxation from the |3−⟩ R state (see discussion in Section III and Eqs. ( 11) and ( 12)), and given by a product of its population (∼ 10 −10 , see Fig. 4) and emission rate taken as Cγ ≈ 2π ×2 10 14 s −1 .This approach yields a similar estimate of coincidence rate of 2π × 2 10 4 s −1 .
While this coincidence rate was estimated with generous assumptions about the collection efficiency, it remains several orders magnitude larger than the rates reported in contributions on the characterisation of statistics of faint emission from nonclassical emitters [48,49].

FIG. 2 .
FIG.2.Eigenvalues obtained as function of the coupling strength within the (a) JCM and (b) QRM.For each |ν⟩ (for the JCM) or |ν⟩ R (for the QRM) eigenstate, we show ∆E: the difference between its eigenenergy (ℏων ) and the ground state energy (ℏω0).
solid green line).The calculations shown in this figure are obtained with the QRM for Γ/γ = 10 −3 .

FIG. 7 .
FIG. 7. Landscape of g (2) (0) as a function of (a) normalized incoherent pumping rate Γ/γ, and (b) coupling parameter η = g/ω0.In (a) the collection of lines represents results for the coupling parameter varied linearly in the range η ∈ [0.1, 1].In (b) we exponentially increase the pumping rates from Γ/γ = 10 −6 to 10.The solid lines are obtained with the QRM, and the dashed orange and blue lines in (b) show the results obtained within the JCM for Γ/γ = 10 −6 and 10, respectively.

FIG. 8 .
FIG.8.Emission spectra from the USC system for the parameter η changing linearly from 0.1 to 0.7, in ∆η = 0.02 steps.The spectra are shifted vertically for clarity.Features describing transitions between eigenstates of the QRM Hamiltonian are traced and marked.The spectra are calculated within the QRM for an intermediate pumping rate Γ/γ = 10 −3 .

FIG. 10 .
FIG. 10. (a)-(d) Elements of the steady-state density matrices calculated for increasing values of the coupling parameter η = g/ω0 denoted at the top.The figure shows the logarithm of the absolute value of the first 6 × 6 elements of the density matrices (see Eq. (B1), where the columns and rows in the figure indicate the different combinations of |µ⟩R and R⟨ν| in Eq. (B1)).The results are color-coded, and the approximated numerical values for the largest terms are explicitly listed in the panels.The density matrices are obtained within the QRM for an intermediate pumping rate, Γ/γ = 10 −3 .