Spectral Classification of Coupling Regimes in the Quantum Rabi Model

# Spectral Classification of Coupling Regimes in the Quantum Rabi Model

## Abstract

The quantum Rabi model is in the scientific spotlight due to the recent theoretical and experimental progress. Nevertheless, a full-fledged classification of its coupling regimes remains as a relevant open question. We propose a spectral classification dividing the coupling regimes into three regions based on the validity of perturbative criteria on the quantum Rabi model, which allows us the use of exactly solvable effective Hamiltonians. These coupling regimes are i) the perturbative ultrastrong coupling regime which comprises the Jaynes-Cummings model, ii) a region where non-perturbative ultrastrong and non-perturbative deep strong coupling regimes coexist, and iii) the perturbative deep strong coupling regime. We show that this spectral classification depends not only on the ratio between the coupling strength and the natural frequencies of the unperturbed parts, but also on the energy to which the system can access. These regimes additionally discriminate the completely different behaviors of several static physical properties, namely the total number of excitations, the photon statistics of the field, and the cavity-qubit entanglement. Finally, we explain the dynamical properties which are traditionally associated to the deep strong coupling regime, such as the collapses and revivals of the state population, in the frame of the proposed spectral classification.

## I Introduction

The well-established Rabi model rabi1936 () describes the simplest class of light-matter interaction, the dipolar coupling between a two-level quantum system (qubit) and a classical monochromatic radiation field (unidimensional harmonic oscillator). In its quantum version, the radiation is specified by a quantized single-mode field, yielding the so-called quantum Rabi model (QRM) braak2011 (); solano2011 (). This model accurately describes the dynamics of a wide variety of physical setups, ranging from quantum optics to condensed matter systems book2013 (). In addition, a plethora of protocols in contemporary quantum information theory nielsen2004 (), with potential applications in future quantum technologies covering from ultrafast gates romero2012 () to quantum error correcting codes kyawprb2015 () or remote entanglement generation felicettiprl2014 (); rossatto2016 (), make use of the QRM as a building block. Therefore, the QRM plays an extremely important role in both theoretical and applied physics.

Typically, the standard experiments on cavity quantum electrodynamics (cavity QED) are restricted to a light-matter coupling strength much smaller than the natural frequencies of the unperturbed parts. Thus, they happen in the realm of the renowned Jaynes-Cummings (JC) model jc1963 (), which is obtained by applying the rotating-wave approximation (RWA) to the QRM knight2005 (). In this scenario, the achievement of the so-called strong coupling (SC) regime, when the coupling strength exceeds all decoherence rates, has driven the field of cavity QED for several decades book2013 (). Therefore, the JC model has served as a theoretical and experimental milestone in the history of quantum physics.

Since the last decade, a new coupling regime of the QRM, in which the coupling strength is a substantial fraction of the natural frequencies of the unperturbed parts, is being theoretically studied ciuti2005 (); bourassa2009 (); beaudoin2011 (); ballester2012 (); pedernales2015 (); todorov2010 () and experimentally reached in diverse solid state systems todorov2010 (); anappara2009 (); gunter2009 (); forndiaz2010 (); niemczyk2010 (); federov2010 (); muravev2011 (); schwartz2011 (); scalari2012 (); geiser2012 (); goryachev2014 (); zhang2016 (); chen2016 (); langford2016 (); braumueller2016 (). In this so-called ultrastrong coupling (USC) regime, the RWA is no longer suitable, such that the counter-rotating terms provide novel counterintuitive physical phenomena and new applications for the QRM emerge ashab2010 (); nataf2012 (); ridolfo2012 (); romero2012 (); ridolfo2013 (); stassi2013 (); felicetti2014 (); garziano2014 (); kyawprb2015 (); felicettiprl2014 (); garziano2015 (); kyaw2015 (); felicetti2015 (); rossatto2016 (); forndiaz2016 (); wang2016 (). When the counter-rotating terms can still be perturbatively treated, as in Refs. todorov2010 (); anappara2009 (); gunter2009 (); forndiaz2010 (); niemczyk2010 (); federov2010 (); muravev2011 (); schwartz2011 (); scalari2012 (); geiser2012 (); goryachev2014 (); zhang2016 (); chen2016 (), the QRM is approximately described by the Bloch-Siegert (BS) Hamiltonian forndiaz2010 (); beaudoin2011 (); klimovbook (). However, a few experiments have recently achieved the non-perturbative USC regime maissen2014 (); gambino2014 (); forndiaz2016-2 (); yoshihara2016 (); yoshihara2016-2 (), for which the full QRM has to be considered.

