Persisting quantum effects in the anisotropic Rabi model at thermal equilibrium

Quantum correlations and nonclassical states are at the heart of emerging quantum technologies. Efforts to produce long-lived states of such quantum resources are a subject of tireless pursuit. Among several platforms useful for quantum technology, the mature quantum system of light-matter interactions offers unprecedented advantages due to current on-chip nanofabrication, efficient quantum control of its constituents, and its wide range of operational regimes. Recently, a continuous transition between the Jaynes-Cummings model and the Rabi model has been proposed by exploiting anisotropies in their light-matter interactions, known as the anisotropic quantum Rabi model. In this work, we study the long-lived quantum correlations and nonclassical states generated in the anisotropic Rabi model and how these indeed persist even at thermal equilibrium. To achieve this, we thoroughly analyze several quantumness quantifiers, where the long-lived quantum state is obtained from a dressed master equation that is valid for all coupling regimes and with the steady state ensured to be the canonical Gibbs state. Furthermore, we demonstrate a stark distinction between virtual excitations produced beyond the strong coupling regime and the quantumness quantifiers once the light-matter interaction has been switched off. This raises the key question about the nature of the equilibrium quantum features generated in the anisotropic quantum Rabi model and paves the way for future experimental investigations, without the need for challenging ground-state cooling.


