Identifying unbound strong bunching and the breakdown of the rotating wave approximation in the quantum Rabi model

PHYSICAL REVIEW RESEARCH 5, 043213 (2023) Identifying unbound strong bunching and the breakdown of the rotating wave approximation in the quantum Rabi model Álvaro Nodar ,1,2,3Ruben Esteban,1,2,*Unai Muniain ,2Michael J. Steel ,3,† Javier Aizpurua ,1,2,‡and Mikołaj K. Schmidt 3,§ 1Materials Physics Center, CSIC-UPV/EHU, 20018 Donostia-San Sebastián, Spain 2Donostia International Physics Center DIPC, 20018 Donostia-San Sebastián, Spain 3School of Mathematical and Physical Sciences, Macquarie University, New South Wales 2109, Australia (Received 30 November 2022; accepted 5 November 2023; published 7 December 2023) 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 ultrastrong coupling regimes of light-matter interaction. We explore this transition using the intensity correlations g(2)(τ) of the emitted light, and we find strong, unbounded bunching of the emission from systems governed by the Rabi Hamiltonian. Surprisingly, this effect is observed not only in the ultrastrong 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 they 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. DOI: 10.1103/PhysRevResearch.5.043213 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 cavityTLS system (CTS) is dictated by the relationships between the coupling strength (g) and the characteristic resonant frequencies (ωσand ωa) and dissipation rates (γand κ)ofthe TLS and the cavity, respectively. In particular, if gis 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 gexceeds 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–5]. Large coupling strengths can lead the system into the socalled ultrastrong coupling (USC) regime, where gbecomes comparable to the resonant frequencies of the cavity and TLS [6]. The phenomenological limit for the onset of the USC is *[email protected] †[email protected] ‡[email protected] §[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. typically defined as g0.1ωaor g0.1ωσ(USC systems are commonly studied in a resonant configuration where ωa= ωσ=ω0). Importantly, systems operating in the USC regime exhibit new characteristics and nontrivial properties, such as the existence of a nonvacuum ground state [7,8], which would be largely missed if attempting to describe them using the typical 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–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 they might delineate a 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 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 longstanding questions about ensuring gauge invariance [14–19], and proposed a complete description of the interaction between a USC system and the environment. This formulation of the gauge-invariant QRM offers an opportunity to carefully 2643-1564/2023/5(4)/043213(13) 043213-1 Published by the American Physical Society ÁLVARO NODAR et al. PHYSICAL REVIEW RESEARCH 5, 043213 (2023) FIG. 1. (a) 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 Sec. II A) and the JCM (dashed lines; Sec. II B). For the strong incoherent pumping rate (/γ =10, orange lines), the results of the two models coincide in the WC regime (η2.5×10−2)and diverge 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. 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 Sec. II). In 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 Sec. II A for details). The value of ηat which the two models deviate depends 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 the more complete quantum Rabi model, and carefully analyzing their relationship. This work is structured as follows: In Sec. II we formally introduce the system under study, the quantum Rabi and Jaynes-Cummings Hamiltonians, and describe the formulation of the excitation, emission, and dissipation of the system. In Sec. III we identify the key mechanism which gives rise to the strong bunching, shown in Fig. 1(b).Finally, in Sec. IV we probe the extent of this effect in the USC, SC, and WC regimes, and we 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 (Sec. IIA) and Jaynes-Cummings (Sec. IIB) 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 ˆ H=¯hωaˆa†ˆa +¯hωσ 2{ˆσzcos[2η(ˆa+ˆa†)] +ˆσysin[2η(ˆa+ˆa†)]}.(1) This form of the Hamiltonian was first derived by Di Stefano et al. in Ref. [14], and it solves the critical problem of the gauge-invariance required from the Hamiltonian. 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 ˆadenotes the cavity photon annihilation, and ˆσz=ˆσˆσ†−ˆσ†ˆσand ˆσy=i(ˆσ†−ˆσ)are the zand yPauli 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 g(2)(τ)=I1(t+τ)I2(t) I1(t+τ)I2(t),(2) where I1(t+τ) and I2(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 043213-2 IDENTIFYING UNBOUND STRONG BUNCHING AND THE … PHYSICAL REVIEW RESEARCH 5, 043213 (2023) FIG. 2. Eigenvalues obtained as a 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 (¯hων) and the ground-state energy (¯hω0). reaches the steady state), this quantity is related to the statistics of photons inside the cavity as [12] g(2)(0) =ˆx† aˆx† aˆxaˆxass ˆx† aˆxa2 ss ,(3) where ˆ Oss denotes the expectation value of operator ˆ Oin the steady state (ss). The operator ˆxa= ν,μ;ων>ωμ|μRR μ|i(ˆa†−ˆa)|νRR ν|(4) 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 |μRbecause the emission from the cavity occurs due to the transitions between the eigenstates of the Hamiltonian (see Refs. [16,27] for an in-depth discussion of this formulation and the role of the secular approximation). Kets |νRand |μRin Eq. (4) are the eigenvectors of the QRM Hamiltonian, and ¯hων>¯hωμ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 Figs. 2(a) and 2(b)]. Note that throughout the text, we follow this naming convention for the eigenstates even after the eigenvalue crossing points (see the 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]: ∂ ∂tˆρ=Lˆρ =1 i¯h[ˆ H,ˆρ]+ 2Dˆx† σ[ˆρ]+γ 2Dˆxσ[ˆρ]+κ 2Dˆxa[ˆρ],(5) with Ldenoting the Lindbladian superoperator, and Dˆ O[ˆρ]=(2 ˆ Oˆρˆ O†−ˆ O†ˆ Oˆρ−ˆρˆ O†ˆ O) denoting the GoriniKossakowski-Sudarshan-Lindblad (GKSL) terms expressed through the dressed operators ˆxa[Eq. (4)] and ˆxσ[15]: ˆxσ= ν,μ;ων>ωμ|μRR μ|(ˆσ†+ˆσ)|νRR ν|.(6) The GKSL terms are introduced to describe the interaction between the TLS, cavity, and the bath following Refs. [15,19,28], including modeling the incoherent pumping of the TLS (with rate ), the dissipation of the TLS (γ) and of the cavity (κ). The dissipation is assumed to be additive, and we neglect any bath-mediated TLS-cavity interaction [7,16], and 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–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 ¯hω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 modeled 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=4g2/(κγ) in our study are C=0.08 (for η=10−3) and C=800 (for η=0.1). This framework is used to calculate the dependence of g(2)(0) [Eq. (3)] on the coupling parameter η=g/ω0, plotted in Fig. 1(b) as solid orange and blue lines, for the case of strong (/γ =10) and weak (/γ =10−6) pumping, respectively. All calculations in this paper have been carried out using the Python package QuTiP [39–41]. 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 ˆσzcos[2η(ˆa+ˆa†)] +ˆσysin[2η(ˆa+ˆa†)] =ˆσz+2ηˆσy(ˆa+ˆa†)+O(η2),(7) and we drop the terms nonlinear in η. Next, we introduce the rotating wave approximation (RWA) by removing the socalled non-number-conserving terms ˆσˆa+ˆσ†ˆa†to find the JCM Hamiltonian ˆ HJC =¯hω0ˆa†ˆa+¯hω0 2ˆσz+i¯hω0η(ˆσ†ˆa−ˆa†ˆσ).(8) 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 ˆaand ˆσoperators (instead of the dressed ˆxaand ˆxσoperators). The correlations arising in 043213-3 ÁLVARO NODAR et al. PHYSICAL REVIEW RESEARCH 5, 043213 (2023) this model can thus be calculated as g(2) JC =ˆa†ˆa†ˆaˆass ˆa†ˆa2 ss ,(9) where the steady state ˆρ(JC) ss is obtained from the standard master equation [42], ∂ ∂tˆρ(JC) =1 i¯h[ˆ HJC,ˆρ(JC)]+ 2Dˆσ†[ˆρ(JC)] +γ 2Dˆσ[ˆρ(JC)]+κ 2Dˆa[ˆρ(JC)],(10) 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 weakly coupled 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 at g(2)(0) <1 (see Appendix Afor 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(b) can be traced to characteristics (decay pathways and population) of the single |3−Reigenstate. To justify the focus on that particular eigenstate of the QRM Hamiltonian, in Fig. 3we plot the exact values of g(2)(0) (solid blue line) obtained for an intermediate pumping /γ =10−3together with approximated results. We start by approximating the steady state density matrix as being diagonal in the basis of the |νReigenstates of the QRM Hamiltonian, ˆρss ≈ ν Rν|νRR ν|.(11) The approximation is justified by noting that the system is driven with an incoherent pumping mechanism, and—much 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 the derivation in Appendix E) g(2)(0) ≈ν,μ Rν|Rμ|ˆxaˆxa|νR|2 ν,μ Rν|Rμ|ˆxa|νR|22,(12) where |μRand |νRare the eigenstates of the QRM Hamiltonian. The intensity correlations calculated by truncating the FIG. 3. Intensity correlations within several approximations: exact values of g(2)(0) (solid blue line); the diagonal approximation truncated to the states {|0R,|1−R,|1+R,|2−R,|2+R,|3−R} in Eq. (11) (solid orange line); and the diagonal approximation considering only the R3−| ˆxaˆxa|1−Rterm in the numerator of g(2)(0) in Eq. (12)(|3−R→|1−R(solid green line). The calculations shown in this figure are obtained with the QRM for /γ =10−3. double sum in the numerator up to |3−Rare 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−2is 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 limiting the double sum in Eq. (12) over νand μto the |ν=|3−Rand |μ= |1−Rterm: g(2)(0) ≈R3−|R1−| ˆxaˆxa|3−R|2 ν,μ Rν|Rμ|ˆxa|νR|22.(13) This approximation, denoted in Fig. 3with the green line, explores the role of the particular, correlated two-photon emission from |3−Rto the |1−Rstate. 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, nonthermal population of this state, expressed in our model as R3−. To analyze the latter effect in detail, we plot in Fig. 4the 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. 4shows that for larger coupling, the population R3−(orange solid line) rapidly grows at a rate proportional to η4, becoming larger than R2+(purple solid line), and eventually approaching R1+(red solid line) for η>0.1. This increase of R3−far exceeds that observed within the JCM (see the orange dashed line denoting the population of |3−). This large population of the |3−Rpolariton can be attributed to the new, direct excitation pathway introduced in theQRM[seeFig.5(a)]. In this model, the |3−Rpolariton can be directly driven from the ground state |0Rthrough the 043213-4 IDENTIFYING UNBOUND STRONG BUNCHING AND THE … PHYSICAL REVIEW RESEARCH 5, 043213 (2023) FIG. 4. Populations of the polaritonic states Rν(see labels) calculated with the QRM (solid lines) and JCM (dashed lines) for the intermediate pumping /γ =10−3. /2Dˆx† σincoherent pumping term introduced in the master equation (5). We illustrate this effect in Fig. 5(b) by plotting the |R3−| ˆx† σ|0R|2matrix element, which quantifies the efficiency of the direct excitation pathway, as a function of η.We identify a clear scaling with η4until η≈0.5, as in R3−. To dominate the statistical properties of the emission, the direct excitation pathway would have to 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 R3−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 R3−∝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 |R3−| ˆx† σ|0R|2, grows rapidly with η[see Fig. 5(b)]. IV. PROBING THE BREAKDOWN OF THE RWA In Sec. III we traced the bunching in the USC regime seen in Figs. 1(b) and 3to the direct excitation and two-photon emission from the |3−Rpolariton. Here, we explore the extent of this new effect, and we identify the threshold for the deviation between the JCM and the QRM. Figures 6(a) and 6(b) compare the values of g(2)(0) calculated with the QRM and the JCM, respectively, for a range of the normalized coupling strengths ηand pumping rates . As we discuss below, these results demonstrate that the deviation between JCM and QRM depends strongly on the latter parameter (), in a manner that can be explained using the formulation laid out in Sec. 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 FIG. 5. (a) Schematic and (b),(c) dependencies of the sequential and direct excitation pathways of the |3−Rpolariton in the QRM. (a) The two mechanisms are characterized by a cubic and linear dependence of the process on the incoherent driving rate , respectively. (b) Efficiency of the direct excitation of the |3−Rstate directly from the |0Rstate, quantified as |R3−| ˆx† σ|0R|2,growing approximately as η4(until η≈0.5), and matching the dependence of R3−showninFig.4. (c) Populations of the |3−Reigenstate as a function of incoherent pumping rate for a range of coupling parameters η(from bottom to top, η=10−3,5×10−3,10−2,5× 10−2,0.1,0.5,1). 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 we find that the intensity correlations generally follow g(2)(0) ∝1 .(14) 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 R3−∝(as discussed in Sec. III). Conversely, the denominator of Eq. (13) is dominated by the contribution from the |1±Reigenstates ∝R1±. Since both |1−Rand |1+Rare populated directly from the ground state, we find R1±∝. Thus we recover the g(2)(0) ∝R3−/(R1±)2∝1/ 043213-5 ÁLVARO NODAR et al. PHYSICAL REVIEW RESEARCH 5, 043213 (2023) FIG. 6. Dependence of g(2)(0) on the pumping rate and the coupling strength ηobtained within the (a) QRM and (b) JCM. The color map indicates in blue and red the antibunching [g(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. dependence, and identify that g2(0) is indeed unbounded as a function of . 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 point 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 the JCM (dashed lines) shown only for /γ =10−6and 10. 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−Reigenstate in the QRM. As we discuss in Sec. 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) 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−6to 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−6and 10, respectively. is determined by the population of the polariton with twophoton 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 ˆ HJC,inEq.(8) we (i) performed the expansion of the QRM Hamiltonian ˆ Hin 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 6shows 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 interpretation lays out an interesting 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 high-cooperativity coupling to the cavities, resulting in sub-ps timescales for the detection 043213-6 IDENTIFYING UNBOUND STRONG BUNCHING AND THE … PHYSICAL REVIEW RESEARCH 5, 043213 (2023) 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. of bunching (see discussion on experimental feasibility 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 a 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] S(1)(ω)∝∞ −∞ ˆx† a(τ)ˆxa(0)ss e−iωτ dτ, (15) where ˆxaand ˆx† aare 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 FIG. 9. Comparison between the emission spectra obtained with the QRM (solid lines) and JCM (dashed lines) for different couplings strengths, as indicated in the figure. Panels (a) and (b) correspond to different pumping rates, /γ =10−3and 10−6, respectively. 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. 8corresponds 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 |1−R→|0Rtransition increases, and two new peaks, corresponding to the |3−R→|2−Rand |2−R→|0Rtransitions (as labeled in the figure), emerge. 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 R3−and R2−, 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)]. To verify the sensitivity of S(1) to η, and the breakdown of the RWA, in Fig. 9we compare the spectra obtained from the JCM (dashed lines) and the QRM (solid lines) also using weaker incoherent pumping and various couplings η=10−3,10 −2, and 0.1 (note that η=10−3and 10−2are 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, which would be likely difficult to identify in experimental settings. On the other hand, as shown in Figs. 6and 043213-7 ÁLVARO NODAR et al. PHYSICAL REVIEW RESEARCH 5, 043213 (2023) 7(b),theg(2)(0) obtained for /γ =10−3and 10−6show qualitative differences for η2.5×10−2due to the direct excitation pathway discussed in Sec. 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−Rand |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 intensity correlations are a far more powerful tool than measuring the one-photon emission spectra for the identification of phenomena caused by the non-number-conserving terms of the QRM Hamiltonian. 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 call 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. ACKNOWLEDGMENTS Á.N., R.E., J.A. and M.K.S. thank A. González-Tudela and D. Burgarth for their useful inputs and valuable discussions. Á.N., R.E., U.M., and J.A. acknowledge the financial support from the Spanish Ministry of Science and Innovation and the Spanish government agency of research MCIN/AEI/10.13039/501100011033 through Project Ref. No. PID2019-107432GB-100, and from the Department of Education, Research and Universities of the Basque Government through Project Ref. No. IT 1526-22. Á.N. acknowledges funding from MCIN/AEI/10.13039/501100011033 and “ESF Investing in your future” through Project Ref. No. BES-2017-080073. M.K.S. acknowledges funding from the Macquarie University Research Fellowship Scheme (MQRF0001036), and the Australian Research Council Discovery Early Career Researcher Award DE220101272. APPENDIX A: CORRELATIONS IN THE JAYNES-CUMMINGS HAMILTONIAN In this Appendix, we derive an analytical expression for g(2)(0) as obtained with the approximated JCM Hamiltonian. We focus here on the same CTS system introduced in the main text, but we consider only the weak-illumination case, where the incoherent pumping of the TLS has a rate much smaller than the TLS losses, γ. We identify the minimum set of operators for which the master equations form an almost closed system (see the discussion below): v= (ˆa†ˆa,ˆσ†ˆσ,ˆa†ˆσ,ˆaˆσ†,ˆa†ˆaˆσ†ˆσ,ˆa†ˆaˆaˆσ†,ˆa†ˆa†ˆaˆσ,ˆa†ˆa†ˆaˆa)T. The equations of motion for the expectation values of these operators can be approximately expressed as d dt v=Mv+b,(A1) with M= ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −κ0−g−g00 00 0−γgg00 00 g−g−1 2(+κ+γ)0 −g000 g−g0−1 2(+κ+γ)−g000 00 0−(γ+κ)gg0 00 0 0 −2g−1 2(+γ+3κ)0 g 00 0 0 −2g0−1 2(+γ+3κ)g 00 0 0 0 −2g−2g−2κ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,(A2) 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) ˆa†ˆa†ˆaˆaˆσˆσ†≈ˆa†ˆa†ˆaˆa, and (iii) ˆa†ˆa†ˆaˆaˆσ†ˆσ≈0. In the steady state d dt vss =0, and we can derive closed expressions for ˆa†ˆass and ˆa†ˆa†ˆaˆass.Weare interested only in the values of ˆa†ˆass and ˆa†ˆa†ˆaˆa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 approximating the exact formulas for ˆa†ˆass and ˆa†ˆa†ˆaˆass assuming κγ: ˆa†ˆass ≈4g2(4g2+κ2) κ(γκ3+16g4+4g2κ2),(A3) 043213-8 IDENTIFYING UNBOUND STRONG BUNCHING AND THE … PHYSICAL REVIEW RESEARCH 5, 043213 (2023) FIG. 10. (a)–(d) Elements of the steady-state density matrices calculated for increasing values of the coupling parameter η=g/ω0denoted 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 |μRand 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. ˆa†ˆa†ˆaˆass ≈32g42 3κ2(γκ3+16g4+4g2κ2).(A4) The 2dependence of ˆa†ˆa†ˆaˆa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 Sec. III. The intensity correlations are found, in the two limits of interest, as g(2)(0) κγg −−−−−−→2 3 γ κ,(A5) g(2)(0) gκγ −−−−−−→2 3.(A6) Note that g(2)(0) is obtained here and in all the paper by neglecting the direct emission from the TLS, which is an invalid approximation for very small gin Eq. (A5) (unless the emission of the TLS is filtered out). If direct emission of the TLS is included, Eq. (A5) should be modified. The clearer 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) =0. APPENDIX B: APPROXIMATING THE STEADY-STATE DENSITY MATRIX In Sec. III, we calculate g(2)(0) considering several approximations of the density matrix of the steady state, such as truncating the Hilbert space to the |3−Rstate, and neglecting the off-diagonal elements. We found that these approximations are convenient for analyzing g(2)(0) beyond the WC regime (i.e., for η>2.5×10−2; see the discussion of Fig. 3). However, the steady-state off-diagonal elements can be crucial in some cases, like evaluating the QRM in the WC regime (see Fig. 3). In Fig. 10 we show the diagonal and off-diagonal Rμ νelements of the steady-state density matrix, ˆρss = μ,ν Rμ ν|μRR ν|,(B1) where |μRand |νRare the eigenstates of the QRM Hamiltonian. In the figure, we show Rμ νfor the first six eigenstates of the system, |μR,|νR= {|0R,|1−R,|1+R,|2−R,|2+R,|3−R}, and four coupling parameters η=0.001, 0.1, 0.25, and 1 (as denoted in the figure). All the density matrices show an approximate blockdiagonal form, with dominant Rn± n±and Rn± n∓terms. For the larger couplings [Figs. 10(b)–10(d)], ˆρss shows some additional coherence terms between eigenstates with the same parity of excitations, i.e., Rm± n±and Rm± n∓terms with mand n (m= n) being both even or odd numbers (e.g., the R3− 1−and R1+ 3−terms). As we increase η, the density matrices experience three main changes: (i) the ratio between the diagonal Rn± n± and off-diagonal Rn± n∓terms (Rn± n±/Rn± n∓) decreases, (ii) the Rn− n− population of the lower |n−Reigenstates becomes higher than the Rn+ n+population of the upper |n+Reigenstates, and (iii) the population of the higher-order eigenstates increases significantly. In particular, the population of |3−Rincreases from R3− 3−≈10−17 at η=0.001 to R3− 3−≈10−6at η=1, and becomes comparable to that of the lowest n=1 eigenstates. 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,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 Rν(Th.) ∝exp Eν kBT−1−1 ,(C1) where Eνis the energy of the νth eigenstate (calculated with respect to the energy of the ground state), kBis the Boltzmann constant, and Tis 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 to 3000 K; for reference, we also include a secondary vertical axis in the figure with kBT/(¯hω0) values]. For each temperature, we show the dependence of 043213-9