When the coupling strength is even stronger, surpassing the natural frequencies of the unperturbed parts, another regime of light-matter interaction appears, with totally different physics than the USC regime casanova2010 (); liberato2014 (). For this so-called deep strong coupling (DSC) regime casanova2010 (), the QRM can be reasonably described by an approximate solution as discussed in Refs. casanova2010 (); feranchuk1996 (); irish2005 (); irish2007 (). And, recently, F. Yoshihara et al. have experimentally achieved such an impressive coupling in superconducting circuits yoshihara2016 (); yoshihara2016-2 ().

Therefore, the advent of the aforementioned remarkable experimental and theoretical achievements has placed the QRM in the scientific spotlight. Nonetheless, the characterization so far established for the coupling regimes of the QRM is not quite universal, and a more specific criterion still remains undetermined. For instance, there are definitions stating that the USC regime is reached when the coupling strength is greater than a critical value related to either dynamical correlation functions wolf2013 () or quantum phase transitions domokos2015 (); jaako2016 (). However, for the latter case there is no consensus whether this transition can be reached in physical setups jaako2016 (); nataf2010 (); viehmann2011 (); ciuti2012 (); bamba2016 (), and this definition does not take into account the properties of the whole model, but only of its ground state. Another attempt was recently proposed in Ref. yoshihara2016-2 (), where the coupling classification is based on unique features exhibited in the transmission spectra of the system for different coupling regions. This approach uses the fact that the selection rules which allow or forbid transitions between eigenstates depend on the coupling value, changing the transmission pattern for different coupling regions. However, similarly to Refs. domokos2015 (); jaako2016 (), this approach does not take into account the properties of the whole model, since Another attempt was recently proposed in Ref. yoshihara2016-2 (), where the coupling classification is based on unique features exhibited in the transmission spectra of the system for different coupling regions. This approach uses the fact that the selection rules which allow or forbid transitions between eigenstates depend on the coupling value, changing the transmission pattern for different coupling regions. However, similarly to Refs. domokos2015 (); jaako2016 (), this approach does not take into account the properties of the whole model, since only the first four eigenstates are considered.

On the other hand, it is also conjectured in the literature that the USC regime is achieved when the coupling strength is just a substantial fraction of the natural frequencies of the unperturbed systems bourassa2009 (); todorov2010 (); forndiaz2010 (); niemczyk2010 (); scalari2012 (). Here, we are interested in quantitatively establishing how substantial this fraction has to be for the system description being significantly affected by the counter-rotating terms. Although the exact analytical solution of the QRM was recently presented for all parameter regimes braak2011 (), it strongly depends on zeros of a transcendental function defined through an infinite power series, making it difficult to extract the fundamental physics of that solution in general. Hence, it is more convenient to use approximate versions of the QRM as far as possible.

In this paper, we show that these approximate solutions are excellent guides to define a quantitative characterization of the coupling regimes of the QRM. In Sec. II, we show that the coupling regimes are naturally divided into three regions, whose boundaries depend not only on the ratio between the coupling strength and the natural frequencies of the unperturbed parts, but also on the energy to which the system can access. In addition, we show in Sec. III that our classification is supported by a completely different behavior of several static physical properties of the QRM, which depends on the region. Section IV provides a connection of our spectral classification with the dynamical properties that yield the traditionally blurry transition between the USC/DSC regimes. Finally, Section V comprises the conclusions of our work and the novel open questions emerging from it.

## Ii Coupling regimes of the quantum Rabi model

The Hamiltonian of the ubiquitous QRM is ()

 HR=ωa†a+Ω2σz+g0σx(a+a†). (1)

Here, are the Pauli matrices for the qubit, with transition frequency ( ground state and excited state), and () stands for the annihilation (creation) operator of a single-mode bosonic field, with frequency . The light-matter coupling is quantified by the vacuum-Rabi frequency .

### ii.1 Perturbative ultrastrong coupling regime