I. INTRODUCTION
The resource theory of quantum information [1,2] demonstrates that quantum correlations found in nonclassical states are of utmost importance for performing quantum processing tasks [3,4].From a theoretical perspective, nonclassical states have proven highly valuable in investigating decoherence [5], witnessing the potential quantum nature of gravity [6,7], studying the quantumto-classical transition [8], enabling remote quantum control of the weak value amplification [9], and serving as a powerful tool for witnessing quantum phase transitions in critical systems [10].From a practical standpoint, quantum correlations have demonstrated their pivotal role as a fundamental element in emerging quantum technologies [11,12], including secure quantum communication [13,14], quantum sensing [15][16][17][18][19][20][21], and quantum simulation [22].Nonetheless, the intrinsic sources of noise and decoherence inherent in quantum dynamics render the production of long-lived quantum states necessary for quantum technologies [23,24].Efforts to advance efficient techniques for generating and protecting nonclassical states against quantum noise [25] encompass approaches such as decoherence-free subspaces [26][27][28], dynamical decoupling [29,30], and reservoir engineering schemes [31][32][33].Therefore, proposing schemes that facilitate the generation of long-lived highly correlated Quantum systems undergoing light-matter interactions stand as diverse platforms for generating nonclassical states and executing quantum information tasks [34,35].The recently achieved ultrastrong-coupling regime (USC) [36][37][38] and deep strong coupling regime (DSC) [39,40] of light-matter interactions have advanced quantum correlation generation beyond the strong coupling regime (SC) [41][42][43][44][45], i.e. when the light-matter coupling strength exceeds both the decay rates and the natural frequencies of the systems.Notably, a clear distinction between USC/DSC and the SC regime is the emergence of nonclassicality in the ground state of the system, manifested through squeezing and entanglement [46,47].Moreover, within the USC regime, thermal photons can exhibit anti-bunching behavior [48], and emissions from the matter subsystem (e.g., a two-level atom) can display photon bunching [49], whereas the behavior is reversed in the SC regime.
The cornerstone model describing light-matter interaction is composed of a two-level atom (qubit) coupled to a quantized single-mode cavity field through dipole interaction, known as the quantum Rabi model (QRM) [50][51][52].The QRM has been the subject of extensive theoretical work in quantum optics [53], for the generation of quantum entanglement [54], as an Otto quantum engine in quantum thermodynamics [55], to exhibit the emergence of quantum phase transitions [56][57][58], and motivated the development of a novel operational criterion for integrability [59,60].Interestingly, the QRM transforms into the Jaynes-Cummings model (JCM) [61] in situations where the natural frequencies of the system significantly surpass the bare light-matter coupling strength.In this scenario, one can invoke the rotating-wave approximation (RWA), permitting the retention of the rotatingwave (RW) terms while omitting the counter-rotatingwave (CRW) terms.However, this approximation becomes invalid once the system transit towards the USC regime [60,[62][63][64][65].
This motivates the introduction of a continuous transition between the JCM and the QRM by varying the RW and CRW coupling strengths, known as the anisotropic quantum Rabi model (AQRM).The AQRM has been explored in both closed systems [66,67] and open systems [68,69], with novel quantum technological schemes such as criticality-enhanced quantum sensing [70] and squeezed vacuum state laser [71].Notably, an exact solution for the AQRM with a biased term has been recently derived [72], and certain analytical solutions have been found via transcendental function extension [72,73] and transformation techniques [74,75].Experimental realizations of the AQRM have been achieved in twodimensional quantum wells [76], cavity quantum electrodynamics [77], and superconducting circuits [72].The flexibility to adjust the RW and CRW coupling strengths renders the AQRM to be highly versatile.Indeed, this grants the ability to explore a plethora of coupling regimes, thereby enabling the creation of highlycorrelated nonclassical states.A pivotal question arises: can we generate long-lived highly-correlated nonclassical states that persist even in the presence of a reservoir at finite temperature?This question is of paramount importance for current experimental platforms.
In this work, we demonstrate that long-lived highlycorrelated nonclassical states of the AQRM can indeed persist even in the presence of a reservoir at finite temperature.To support our theoretical findings, we use a quantum Markovian master equation to describe the evolution of the AQRM.This approach remains applicable across all coupling regimes, including degenerate points in the energy spectrum, i.e. without resorting to secular approximation.Consequently, such a quantum open system genuinely approaches a thermal Gibbs state as its steady state.Furthermore, through a series of complementary analyses, we characterize the quantumness of the steady state in the finite-temperaturecoupling strength phase diagram, showing that highlycorrelated nonclassical states can be achieved across a wide parameter region.
The rest of the paper is structured as follows: In Sec.II we present the open system description of our model.In Sec.III we briefly defined the measures and witnesses of nonclassicality thorughout our work.In Sec.IV, we outline our results, including the effects of temperature in the system.Finally, we present our concluding remarks in Sec.VI.
One key aspect of the AQRM in Eq. ( 1) is that preserves the parity symmetry Z 2 as in the isotropic model [38].This fact can be evidenced by considering the total number of excitations operator n = â † â + σ+ σ− and the parity operator π = exp (iπn).One can readily notice that, as opposed to the JCM, for the AQRM case [ Ĥ, n] ̸ = 0 and [ Ĥ, π] = 0.The above has two immediate consequences: (i) even in the ground state, the expected mean number of qubit and boson excitations is non-zero, and (ii) it implies that π possesses two eigenvalues, ⟨π⟩ = ±1, depending on whether the total number of excitations are even or odd.This in turns enable to solve the AQRM analytically by similar techniques used for the isotropic case, for instance, the Bargmann-Fock space of analytic functions and the Bogoliubov operator approach.
To gain a clearer understanding of the energy levels of the AQRM, in Fig. 1, we plot the energy spectrum in relation to the ground state energy as a function of the light-matter coupling λ 1 (λ 2 ) for several coupling ratios λ 2 /λ 1 (λ 1 /λ 2 ).As seen from the figure, the spectrum exhibits level crossings where the ground state parity changes sign (see green circles), leading the system to undergo infinite first-order quantum phase transitions as the coupling strength increases.In addition to the above degenerate crossing, several other non-degenerate and quasi-degenerate points can be observed.The analysis of the AQRM spectrum, as shown in Fig. 1, is important for understanding the physics underlying the emergence of thermal correlations, a study that we will address in later sections.
As stated before, the central objective of our work is to study the quantum correlations and quantum nonclassicality, which remain present even at thermal equilibrium -the so-called thermal correlations and thermal nonclassicality of the AQRM, respectively.See quantum thermalization and thermal entanglement in the open JCM [78] and QRM [79] as particular cases of the AQRM.To do so, we obtain the AQRM steady-state when such a system interacts unavoidably with a reservoir at temperature T .It is essential to point out that, even though 1. Energy level differences are plotted as functions of the light-matter coupling strengths λ1 (λ2) for different coupling ratios λ2/λ1 (λ1/λ2).In (c) and (h), the green circles represent the (degenerate) crossing of the energy levels and the plus and minus signs denote the change of parity of the energy levels.Other system parameters are specified as ∆ = 1 and ω = 1.
the standard Born-Markov master equation accurately describes the quantum state within the weak coupling regime λ i ≪ {ω, ∆}, it breaks down for stronger coupling regimes such as the USC and DSC regimes [37,38,80].Consequently, the steady-state at thermal equilibrium cannot be fully captured using the standard master equation in those regimes.To ensure that the steady-state is the actual quantum state at thermal equilibrium for an arbitrary set of coupling strengths λ i , one needs to switch to a Born-Markov master equation in a dressed picture [81], see Refs.[82,83] for a similar approach in the JC case.In this more suitable representation, the quantum jumps occur between the dressed eigenstates ladder of the full Hamiltonian as opposed to the standard master equation where the upward and downward jumps only consider the free energy of the Hamiltonian.We emphasize that, the AQRM spectrum (see Fig. 1) shows non-degenerate and degenerate eigenvalues which are appropriately modeled by the secular approximation in the derivation of the Markovian master equation, but it fails in the quasi-degenerate case.Following the steps of Ref. [84], we can derive a valid master equation for the entire spectrum of the Hamiltonian Eq. ( 1), as long as we consider weak damping, a high bath cut-off frequency, and a flat spectral density, these are sufficient conditions to ensure positivity and Markovianity.Indeed, the dissipative dynamics undergoes the following dressed master equation with Lindbladian superoperator term defined as In Eq. ( 2), |ϕ k ⟩ is the eigenvector of the AQRM, namely where E k is the kth eigenenergy associated to the eigenstate |ϕ k ⟩.We also define the dissipation rates with two explicit contributions, (i) The spectral function: where λ k,u accounts for the uth thermal bath coupled to a single-mode boson field with the frequency ω k , and energy gap (ii) And the transition coefficients: To satisfy the necessary conditions for the validity of the above master equation, we consider the Ohmic case: where α is the coupling strength between the system and the environment and ω c is the cutoff frequency of thermal baths.
Finally, the temperature of the bath is encoded in the Bose-Einstein distribution (k B = 1) It is important to emphasize that, when both reservoirs have the same temperature T a = T σ − (or only one of the subsystems is coupled to the reservoir) the system reaches a thermal equilibrium.Consequently, the dressed master equation steady-state solution of Eq. ( 2) results in the density matrix of the canonical ensemble [a straightforward numerical simulation proves that the Gibbs state is indeed the steady-state solution of Eq. ( 2)] where Z = k e −E k /T is the partition function and the steady-state population is Notice that, for unequal reservoir temperatures T a ̸ = T σ − , the steady solution of the dressed master equation cannot be written in terms of the canonical ensemble as in Eq. (12).In what follows, we assume equal temperatures reservoir such that we can always extract the statistical properties of the quantum state at thermal equilibrium from Eq. ( 12).

III. QUANTUM CORRELATION AND NONCLASSICALITY MEASURES
The concept of quantumness (nonclassicality) is related to the impossibility of describing physical phenomena by a deterministic or probabilistic classical theory.To understand the impact of the system's lightmatter anisotropy on the thermal quantum correlations and thermal nonclassicality, i.e. quantum correlations and nonclassicality that are present at thermal equilibrium, we thoroughly study the following quantities: (i) the zero-delay second-order correlation function g (2) (0), (ii) the bosonic field quadrature squeezing ζ 2 , (iii) a phase-space interference macroscopicity measure I(ρ) for the quantum state ρ, (iv) the negativity N (ρ), and (v) the quantum discord D(ρ).In the subsequent sections, we briefly explain the above quantities.

A. Zero-delay second-order correlation function
While the Poissonian and the super-Poissonian photon statistics of a light beam can be entirely explained in terms of a classical theory of light, the occurrence of sub-Poissonian photon statistics characterizes the quantumness of photonic states without a classical counterpart.It is also relevant to point out that, although this property does not consistently manifest in all quantum states of a field mode, a state can be classified as nonclassical when it is present.In this context, sub-Poissonian photon statistics serve as an authentic signature of the quantum nature of light.
Alternatively, one can classify the quantumness of a light beam by the study of the probabilities in measuring photons at a detector in a defined time interval t 2 − t 1 ≡ τ , the so-called second-order correlation function g (2) (τ ).This definition classifies the light beam in a threefold fashion, namely: (i) g (2) (0) > 1 bunched light where photons populate themselves together (classically a chaotic description), (ii) g (2) (0) = 1 random photon stream (classically a coherent description), and (iii) g (2) (τ ) < 1 antibunched light where photons distribute separately (with no classical analogy) [85,86], where we have considered an infinitesimal zero-delay time window τ → 0. Indeed, thermal photons emitted from non-interacting quantum modes or in the strong coupling regime is bunched (g (2) (0) = 2), being the standard example of an incoherent source of light [85][86][87].
The conventional definition of the normalized zerodelay second-order correlation function is [88] g (2) This quantity describes the probability of detecting two photons simultaneously (τ → 0), which is normalized by the probability of detecting two photons at once within a random photon source.Nonetheless, this definition holds for weak light-matter couplings, where the intracavity photons, described by â, suffice to explain the observed photons correlation.On the other hand, in the USC regime, where the qubit system strongly dresses the bosonic field, the second-order correlation function g (2) (0) is derived from the input-output formalism as [89,90] where with X− = ( X+ ) † , ∆ kj = E k − E j being the energy gap, and Here, X + jk describes the transition from the higher eigenstate |ϕ k ⟩ to the lower one |ϕ j ⟩.Notice that, in the weak light-matter interaction limit (i.e.λ i ≪ 1), the operator X+ is reduced to X+ = −iωâ.Thus, the correlation function in Eq. ( 15) simplifies to the conventional case.
One fundamental observation regarding the secondorder correlation function in Eq. ( 15) is that, in the eigenbasis of the AQRM Hamiltonian shown in Eq. ( 1), the dressed light-matter jump operator X provides the accurate expression for the average number of excitations in the ground state: ⟨ϕ 0 | X− X+ |ϕ 0 ⟩ = 0.This is in contrast to the seemingly incorrect result ⟨ϕ 0 |â † â|ϕ 0 ⟩ ̸ = 0. Note that as the dressed ground state can be spanned in the bare light-matter basis {|g⟩, |e⟩, |n⟩}, the dressed ground state will be composed of certain amount of virtual excitations [91].To convert the virtual excitations from the dressed picture into measurable photonic excitations, one can switch the interaction coupling strength on and off [91].
Indeed, consider a quantum state in the dressed basis |ϕ k ⟩.Once the interaction is switched off, the state will still contain the excitations that can now be detected using the intracavity photon operator â.Hence, ⟨ϕ k |â † â|ϕ k ⟩ ̸ = 0 will accurately describe the number of photonic excitations.It is worth noting that the nonadiabatic switch on/off needs to be of the order of the inverse of the relevant frequencies in the system, here ω, ∆ [91].
The conversion between virtual excitations and real photons via switching the interaction on and off is a consistent way to reconcile the correlation functions and other quantities that describe the thermal quantum correlations and nonclassicality in our work.Let us point out that the photons in the ground state of the USC models are virtual and cannot be detected [33], unless the coupling is suddenly switched off.
Throughout the paper, whenever a quantifier is defined using standard bosonic and spin operators, we refer to a scenario in which the coupling has been turned off after reaching thermal equilibrium.

B. Quadrature squeezing
The second-order correlation function fails as a nonclassical quantifier in the super-Poissonian regime.For instance, nonclassical Gaussian states, such as squeezed states may show either photon bunching or antibunching depending on whether the amplitude fluctuations are increased or reduced.To reveal the nonclassicality present in the steady state of the AQRM field mode, we compute the degree of squeezing.To do so, we introduce the generalized rotated field quadrature √ 2x θ = âe −iθ + â † e iθ , where any two operators that differ by θ = π/2 form a conjugate pair which satisfies the position (θ = 0) momentum (θ = π/2) commutator relation.We consider the following squeezing parameter [92]: Only a value of ζ 2 < 1 indicates bosonic squeezing, whereas ζ 2 = 1 corresponds to the photonic coherent state of radiation field.The minimization in the above definition can be easily obtained where we have assumed that once the steady-state is reached, the light-matter interaction is switched off.Note that, Refs.[93,94] studied the squeezing of output quadratures in the USC/DSC regime through the dressed light-matter jump operators for a vacuum input field, that is the thermal equilibrium case considered here.Interestingly, the authors found that any open system in its ground state cannot produce output squeezing, even if its ground state is a squeezed state.Since most of the squeezing in the AQRM is present in the ground state, these generalized measures leads to zero degree of squeezing in the thermal state if the coupling is on.For this reason we only focus in evaluating the degree of squeezing once the light-matter interaction is switched off, i.e., ζ 2 .

C. Phase space interference measure
To further study the nonclassicality of the AQRM steady-state, we consider a phase-space interference measure that quantifies the degree of macroscopicity of the state [95], defined as: where W (q, p) is the Wigner function [53].The macroscopicity I takes the value of 0 for classical states, while for pure quantum states such as superpositions of coherent states, NOON states, and Fock states, it corresponds to the average number of photons ⟨n⟩.In the case of the mixed state obtained by Eq. ( 12), the quantifier is constrained by the inequality I(ρ ss ) < Tr(ρ ss â † â).

D. Negativity
To explore the quantum correlations between the lightmatter constituents, we commence by calculating the quantum entanglement within bipartite systems.Among the various entanglement quantifiers available, we opt to use the negativity N (ρ) [96][97][98][99][100] where ρ T A is the partial transpose of the quantum state ρ with respect to subsystem A, ∥Y ∥ 1 = Tr|Y | = Tr √ Y † Y is the trace norm or the sum of the singular value the operator Y .Equivalently, the negativity can be computed as N (ρ) = 1/2 i (|ε i | − ε i ), where ε i are the eigenvalues of the partially transposed light-matter density matrix ρ.Note that N (ρ) = 0 corresponds to separable (not entangled) quantum states.

E. Quantum discord
Quantum entanglement is not the sole manifestation of quantum correlations.In fact, quantum discord (QD) represents another measure of potential quantum correlations within a quantum system.This measure of correlation emerges from the observation that two classically equivalent methods of defining mutual information yield different outcomes within the quantum domain [101].The QD of two subsystems A and B can be expressed as [101]: where S(ρ i ) = −Trρ i log ρi is the von Neumann entropy for the reduced density matrix, and S(ρ B|{Π A j } is the entropy conditioned through performing measurements on the A system, defined as: In the above, Π A j is a von Neumann projection operator on subsystems A and p j = Tr({Π A j } ⊗ I B ρAB ) is the probability with measurement outcome j.Note that, one of the most distinct differences between QD from entanglement is that QD can be non-zero for certain separable states [9,102,103].
g (2) (0) approaches unity regardless of the choice of coupling ratios.In other words, the thermal emission now displays the same statistical behavior as a coherent state.Note that in computing g (2) (0) using the bare bosonic field, we make the assumption that the coupling strength has been switched off [91] once the system has reached its thermal equilibrium.Notably, there is a significant distinction between the scenario in which the correlation function is computed from G (2) (0) and the one obtained when the light-matter coupling is switched off, resulting in g (2) (0).This subtlety demonstrates the intricate impact of virtual photonic excitations on the thermal state properties of the AQRM.
We can gain a clearer understanding of the secondorder correlation function G (2) (0) by considering the energy level crossing and its associated change in parity, as similarly discussed in Ref. [104].
Consider the particular case of G (2) (0) as a function of λ 1 for a fixed coupling ratio λ 2 /λ 1 = 0.5, as redrawn in Fig. 3(c) for clarity.The non-analytical (sharp) peaks result from the energy level crossing depicted in Fig. 3(a) (see Fig. 1 for other choices of parameters).Indeed, the closing of the energy gap between higher excited states, namely the parity change between the second (E 2 ) and the third (E 3 ) energy levels are responsible for the antibunching peaks shown in Fig. 3(c).Furthermore, the (0) 3. Relationship between energy level crossing and the behavior of the second-order correlation function G (2) (0).In panels (a) and (b), we redraw the energy gap En − E0 of the AQRM as a function of λ1 (λ2) for λ2/λ1 = 0.5 (λ1/λ2 = 0.5), see more cases of coupling ratios in Fig. 1.In panels (c) and (d), we show the particular case of G (2) (0) as a function of λ1 (λ2) for λ2/λ1 = 0.5 (λ1/λ2 = 0.5).The vertical gray lines evidence the correspondence between sharp peaks of G (2) (0) and the closing of energy gaps.Other system parameters are chosen as in Fig. 1.
closing of the energy gap between the ground state (E 0 ) and its first energy level (E 1 ) gives rise to the pronounced photonic bunched peak in Fig. 3(c) (around the value of λ 1 ∼ ω).Remarkably, the closing of the energy gap between E 1 and its ground state E 0 reopens as the coupling strength increases.This cycle of opening, closing, and reopening dynamics is responsible for the transition between antibunched-bunched-antibunched behavior of G (2) (0) shown in Fig. 3(c).A similar analysis can be conducted for Fig. 3(d), in which we present G (2) (0) as a function of λ 2 for a fixed coupling ratio λ 1 /λ 2 = 0.5.In contrast to the previous scenario, no non-analytic behavior is observed.This is because although a crossing of energy levels between the second and third eigenenergies occurs, the energy gap between the ground state and the first excited state monotonically vanishes as the coupling strength increases.As a result, a smooth transition towards a photonic bunching effect takes place.
To investigate the degree of photonic squeezing, in Fig. 4(a) [Fig.4(c)], we plot the squeezing parameter ζ 2 as a function of λ 1 [λ 2 ] for several coupling ratios λ 2 /λ 1 [λ 1 /λ 2 ] for a given temperature T = 0.1ω.As the figures show, there is a clear degree of squeezing ζ 2 < 1 of the AQRM steady-state (after the interaction has been switched off) for certain coupling ratios and a specific window of light-matter coupling.Note that for the JCM (λ 2 /λ 1 = 0), no squeezing is achieved, i.e. ζ 2 ≥ 1, for all values of λ 1 .Conversely, for the QRM (λ 2 = λ 1 ), the degree of squeezing tends to a photonic coherent field as the light-matter coupling strength goes towards the USC/DSC regime.Notably, no photonic squeezing occurs for the high anisotropy regime (λ i /λ j ≪ 1, i ̸ = j).For instance, consider λ 2 /λ 1 = 0.1 and λ 1 /λ 2 = 0.1 in Figs.4(a) and 4(c), respectively.As seen from the figures, no relevant photonic squeezing is achieved with the increasing of the light-matter coupling.This nonclassicality measure is in stark contrast with the second-order correlation function g (2) (0) (we recall that both cases are evaluated by turning the coupling off) for which it presents high degree of nonclassicality in the region of high anisotropy.
To further evidence quantum features of the steadystate AQRM field that may have not been captured by previous quantifiers, in Fig. 4(b) [Fig.4(d)], we show the macroscopicity interference-based measure I(ρ) as a function of λ 1 [λ 2 ] for different coupling ratios λ 2 /λ 1 [λ 1 /λ 2 ] and T = 0.1ω.As the figures show, the higher degree of interference-based measure generally occurs in the coupling region 1 < λ i /ω < 2 (i = 1, 2) for most of the coupling ratios of anisotropies, where quantum squeezing did not capture such a nonclassical behavior, yet it appears to be closer to capturing the second-order correlation quantifier.Interestingly, as shown in Fig. 4(b) I(ρ) increases with increasing anisotropy, thus exhibiting clear quantum signatures in the DSC regime.In contrast, in Fig. 4(d) I(ρ) is maximum for λ 1 /λ 2 = 0.5, which is a maximum of nonclassicality that is hidden in the previous measures.
Now, we study the thermal quantum correlations between the photons and the qubit through the quantifiers of negativity N (ρ) and quantum discord D(ρ).In Fig. 5(a), we plot the steady-state negativity N (ρ) as a function of λ 1 for several coupling ratios λ 2 /λ 1 for T = 0.1ω.As the figure shows, at some coupling ratios such as λ 2 /λ 1 = 0.1, 0.3, 0.5, 0.7, the negativity shows a clear dip near to zero between two maxima.This minima near zero corresponds to the energy level crossing of the ground-state and the first excited state, see Figs. 1(b)-(d).In the case of the JCM (λ 2 /λ 1 = 0), the steadystate negativity N (ρ) is nearly zero for coupling strengths λ 1 ≤ 0.7ω.However, as the coupling strength λ 1 increases, the negativity gradually rises until it reaches the highest degree of entanglement compared to any other coupling rate.In the case of the QRM (λ 2 /λ 1 = 1), N (ρ) exhibits a single maximum with the highest degree of entanglement for coupling strength λ 1 < ω.In Fig. 5(c), we depict the steady-state negativity N (ρ) as a function of λ 2 for different coupling ratios λ 1 /λ 2 at a fixed temperature of T = 0.1ω.As observed in the figure, the N (ρ) exhibits a smooth behavior without a sudden drop near zero.This is because, in contrast to the previous scenario, there is no energy level crossing between the ground state and the first excited state, as shown in Figs.1(g)-(i).While the steady-state negativity of Fig. 5(c) is qualitatively similar to the steady-state quantum discord D(ρ) of Fig. 5(d), the situation differs for Figs.5(a)-(b).Indeed, in Fig. 5(b), we plot the quantum discord D(ρ) as a function of λ 1 for several coupling ratios λ 2 /λ 1 for T = 0.1ω.As the figure shows, its overall profile is similar to that of Fig. 5(a), this means that almost all non-local quantum correlations are due to entanglement.However, unlike negativity shown in Fig. 5(a), at the crossing of the ground state and the first excited state the quantum discord D(ρ) does not decrease near to zero for some coupling strengths (e.g.λ 2 /λ 1 = 0.1, 0.3, 0.5, 0.7).Thus, we can infer: (i) all non-local quantum correlations are in the form of entanglement; (ii) the quantum correlations increase with anisotropy, and (iii) for high anisotropy they increase with the coupling strength presenting its maximum in the DSC regime.
It is worth noting that the system's behavior for increasing light-matter coupling, e.g., λ i ≳ 3ω (i = 1, 2), qualitatively resembles the system's behavior for the decoupled light-matter case, i.e., λ i = 0.Such a qualitative correspondence implies a trivial behavior of the quantifiers for very large light-matter interaction strengths, λ i ≳ 3ω (i = 1, 2), and for anisotropies with values exceeding λ i /λ j ≥ 0.3, where i ̸ = j, see Fig. 1.This behavior can be understood as, in the depths of the DSC regime, the spectrum becomes harmonic again (see Fig. 1 for λ i > 3ω), with the dressed states forming product states between displaced Fock photonic states and atomic states that are x-polarized, leading to G (2) (0) = 2.The nearly equal energy gaps (quasi-harmonic spectrum) and the eigenstates being product states in the DSC regime suggest that: i) no degree of squeezing and macroscopicity, ii) the vanishing of quantum entanglement and quantum discord at all temperatures should be expected (as shown in Figs. 2 to 5).