Whenever , with , , and , the QRM is well described by the JC model using the RWA knight2005 (),

 HJC=ωa†a+Ω2σz+g0(aσ++a†σ−), (2)

with . Paradigmatic examples of the intuitive physics behind the JC dynamics are the Rabi oscillations in the JC doublets, Eqs. (9) and (10) with , as a consequence of the conservation of the total number of excitations, and the collapses and revivals of the population inversion of the qubit knight2005 ().

When the counter-rotating terms can still be perturbatively treated, it is convenient to use the unitary transformation

 U=exp(Λ(aσ−−a†σ+)+ξ(a2−a†2)σz), (3)

with and . To second order in , this yields the Bloch-Siegert (BS) Hamiltonian forndiaz2010 (); beaudoin2011 (); klimovbook ()

 U†HRU≈H(2)BS=ω%BSσza†a+ωBSσz2−ωBS2+HJC, (4)

in which is the BS shift. From Eq. (4), it is straightforward to note that while the BS Hamiltonian provides the second-order correction in , the JC provides the zeroth-order one. Hence, the JC model is recovered from the BS Hamiltonian by imposing .

The energy spectrum of the BS Hamiltonian is

 EBS0 =−Ω2−ωBS, (5) EBSn,± =(n−12)ω−ωBS±12√(ΔBSn)2+4g20n, (6)

with and . The eigenstates are

 |ϕBS0⟩ =U|g,0⟩, (7) |ϕBSn,±⟩ =U|±,n⟩, (8)

with , where is the Fock state, and

 |+,n⟩=cosθn2|e,n−1⟩+sinθn2|g,n⟩, (9) |−,n⟩=sinθn2|e,n−1⟩−cosθn2|g,n⟩, (10)

in which the BS mixing angle is

 θn=arctan(2g0√nΔBSn). (11)

For the sake of simplicity, we will consider the resonant case () hereafter, but it is worth stressing that the following discussion is also suitable for the general case. In Fig. 1, we observe that the BS Hamiltonian provides an energy spectrum in excellent agreement with the one of the full QRM fidBS (), surprisingly up to the first energy-level crossings (the so-called Juddian points braak2011 ()).

Therefore, the use of the first Juddian points is an excellent attempt to define a boundary for a coupling regime. In this case, since the BS Hamiltonian perturbatively takes into account the breakdown of the RWA, we can establish the perturbative USC regime (pUSC) of the QRM as the region before the first Juddian points (), which is obtained by imposing . By squaring both sides of up to the elimination of the square roots, we end up with an eighth-degree polynomial in . Since the BS Hamiltonian is valid for perturbative values of , we truncate this polynomial up to second order, so that its non-negative solution yields the first Juddian points for each ,

 g×pUSCω≃1√2(2n+1), (12)

such that .

We can also notice from Fig. 1 that the more energetic the eigenenergies, the smaller the coupling values of first the Juddian points. This indicates that the importance of the counter-rotating terms depends not only on the ratio , but also on the energy to which the system can access, showing that the properties of the ground state are not sufficient to fully classify the coupling regimes of the QRM. Thus, the definition of the pUSC coupling regime is also connected to the energy to which the system can access.

Let us enlarge upon this point for the sake of clarity. The question we want to answer is whether a quantum state evolving under with a given will show features corresponding to the pUSC regime. This state is not necessarily an eigenstate, but it may be expanded in terms of eigenstates of the QRM. Thus, for a given , can have contributions from eigenstates in the pUSC region and from the region beyond that. Therefore, we take as a natural qualitative quantifier the mean energy of the state and choose the criterion that this state is in the pUSC regime when its energy is below the curve shown in the following.

If we invert Eq. (12) and replace in , assuming it as a continuous parameter, we can define the boundary of the pUSC regime as

 EpUSCω ≃14(ωg0)2[1−2(g0ω)4]−1 +14√[5−2(g0ω)2][1−2(g0ω)2]. (13)

This boundary is illustrated as the dotted line in Fig. 1, with the shaded area standing for the region where the perturbative USC regime is valid, i.e., when

 g0≲g×pUSCand¯E≲EpUSC. (14)

It is worth stressing that, besides the BS approach, there are other ones that can describe more accurately the QRM in a perturbative way feranchuk1996 (); irish2005 (); irish2007 (); jia2012 (). However, these methods result in a much more complicated solution for . We have also checked that the BS Hamiltonian expanded to third order in forndiaz2016 () provides more accurate eigenenergies, which also diverge from the exact calculated ones beyond the first Juddian points (see Appendix A). This indicates that the proposed definition for the pUSC region is not a simple consequence of the second-order term, but something deeper related to the breaking of the assumptions for the adiabatic expansion and the point from which the total number of excitation is no longer preserved, as discussed in Sec. III and shown in the left panel of Fig. 4(a).

### ii.2 Perturbative deep strong coupling regime

Analogously to the previous case, we can also employ the same ideas for the coupling regime at the other end, i.e., when the interaction term is no more a mere perturbation, but the main driver of the dynamics (DSC regime). In order to visualize the essence of this regime, it is convenient to rewrite in terms of the parity operator casanova2010 (), a conserved quantity of the QRM besides the total energy braak2011 (),

 HR=ωb†b+g0(b+b†)−Ω2(−1)b†bΠ, (15)

with . Since has eigenvalues and , there exists an independent Hamiltonian describing a perturbed displaced harmonic oscillator for each parity chain (), whose perturbation is given by the qubit term as an energy shift proportional to casanova2010 (). Thus, a perturbative approach up to first order in provides the following energy spectrum and the zeroth-order eigenstates (adiabatic approximation) irish2005 (); feranchuk1996 (); casanova2010 ()

 EpDSCn,± =(n−α2)ω±Ω2e−2α2Ln(4α2), (16) |ϕpDSC±,n⟩ =1√2[|+⟩⊗D(−α)|n⟩±|−⟩⊗D(α)|n⟩], (17)

in which is the Laguerre polynomial, , with , and . The energy of the ground state for this case is given by . A more refined approximation improves only marginally the accuracy of the eigenenergies and does not reveal new physical behavior wolf2013 ().

The DSC regime has a typical dynamical feature, which is the appearance of photon-number wave packets that bounce back and forth along a defined parity chain, yielding collapses and revivals of the initial population, even when the field and the qubit are initially in the vacuum and ground state, respectively casanova2010 (). This feature appears only for sufficient large values of , and it is more prominent after the last Juddian points (last energy-level crossings), when the adjacent eigenenergies asymptotically approach, becoming quasi-degenerate. Notice that the spectrum and the eigenstates of the QRM are described by Eqs. (16) and (17) with high fidelity wolf2013 (). Therefore, it is straightforward to note that the collapse-revival phenomenon is strictly related to the Schrödinger-cat-like states given by Eq. (17), which makes the perturbative solution an excellent attempt to define a boundary for a coupling regime.

In Fig. 2, considering , we notice that the energy spectrum given by Eq. (16) strongly agrees with the one of the full QRM beyond the last Juddian points, when the adjacent eigenenergies become quasi-degenerate. Thus, we use this fact to establish the boundary delimiting the perturbative DSC (pDSC) region. The boundary also connects with the appearance of the collapse and revivals of the initial population. For this purpose, let us consider the set of , with and , which are solutions of the equation

 1ω∣∣EpDSC+,n−EpDSC−,n∣∣=e−2α2∣∣Ln(4α2)∣∣≡δ, (18)

where is the maximum allowed energy difference close to the quasi-degenerate-energy region, which is related to the minimum fidelity allowed between the exact solution of the QRM and the perturbative states . Therefore, the limit of the pDSC region is given by the set of largest solutions of the above transcendental equation, .

For our calculations, we have chosen , value for which we have numerically observed better than fidelity between and the corresponding exact eigenstates of the QRM. The numerical solutions of Eq. (18) corresponding to the lowest values of are provided in Table 1. Inserting into , these points can be accurately fitted by the second-order equation

 EpDSCω+(g0ω)2≃a(g0ω)2+b(g0ω)+c, (19)

with , and , which is our definition of the boundary of the pDSC regime. Such a boundary is illustrated as the dotted line in Fig. 2, with the shaded area indicating the region where the perturbative DSC regime is valid, i.e., when

 g0≳g×pDSCand¯E≲EpDSC. (20)

If we change , then the change in the solution of the transcendental equation is approximately given by

 Δα=−e2α2δΔδ4αδ(Ln(4α4δ)+2L1n−1(4α2δ)), (21)