V. ROBUSTNESS OF THERMAL QUANTUMNESS AT HIGHER TEMPERATURES
As temperatures rise, thermal fluctuations increase, thereby modifying the AQRM steady-state as described in Eq. ( 12).On one hand, higher temperatures result in an increased expansion of states in the polariton basis, which may potentially lead to more frequent energy level crossings among higher AQRM eigenstates.On the other hand, these thermal fluctuations may cause the system to lose its coherence, potentially having a detrimental impact on the AQRM steady-state.To explore this non-trivial scenario, which involves a trade-off between increased energy level interplay and the influence of thermal fluctuations, we investigate how higher finite thermal bath temperatures affect the generation of nonclassical and quantum correlations in the AQRM, while considering all the above quantifiers.
In Fig. 6, we illustrate the second-order correlation function G (2) (0) as a function of the light-matter cou-FIG.6. Second-order correlation function G (2) (0) as a function of the light-matter coupling λi (i = 1, 2) for several choices of coupling ratios λi/λj (i ̸ = j) and different temperatures T .The black dashed horizontal line corresponds to G (2) (0) = 2.Other system parameters are as in Fig. 1.
pling λ i (i = 1, 2) for several choices of coupling ratios λ i /λ j (i ̸ = j) and different temperatures T .As evident from the figures, the higher the anisotropy, the more robust the quantum correlations are against temperature effects.Moreover, it exhibits a larger quantum region as a function of the coupling strengths.Notably, the correlation function G (2) (0) is prone to thermal decoherence, as signatures of nonclassicality (antibunching) survives only at temperatures as high as T ∼ 0.15ω, as depicted in all panels of Fig. 6.For higher temperatures, the antibunching effect is largely suppressed with increasing temperature for the same values of the coupling ratio, suggesting that thermodynamic fluctuations consistently destroy nonclassical features in high-temperature regions.
In Fig. 7, we show the temperature effects on the four quantifiers, namely ζ 2 , I(ρ), N (ρ), and D(ρ), as functions of the light-matter coupling λ i (i = 1, 2) and the temperature T for different coupling ratios λ i /λ j (i ̸ = j).
In the first row of Fig. 7, we plot the squeezing parameter ζ 2 as functions of λ i (i = 1, 2) and T for different λ i /λ j (i ̸ = j) values.Specifically, as observed in panels (1) and (4) of Fig. 7, no squeezing effect is achieved in the specific cases of the JCM (λ 2 = 0) and the AJCM (λ 1 = 0).However, for other coupling ratios, significant squeezing effects are observed at temperatures around T ≈ 0.3ω.Such a high degree of squeezing is attained by decreasing the anisotropy, making the Rabi model the ideal scenario for obtaining this resource.
In the second row of Fig. 7, we depict the macroscop-icity interference-based measure I(ρ) as functions of λ i (i = 1, 2) and T for different λ i /λ j (i ̸ = j) values.Unlike in the case of the squeezing parameter, the macroscopicity interference-based measure for the JCM and AJCM cases is non-zero, as shown in panels ( 8) and (11) of Fig. 7. It's worth noting that macroscopicity takes non-zero values over a broader range of coupling strengths compared to other ratios at low temperatures T < 0.1ω.However, as the figure shows, this quantifier is more susceptible to thermal fluctuations originating from the thermal reservoir and practically disappears for temperatures T > 0.1ω.Lastly, in the third and fourth rows of Fig. 7, we depict steady-state quantum entanglement N (ρ) and quantum discord D(ρ).Both quantifiers exhibit qualitatively similar behaviors.Similar to the macroscopicity interferencebased measure, both quantifiers are more pronounced over a broader range of coupling strengths and temperatures.As the ratio λ 2 /λ 1 (λ 1 /λ 2 ) increases, the range in which both quantifiers appear gradually narrows.In the case of 0 < λ 2 /λ 1 < 1, as λ 1 changes, quantum entanglement N (ρ) appears in two regions, with a value of zero between these two regions.This result is due to the energy level crossing between the ground state and the first excited state.However, for quantum discord D(ρ), only one region appears.Interestingly, in the case of high anisotropy, entanglement persists up to temperatures T ≈ 0.3ω, while quantum discord demonstrates greater robustness against thermal fluctuations up to temperatures T ≲ 0.4ω.
Our theoretical proposal to reach nonclassical equilibrium states is not just of academic interests and also meets the present-day nontrivial experimental challenges, mainly involving the coupling strength and low temperature environments.On the one hand, the physical model considered is general and has been implemented in a large number of systems, such as cavity quantum electrodynamics [77], and superconducting circuits [72].Among them, the USC regime has been experimentally demonstrated for more than a decade [36] (λ/ω > 0.1).Since then, the USC regime has been achieved in several other systems [38].On the other hand, the DSC regime is within the reach of experimental techniques, notably in Ref. [39] in circuit quantum electrodynamics setup implementing the QRM, for the highest reported strength of light-matter interaction of λ/ω ≈ 1.34 with temperature of approximately 45 mK (k B T /ℏω ≈ 0.16), where thermal entanglement was also confirmed.

VI. CONCLUSION
In this work, we consider the steady-state of the AQRM by solving a dressed master equation at thermal equilibrium.Such a master equation in the dressed basis is valid for any light-matter coupling strength, and where the resulting steady-state is indeed a Gibbs state.Our findings are two-fold: (i) we investigated the generation of quantum correlations and nonclassicality in the AQRM steady-state at thermal equilibrium.Through a comprehensive analysis of various measures of quantum correlation and nonclassicality, we demonstrate that quantum effects persist across a broad range of anisotropies and coupling strengths, even at moderate thermal equilibrium temperatures; and (ii) such quantum features significantly emerge in the USC regime and at the onset of the DSC regime, while practically vanishing in the deeper DSC regime.In the dressed picture, this results in virtual field excitations within the system.Notably, we show a significant difference in the second-order correlation function when computed using dressed jump operators compared to using bare field operators (i.e., when the light-matter coupling is assumed to be switched off after the system reaches thermal equilibrium).Both the difference between these quantifiers and the persisting quantum features of the AQRM steady-state at high temperatures, could allow for experimental approaches to study the impact of virtual excitations in the USC regime.In particular, as the quantum features emerging from the system avoid the need for demanding ground-state cooling.