where is the solution for .

It is noteworthy to mention that we have numerically observed that the last maximum of the function is monotonically decreasing with . As a consequence, it could be that, for any fixed , there exists a value of such that the last solution to Eq. (18) could be placed before the last Juddian point. In any case, even in the situation in which this does not hold or when we want above this threshold, there are several strategies to overcome this complication. The most straightforward approach is to consider the analytical extension of the curve fitted for smaller . This actually works since this region is indeed perturbative. A second more complicated approach is introducing , with the fidelity of the perturbative eigenstates given by Eq. (17) in comparison with the exact eigenstates of the QRM. The problem with this approach relies on the lack of an analytically simple expression for the exact eigenstates of the QRM, and thus the fidelity can only be calculated numerically.

### ii.3 Non-perturbative ultrastrong/deep strong coupling regime

According to the aforementioned results, we are able to classify the coupling regimes of the QRM into three regions, whose boundaries depend not only on the relation between the coupling strength and the natural frequencies of the unperturbed parts, but also on the mean energy that the system can access, as summarized in Fig. 3. The pUSC regime belongs to the region right before the first Juddian points, whose physics are well described by the BS Hamiltonian, which still considers the interaction term as a perturbation. The BS Hamiltonian includes the well-known JC model, i.e., the QRM under the RWA. On the other hand, the pDSC regime belongs to the region beyond the last Juddian points, where there is a role interchange, since the interaction Hamiltonian becomes the main driver of the dynamics, while the bare Hamiltonian is the perturbative term. Finally, between these two coupling regimes, there is a region in which all parts of the Hamiltonian contribute on an equal footing to the dynamics. Then, we can establish this region as the non-perturbative USC (npUSC) regime, or as the non-perturbative DSC (npDSC) regime.

## Iii Static properties

Although our classification seems originally related to the validity of approximate mathematical models, it is associated to physical properties of the QRM which change their behavior qualitatively from region to region. In this section, we will focus on three relevant static properties, namely, the total number of excitations in the system, the photon statistics of the field, and the cavity-qubit entanglement.

In Fig. 4(a), we show the total number of excitations () for each eigenstate of the QRM as function of . In the left panel, we observe that the remains almost unchanged for the lower-coupling region of the pUSC regime (), which is expected since the BS Hamiltonian which governs the dynamics in this region commutes with . As we enter in the npUSC/npDSC regime, has a non-trivial oscillatory dependency with the coupling strength, which ceases as we approach the pDSC region (), as depicted in the right panel of Fig. 4(a). In the pDSC region, the total number of excitations becomes quasi-degenerate and increases with for the higher-coupling region of this region, with as shown in Appendix B.

Another physical property with characteristic behavior in each coupling region is the photon statistics of the field. Using the Fano-Mandel parameter given by Eq. (28) knight2005 (), we distinguish sub-Poissonian ( genuine nonclassical statistics), Poissonian (), and super-Poissonian () statistics. Such characteristic behavior is illustrated in Fig. 4(b), where we note that, except for the first two eigenstates that do not have an energy crossing, the field always exhibits a sub-Poissonian and a super-Poissonian photon statistics in the pUSC and pDSC regimes (see Appendix C), respectively, while all kind of photon statistics can be observed in the npUSC/npDSC regime. Moreover, we show in Appendix C that there are transitions in the photon statistics only in the npUSC/npDSC regime.

In Fig. 4(c), we observe that the entanglement between the qubit and the field (von-Neumann entropy of each subsystem knight2005 ()) also shows a peculiar behaviour, with the minima only appearing in the npUSC/npDSC regime. In addition, each minimum is always localized between two Juddian points and after the last one. The approximate analytical expressions of the qubit-field entanglement for the pUSC and pDSC regimes are given in Appendix D.

Our last remark is related to the decomposition of the field state in the Fock basis {} for the higher-coupling region of the pDSC regime, i.e., for . Using the terminology of the pDSC regime just to label the eigenstates, this decomposition is given by

 P(±,n)m=Tr(1q⊗|m⟩⟨m|⊗|ϕ±,n⟩⟨ϕ±,n|), (22)

in which and is the trace operation. We display in Fig. 5 considering , in which we firstly recognize that and tend toward the same multimodal distribution centered at . This can be confirmed by using Eq. (17), which predicts (red solid line in Fig. 5)

 P(±,n)m=α2|m−n|eα2min(m,n)!max(m,n)!(L|m−n|min(m,n)(α2))2, (23)

in which is the generalized Laguerre polynomial and . Secondly, we can also see that the number of modes of , , seems to be related to the number of energy crossings between and , which is .

It is worth to emphasize that the boundaries of our spectral classification does not imply an abrupt change in the physical properties of the QRM, as noticed in Fig. 4. Actually, such change gradually occurs around the boundaries of the pUSC and pDSC regions, even for dynamical properties as we will see in the next section.

## Iv Connection with dynamical properties

As already mentioned in Sec. II.2, the traditional characteristic signature of the DSC regime is not a static property, but a dynamical one. Namely, the appearance of photon-number wave packets that bounce back and forth along a defined parity chain, yielding collapses and revivals of the initial population. In this section, we show how the appearance of this phenomenon is related to our spectral classification.

In the upper panel of Fig. 6, we show the spectral classification together with the mean energy of two initial states, (solid line) and (dashed line), as function of . Considering the values of pointed out by the arrows, we computed the initial population . In the pUSC region, remains almost constant since is basically the ground state in that region [Fig. 6(a)], while exhibits Rabi oscillations due to the conservation of the total number of excitation [Fig. 6(b)]. As we enter in the npUSC/npDSC region, the Rabi oscillations pattern is lost, since the counter-rotating terms introduce a non-trivial oscillatory behavior in the initial population, as shown in Fig. 6(c)(d). Finally, as we approach the pDSC region [Fig. 6(e)(f)], the initial population start to present the collapse-revival pattern, which become more prominent as we go inside that region [Fig. 6(g)(h)].

## V Conclusion

In summary, we have introduced a spectral classification of the coupling regimes of the quantum Rabi model based on the validity of different perturbative approximations, showing that such regimes depend not only on the ratio between the coupling strength and the natural frequencies of the unperturbed parts, but also on the mean energy accesible by the system. Our classification is comprised by three coupling regions, namely the perturbative ultrastrong, the non-perturbative ultrastrong/deep strong and the perturbative deep strong coupling regimes. Remarkably, we have shown that the spectral classification is supported by a clearly divergent behavior of several relevant static physical properties in different coupling regimes. Additionally, we have also tested the suitability of our classification for the usual dynamical properties studied in the literature, which yield the traditional vague USC/DSC division. Therefore, our results clearly answer the long-standing question of providing a founded comprehensible classification of the coupling regimes in the QRM. Moreover, our results also open novel questions which motivate further studies of the mathematical and physical properties of these coupling regions, such as the physical role of the Juddian points in the QRM.

###### Acknowledgements.
We thank D. Braak, S. Felicetti, G. Romero, J. Casanova, and P. Forn-Díaz for fruitful discussions. This work was supported by the São Paulo Research Foundation (FAPESP) Grants No. 2013/04162-5, 2013/23512-7, and 2014/24576-1, Brazilian National Institute of Science and Technology for Quantum Information (INCT-IQ), CNPq, Spanish MINECO/FEDER Grant FIS2015-69983-P, Basque Government Grant IT986-16, and UPV/EHU UFI 11/55.

## Appendix A Influence of the higher orders of the BS approximation

In this Appendix, we discuss the influence of the third order of the BS expansion of the QRM in the definition of the pUSC region. First, let us consider the BS Hamiltonian expanded up to the third order in  forndiaz2016 ()

 H(3)BS =ωa†a+ω2σz−ωBS(σza†a+12) +g(^n)(a†σ−+aσ+), (24)

where

 g(^n)=g0(1−a†aωBS2ω) (25)

is the photon-dependent coupling strength. Notice that the Hamiltonian given by Eq. (A) preserves the number of excitations, i.e., . In Fig. 7, the exact eigenenergies of the QRM are depicted and compared with both second and third orders of the BS expansion. One can observe that, even though the third order is more accurate, it still diverges from the correct eigenenergies also after the first Juddian points. Thus, as mentioned in Sec. IIA, the proposed definition for the pUSC region is not a simple consequence of the second order term, but something deeper related to the breaking of the assumptions for the adiabatic expansion and the point from which the number of excitations starts to be not preserved, as we have seen in the left panel of Fig. 4(a).

## Appendix B Total excitations in the perturbative regimes

The total number of excitations is given by the operator . The fact that the number of excitations is preserved in the pUSC regime is a direct consequence of the commutation of this operator with the BS Hamiltonian (note that this also holds for ).

Let us now compute the mean value of the operator in the pDSC regime, i.e. . Hence,

with . In order to compute the first term, we use that . The second term is given by Eq. (D.2), so that , in which is the Kummer’s confluent hypergeometric function. By using that these functions are real valued, we obtain that

 ⟨^ne⟩pDSC−α2 =n+12±e−2α21F1(−n,1;4α2) =n+12±e−2α2Ln(4α2)≤δ, (27)

where we have made use of Kummer’s transformation  nist (). Therefore, effectively, the variation in is exponentially suppressed when , and upper-bounded by in the pDSC region.

## Appendix C Fano-Mandel parameter in pUSC and pDSC

In this Appendix, we compute the photon statistics of the QRM eigenstates in the pUSC and pDSC regimes through the Fano-Mandel parameter

 Q=⟨^n2⟩−⟨^n⟩2⟨^n⟩−1. (28)

Let us recall that the photon distribution is classified as sub-Poissonian ( genuine nonclassical statistics), Poissonian (), and super-Poissonian ().

### c.1 Perturbative USC regime

First, we must compute the photon distributions of the eigenstates , which is defined by . In order to perform the calculation, it is useful noticing that , with given by Eq. (3). By using the Baker-Campbell-Hausdorff formula to second order,

 U†amU =am+[am,H(α)] +12[[am,H(α)],H(α)]+O(α3), (29)

with,

 H(α)=(α/2)(aσ−−a†σ+)+(α2/4)(a2−a†2)σz, (30)

and .

It is straightforward to prove the useful expressions and , which may be used to compute, to the second order, the commutator

 [am,H(α) ]=−α2mam−1σ† −α24mσz(a†am−1+am−1a†)+O(α3). (31)

Let us now compute the second commutator of Eq. (C.1) to the second order

 [[am,H(α)],H(α)]=−α24mamσz+O(α3). (32)

Therefore, by replacing Eqs. (C.1) and (32) into Eq. (C.1), we obtain

 U†amU =am−α2mam−1σ†−α24mamσz −α24mσz(a†am−1+am−1a†)+O(α3). (33)

Now, we have to compute , also to the second order in , i.e., , which yields after normalization

 U†|g,0⟩ =(1−α28)|g,0⟩+α2|e,1⟩ −α2√24|g,2⟩+O(α3). (34)

By using this together with Eq. (C.1), we obtain

 U†|g,m⟩ =[1−(m+1)α28]|g,m⟩+α2√m+1|e,m+1⟩ −α24(√(m+1)(m+2)|g,m+2⟩ −√m(m−1)|g,m−2⟩). (35)

Analogously,

 U†|e,m⟩ =(1+mα24)|e,m⟩−α2√m|g,m−1⟩ +α24(√(m+1)(m+2)|e,m+2⟩ −√m(m−1)|e,m−2⟩). (36)

The scalar products of these states with respect to the state given by Eq. (9) yields

 ⟨g,m|U|+,n⟩ =(1+mα24)cos(θm+12)δn,m+1 −α2√msin(θm−12)δn,m−1 +α24[√(m+1)(m+2)cos(θm+32)δn,m+3 −√m(m−1)cos(θm−32)δn,m−3]. (37) ⟨e,m|U|+,n⟩ =[1−(m+1)α28]sin(θm2)δn,m −α24[√(m+1)(m+2)cos(θm+22)δn,m+2 −√m(m−1)cos(θm−22)δn,m−2]. (38)

The sum of the squares of these elements gives the photon distributions , for which we need to use the expressions and .

Now, we want to compute the first and second moments of the distribution

 ⟨^n⟩ =∞∑m=0mP+,nm=(n−12)−α√n5 +18(5−7n+2n2)+O(α3), (39) ⟨^n2⟩ =∞∑m=0m2P+,nm=(n2−n+12)+α4(1−2n)√n +α28(2n3−10n2+17n−5)+O(α3). (40)

In order to prove that the distribution is sub-Poissonian, it is sufficient to study the sign of

 ⟨^n2⟩−⟨^n⟩2−⟨^n⟩ =(34−n)+α√n4 −α216(4n3−8n2−13n+10), (41)

which can be straightforwardly proven to be negative in the pUSC regime, i.e., assuming that . The cubic polynomial is negative when and positive when . The first case can be directly checked. In the second case, , which finally proves that the photon distribution of the states is sub-Poissonian. In order to extend it to the states , it is only necessary to apply the substitutions and in Eqs. (37) and (38) and proceed analogously. This yields

 ⟨^n2⟩−⟨^n⟩2−⟨^n⟩ =(34−n)−α√n4 −α216(4n3−8n2−13n+10), (42)

which is also negative for . This concludes the proof.

### c.2 Perturbative DSC regime

Here, we prove that the photon distribution of the eigenstates of the QRM in the DSC regime is super-Poissonian. To achieve it, we proceed similarly to the previous subsection, assuming that in DSC the eigenstates are correctly described by Eq. (17). It is straightforward to see that

 ⟨g,m|ϕpDSC±,n⟩ =12(⟨m|D(−α)|n⟩∓⟨m|D(α)|n⟩), (43) ⟨e,m|ϕpDSC±,n⟩ =12(⟨m|D(−α)|n⟩±⟨m|D(α)|n⟩), (44)

where

 ⟨m|D(α)|n⟩=√m!n!e−12α2αm+n ×min(n,m)∑k=0(−1)n−k(m−k)!(nk)α−2k (45)

Let us start by computing the photon distribution for , which means that is given by

 Pm(ϕpDSCn,+)=∣∣ ∣∣√m!n!e−12α2αm+n2n∑k=0(−1)n−k(m−k)!(nk)α−2k∣∣ ∣∣2 ×((1−(−1)m+n)2+(1+(−1)m+n)2)=4=|⟨m|D(α)|n⟩)|2

Taking this into account, the computation of the first and second moments is straightforward, since

 ⟨^n⟩ =∞∑m=0mP+,nm=∞∑m=0⟨n|D†(α)|m⟩m⟨m|^n=a†aD(α)|n⟩ =n+α2, (46) ⟨^n2⟩ =∞∑m=0m2P+,nm=∞∑m=0⟨n|D†(α)|m⟩m2⟨m|^n2=(a†a)2D(α)|n⟩ =n2+α4+α2(4n+1) (47)

Therefore, in order to prove that the distribution is super-Poissonian for , we have to study the sign of

 ⟨^n2⟩−⟨^n⟩2−⟨^n⟩=n(2α2−1)>0, (48)

which proves it, since is always bigger than in the pDSC regime (see Table 1). For the case of the eigenstates , we only need to notice that the photon distribution is exactly the same, hence Eqs. (C.2)-(48) also hold, which concludes the proof.

## Appendix D Cavity-qubit entanglement

In this Appendix, we compute the cavity-qubit entanglement via the von-Neumann entropy in the pUSC and pDSC regimes, which allows us to analytically prove the numerical observations in Sec. III.

### d.1 Von-Neumann entropy in the pUSC regime

We have to compute the reduce density matrix for the qubit system. By using the Baker-Haussdorff-Campbell formula,

 U|±, n⟩⟨±n|U†=|±,n⟩⟨±n|+[H(α),|±,n⟩⟨±n|] +12[H(α),[H(α),|±,n⟩⟨±n|]]+O(H(α)3), (49)

in which with given by Eq. (30). The reduced density matrix is obtained by tracing out the bosonic degrees of freedom, i.e., . As usual, let us first consider the states . Then, the contribution to the reduced density matrix due to the first term in Eq. (D.1) is

 ρ(1) =Trcav(|+,n⟩⟨+n|) =cos2θn2|e⟩⟨e|+sin2θn2|g⟩⟨g|. (50)

For the second term, it is straightforward to prove that , so it will not be considered.

Finally, for the third term we must only consider the influence of the Hamiltonian term , since we are working in . Let us notice that the double commutator can be rewritten as , so let us compute both terms separately. The anti-commutator yields

 Trcav({|+,n⟩⟨+,n|,H (α)2})=−α24