# Spin fluctuations in the dissipative phase transitions of the quantum Rabi
model

Jiahui Li  Beijing Computational Science Research Center, Beijing 100193,
China  Rosario Fazio  The Abdus Salam International Center for Theoretical
Physics (ICTP), Strada Costiera 11, 34151 Trieste, Italy  Dipartimento di
Fisica, Università di Napoli Federico II, Monte S. Angelo, I-80126 Napoli,
Italy  Yingdan Wang<EMAIL_ADDRESS>Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Stefano Chesi<EMAIL_ADDRESS>Beijing Computational Science Research
Center, Beijing 100193, China  Department of Physics, Beijing Normal
University, Beijing 100875, China

###### Abstract

We investigate the dissipative phase transitions of the anisotropic quantum
Rabi model with cavity decay and demonstrate that large spin fluctuations
persist in the stationary state, having important consequences on the phase
diagram and the critical properties. In the second-order phase transition to
the superradiant phase, there is a significant suppression of the order
parameter and the appearance of non-universal factors, which directly reflect
the spin populations. Furthermore, upon entering a parameter regime where
mean-field theory predicts a tricritical phase, we find a first-order phase
transition due to the unexpected collapse of superradiance. An accurate and
physically transparent description going beyond mean-field theory is
established by combining exact numerical simulations, the cumulant expansion,
and analytical approximations based on reduced master equations and an
effective equilibrium theory. Our findings, compared to the conventional
thermodynamic limit of the Dicke model, indicate a general tendency of forming
extreme non-equilibrium states in the single-spin system, thus have broad
implications for dissipative phase transitions of few-body systems.

##  I Introduction

Coupling qubits to discrete bosonic modes is at the heart of quantum
information processing and quantum simulators. It allows to realize tunable
long-range interactions, useful to implement quantum gates or desired many-
body Hamiltonians. Besides cavity electrodynamics [1, 2], remarkable
collective phenomena such as dissipative phase transitions [3, 4], many-body
localization [5, 6], time crystals [7], and dynamical phase transitions [8, 9]
were explored in trapped ions [10] and circuit QED [11] platforms, where qubit
interactions are mediated by vibrational or cavity modes. Coupling multiple
spin qubits to superconducting resonators has also recently witnessed
substantial progress [12, 13].

Among these spin-boson systems, the Dicke model is one of the most actively
investigated paradigms of light-matter interactions [14, 15]. It describes $N$
qubits collectively coupled to a single-mode cavity which, for a sufficiently
strong interaction, undergoes a transition to a symmetry-broken superradiant
state. This superradiant phase transition can occur both in equilibrium [16,
17] and in the driven-dissipative case [18, 19], and has been observed in
various cold atoms and trapped ions experiments [20, 21, 22, 23, 24]. In
recent years, the strictly related concept of a superradiant phase transition
of the Rabi model [25, 26, 27, 28] (with $N=1$) has also attracted
considerable attention [29, 30, 31, 32]. In the absence of a conventional
thermodynamic limit, such transition takes place for a large value of the
ratio between the atomic and optical frequencies. The idea has been
generalized to a variety of setups, including models with multiple sites and
central-spin systems [33, 34, 35, 36, 37, 38, 39, 40, 41]. Remarkably, the
single-atom phase transition has been recently implemented with a trapped ion
coupled to its vibrational motion [30, 42, 43], and has been shown to have
great value for quantum metrology applications [44, 45, 46, 47].

In this article, we take a fresh look at the dissipative phase transition of
the Rabi model [32] and demonstrate the crucial role played by spin
fluctuations in correctly understanding the properties of the system. This is
an unexpected finding within this class of models, for which the semiclassical
treatment usually becomes exact in the appropriate thermodynamic limit. Here,
instead, spin fluctuations survive in the stationary state and significantly
suppress the order parameter of the superradiant phase. Their effect is
already large in the isotropic case of balanced co- and counter-rotating
interactions, when the spin is driven to an infinite-temperature state, but
can be further enhanced if counter-rotating terms dominate. Most dramatically,
we observe a sudden collapse of superradiance in a tristable region of the
phase diagram, where a superradiant state would be allowed within mean-field
theory. The disappearance of superradiance is due to an effective incoherent
pumping process of the spin, which ultimately leads to population trapping in
the normal state and gives rise to a first-order transition between the normal
and superradiant phases. Thus, spin fluctuations induce profound modifications
in the structure of the phase diagram.

The need of a careful treatment of quantum fluctuations is already clear from
mean-field theory, where multi-stability is found in the whole phase diagram.
For the Dicke model, bistability is only allowed in limited regions of
parameters [48, 49, 50], in the experimentally relevant [22] case of finite
anisotropy (i.e., unbalanced co- and counter-rotating interactions [18, 51,
52, 53, 54, 31, 49, 50]). Here, instead, the normal state with an excited
qubit is always stable in the appropriate ‘thermodynamic’ limit. Bistability
persists in the isotropic model and a tristable phase appears at large
anisotropy. The fate of multistablity is determined by quantum jumps between
mean-field solutions, thus requires an analysis of the stationary state based
on the full master equation.

After first establishing the above picture within mean-field theory, we
perform a preliminary study of quantum fluctuations based on the cumulant
expansion in the normal phase. We then provide a comprehensive and physically
transparent treatment of the stationary state and derive the corresponding
phase diagram. The stationary state of the normal phase is already non-
trivial, as it is determined by a competition between effective spin
relaxation and incoherent pumping. As mentioned, these processes bring the
qubit to an infinite-temperature state in the balanced model, and lead to
population inversion when counter-rotating terms dominate. In the superradiant
phase, the relaxation and pumping processes become more involved, but are
still amenable of analytical treatment and their effect can be summarized by
relatively simple transition diagrams (see Figs. 4 and 8). The dissipative
dynamics leads to a significant suppression of the order parameter from an
ideal supperradiant state (see, e.g., Fig. 5), due to the finite population in
the high-energy normal state. Identifying the structure of relevant
transitions between semiclassical states not only allows us to describe
accurately the superradiant order parameter, but also leads to a clear
explanation for the absence of superradiance in the tricritical region (see
Fig. 8).

Furthermore, we show that at the second-order phase transition an effective
equilibrium theory originally developed for the isotropic Dicke model [55] can
be extended to the anisotropic Rabi model. This approach allows us to derive
accurate analytic formulas for the scaling functions and the non-universal
factors relating the two (Rabi and Dicke) models. In particular, for the
isotropic case the exact value of the non-universal factor is 1/2, which
agrees with previous numerical findings [32]. While the precise origin of
these non-universal factors has not been fully appreciated, our analysis shows
that they coincide with the spin populations. This simple interpretation links
in a direct manner these prefactors of the scaling behavior to the nontrivial
spin dynamics taking place in the Rabi model.

The structure of our paper is as follows: In Sec. II we define the model and
in Sec. III we obtain the phase diagram within mean-field theory. In Sec. IV
we treat the leading quantum corrections induced by a finite value of $\eta$.
In particular, in Sec. IV.1 we use the cumulant expansion to capture quantum
fluctuations in the normal phase. An alternative approach based on an
effective master equation is introduced in Sec. IV.2, which is later
generalized in Sec. IV.3 to the superradiant phase and in Sec. IV.4 to the
tricritical phase. This allows us to characterize the stationary state in
general, and derive the corresponding phase diagram. Finally, in Sec. V we
show that a modified effective equilibrium theory is able to describe
accurately the critical scaling of the second-order phase transition. Our
conclusions are presented in Sec. VI.

##  II Model

The Hamiltonian of the anisotropic Rabi model reads:

$\displaystyle\hat{H}=$
$\displaystyle\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}-\lambda_{-}\left(\hat{a}\hat{\sigma}_{+}+\hat{a}^{\dagger}\hat{\sigma}_{-}\right)$
$\displaystyle-\lambda_{+}\left(\hat{a}\hat{\sigma}_{-}+\hat{a}^{\dagger}\hat{\sigma}_{+}\right),$
(1)

where $\hat{a}$ is the annihilation operator of the photon field,
$\hat{\sigma}_{x,y,z}$ are the Pauli matrices,
$\hat{\sigma}_{\pm}=(\hat{\sigma}_{x}\pm i\hat{\sigma}_{y})/2$, $\omega_{c}$
is the frequency of the bosonic mode, and $\Omega$ is the energy splitting of
the two-level system (we set $\hbar=1$). Here the interaction is written in
terms of the couplings $\lambda_{-}$ and $\lambda_{+}$ corresponding to co-
rotating and counter-rotating terms, respectively. Alternatively, the
notation:

$\lambda_{\pm}=\lambda_{x}\pm\lambda_{y},$ (2)

is sometimes preferred, which emphasizes the coupling strength $\lambda_{x}$
($\lambda_{y}$) of $\hat{J}_{x}$ ($\hat{J}_{y}$) to the
$\hat{a}+\hat{a}^{\dagger}$ ($i\hat{a}-i\hat{a}^{\dagger}$) quadrature of the
field [51, 49]. The imbalance of the corotating and counter-rotating terms can
be altered by the power of a bichromatic laser drive in trapped ion setups, or
depends on the relative strength of the electric and magnetic couplings in
circuit-QED implementations [56, 30, 57, 58]. Considering an open environment
causing photon decay, the model is described by the master equation:

$\displaystyle\dot{\hat{\rho}}=-i[\hat{H},\hat{\rho}]+\kappa\mathcal{D}[\hat{a}]\hat{\rho},$
(3)

where
$\mathcal{D}[\hat{a}]\hat{\rho}=2\hat{a}\hat{\rho}\hat{a}^{\dagger}-\hat{a}^{\dagger}\hat{a}\hat{\rho}-\hat{\rho}\hat{a}^{\dagger}\hat{a}$.
The above quantum evolution can be alternatively expressed as:

$\displaystyle\frac{d}{dt}\hat{a}=-\left(\kappa+i\omega_{c}\right)\hat{a}+i\left(\lambda_{-}\hat{\sigma}_{-}+\lambda_{+}\hat{\sigma}_{+}\right)+\hat{\mathcal{F}},$
$\displaystyle\frac{d}{dt}\hat{\sigma}_{+}=i\Omega\hat{\sigma}_{+}+i\left(\lambda_{-}\hat{a}^{\dagger}+\lambda_{+}\hat{a}\right)\hat{\sigma}_{z},$
$\displaystyle\frac{d}{dt}\hat{\sigma}_{z}=2i\left[\lambda_{-}\left(\hat{a}\hat{\sigma}_{+}-\hat{a}^{\dagger}\hat{\sigma}_{-}\right)-\lambda_{+}\left(\hat{a}\hat{\sigma}_{-}-\hat{a}^{\dagger}\hat{\sigma}_{+}\right)\right],$
(4)

where the force $\hat{\mathcal{F}}(t)$ is a stochastic operator satisfying
$\langle\hat{\mathcal{F}}(t)\hat{\mathcal{F}}^{\dagger}(t^{\prime})\rangle=2\kappa\delta(t-t^{\prime})$
and
$\langle{\hat{\mathcal{F}}}^{\dagger}(t^{\prime})\hat{\mathcal{F}}(t)\rangle=0$.

In the rest of the paper, we will assume $\omega_{c}>0$ and, unless explicitly
stated, $\Omega,\lambda_{\pm}>0$. To see that the latter restriction is not
important, we first consider
$\hat{U}_{1}=\exp[i\frac{\pi}{4}(\hat{\sigma}_{z}-2\hat{a}^{\dagger}\hat{a})]$,
which transforms the couplings as $\lambda_{\pm}\to\pm\lambda_{\pm}$. Instead,
$\hat{U}_{2}=\exp[i\frac{\pi}{2}\hat{\sigma}_{z}]$ leads to a change of sign
of both couplings, $\lambda_{\pm}\to-\lambda_{\pm}$. Through $U_{1}$ and
$U_{2}$ we can always map the system to the region $\lambda_{\pm}>0$.
Furthermore, $\hat{U}_{3}=\exp[i\frac{\pi}{2}\hat{\sigma}_{x}]$ changes the
sign of the atomic level splitting, $\Omega\to-\Omega$, while interchanging
the role of co-rotating and counter-rotating terms,
$\lambda_{\pm}\to\lambda_{\mp}$. Therefore, we can also assume $\Omega>0$. We
will also often use dimensionless couplings
$\tilde{\lambda}_{\pm}=\lambda_{\pm}/\lambda_{c}$, where:

$\lambda_{c}=\frac{1}{2}\sqrt{\omega_{c}\Omega}$ (5)

is the critical coupling to the superradiant phase in the closed ($\kappa=0$)
and isotropic ($\lambda_{-}=\lambda_{+}$) limit, both for the Dicke and Rabi
model. However, while the usual quantum phase transition (QPT) of the Dicke
model takes place for $N\to\infty$, for the Rabi model one should consider the
following frequency ratio:

$\eta=\frac{\Omega}{\omega_{c}}.$ (6)

and take the limit $\eta\to\infty$, which plays the same role of the
conventional thermodynamic limit [25, 26, 27, 28, 29, 30, 31, 32].

##  III Mean-field analysis

We fist examine the dissipative phase transition through standard mean-field
theory, which is normally expected to become accurate in the thermodynamic
limit $\eta\to\infty$. Here, however, since we are also interested in the
stationary state of the system, the regime of validity should be carefully
specified. As we will see, the problem arises due to the presence of
multistability at arbitrary values of $\lambda_{\pm}$. In short, mean-field
theory is valid when taking the limit $\eta\to\infty$ from the onset, i.e.,
before considering $t\to\infty$. This effectively suppresses the transitions
between the stationary mean-field solutions, making them infinitely long-
lived. A more precise discussion is given at the end of this section.

From Eq. (II), mean-field equations for the expectation values
$a=\langle\hat{a}\rangle$, $s_{\pm}=\langle\hat{\sigma}_{\pm}\rangle$, and
$s_{z}=\langle\hat{\sigma}_{z}\rangle$ are obtained by neglecting quantum
fluctuations, i.e., assuming the factorization
$\langle\hat{A}\hat{B}\rangle=\langle\hat{A}\rangle\langle\hat{B}\rangle$. The
system conserves the pseudo-angular momentum, thus we set
$4|s_{+}|^{2}+s_{z}^{2}=1$. Here, as we are interested in the limit
$\eta\to\infty$ instead of $N\to\infty$, we find that the following rescaling
of the photon field and coupling strengths are particularly appropriate:

$c=\frac{\langle\hat{a}\rangle}{\sqrt{\eta}},\qquad\tilde{\lambda}_{\pm}=\frac{\lambda_{\pm}}{\lambda_{c}},$
(7)

where $\eta$ is defined in Eq. (6) and $\lambda_{c}$ is given in Eq. (5).
After also rescaling the time as $\tau=\omega_{c}t$, we find the mean-field
equations:

$\displaystyle\frac{dc}{d\tau}=-\left(\frac{\kappa}{\omega_{c}}+i\right)c+\frac{i}{2}\left(\tilde{\lambda}_{-}s_{-}+\tilde{\lambda}_{+}s_{+}\right),$
(8) $\displaystyle\frac{1}{\eta}\frac{ds_{+}}{d\tau}=i\,{\rm
sgn}(\Omega)\,s_{+}+\frac{i}{2}\left(\tilde{\lambda}_{-}c^{\ast}+\tilde{\lambda}_{+}c\right)s_{z},$
(9)
$\displaystyle\frac{1}{\eta}\frac{ds_{z}}{d\tau}=i\tilde{\lambda}_{-}\left(cs_{+}-c^{\ast}s_{-}\right)-i\tilde{\lambda}_{+}\left(cs_{-}-c^{\ast}s_{+}\right).$
(10)

The stationary states of these semiclassical equations are independent of
$\eta$ and, actually, coincide with the stationary states of the Dicke model
[49, 50]. As we will see, however, an important difference appears in the
stability analysis. Since here we are interested in the limit $\eta\to\infty$,
the left-hand-side of Eqs. (9) and (10) vanish, which allows us to solve these
equations explicitly. This reflects the fact that the spin follows
adiabatically the instantaneous value of $c$ when
$\Omega\gg\omega_{c},\lambda_{\pm}$. Using the normalization condition
$s_{z}^{2}+4|s_{+}|^{2}=1$ and assuming $\Omega>0$ we obtain:

$\displaystyle
s_{z}=\pm\frac{1}{\sqrt{1+|\tilde{\lambda}_{-}c^{*}+\tilde{\lambda}_{+}c|^{2}}},$
(11) $\displaystyle
s_{+}=-\frac{1}{2}(\tilde{\lambda}_{-}c^{*}+\tilde{\lambda}_{+}c)s_{z},$ (12)

where the $\pm$ sign is determined by the initial condition of $s_{z}$. In
this limit, spin states in the upper (lower) half of the Bloch sphere never
cross to the lower (upper) side. Finally, Eqs. (8–10) can be written as an
equation of motion for the cavity field alone:

$\frac{dc}{d\tau}=-\left(\frac{\kappa}{\omega_{c}}+i\right)c\mp\frac{i}{4}\frac{(\tilde{\lambda}_{-}^{2}+\tilde{\lambda}_{+}^{2})c+2\tilde{\lambda}_{-}\tilde{\lambda}_{+}c^{*}}{\sqrt{1+|\tilde{\lambda}_{-}c^{*}+\tilde{\lambda}_{+}c|^{2}}}.$
(13)

Figure 1:  Mean-field phase diagram at $\eta\to\infty$. Only stable stationary
solutions are indicated. In the coexistence region C the three states (NP↑,
NP↓, and SP↓) are all stable. The black line indicates the isotropic case
($\tilde{\lambda}_{-}=\tilde{\lambda}_{+}$), where the asterisk marks the
critical point $\tilde{\lambda}_{c}^{*}=\sqrt{1+\kappa^{2}/\omega_{c}^{2}}$.
The empty circle marks one of the tricritical points, with
$\lambda_{-}=\tilde{\lambda}_{c}^{\rm T}$ [see Eq. (21)]. Here we used
$\kappa/\omega_{c}=0.5$.

Based on the mean-field equations, we obtain the phase diagram shown in Fig.
1. While at finite $\eta$ the long-time dynamics may be characterized by an
oscillatory or chaotic behavior, resulting in a complex phase diagram [50],
the mean-field treatment becomes much simpler when $\eta\to\infty$. Here, we
only need to consider the normal and superradiant fixed points. In the normal-
phase (NP), the bosonic field has zero mean-field amplitude while the qubit
can be either in the $|\downarrow\rangle$ state

$\displaystyle{\rm{NP}_{\downarrow}}:s_{z}=-1,\,\,\,s_{+}=0,\,\,\,c=0,$ (14)

or fully populated in the up-state $|\uparrow\rangle$:

$\displaystyle{\rm{NP}_{\uparrow}}:s_{z}=1,\,\,\,s_{+}=0,\,\,\,c=0.$ (15)

The stability analysis of the NP states can be performed by linearizing Eq.
(13) around $c=0$. It is then found that the NP↑ state is always stable, while
NP↓ becomes unstable in the purple region of Fig. 1. It is convenient to
specify the stability condition for a fixed ratio
$r=\tilde{\lambda}_{+}/\tilde{\lambda}_{-}$. Then, NP↓ is unstable for

$\tilde{\lambda}_{c}^{-}<\tilde{\lambda}_{-}<\tilde{\lambda}_{c}^{+},$ (16)

where we defined:

$\tilde{\lambda}_{c}^{\pm}=\frac{2\sqrt{1+r^{2}\pm\sqrt{4r^{2}-(1-r^{2})^{2}\kappa^{2}/\omega_{c}^{2}}}}{|1-r^{2}|}.$
(17)

In Fig. 1, the solid phase boundary of the purple region is at
$\tilde{\lambda}_{-}=\tilde{\lambda}_{c}^{-}$ (while
$\tilde{\lambda}_{+}=r\tilde{\lambda}_{c}^{-}$) and the two dotted boundaries
are at $\tilde{\lambda}_{-}=\tilde{\lambda}_{c}^{+}$. Note that in Eq. (17)
$r$ should satisfy $r_{-}\leq r\leq r_{+}$, where:

$r_{\pm}=\sqrt{1+\omega_{c}^{2}/\kappa^{2}}\pm\omega_{c}/\kappa.$ (18)

When $r<r_{-}$ or $r>r_{+}$, both NP↓ and NP↑ are stable fixed points for
arbitrary values of $\tilde{\lambda}_{-}$.

Besides NP states, the following superradiant fixed-point plays an important
role:

$\displaystyle{\rm{SP}}_{\downarrow}:$
$\displaystyle\,\,s_{z}=-\cos{\theta_{\rm SP}},\quad
s_{+}=\pm\frac{1}{2}e^{i\phi_{\rm SP}}\sin\theta_{\rm SP}$ $\displaystyle
c=\frac{\tilde{\lambda}_{-}s_{-}+\tilde{\lambda}_{+}s_{+}}{2(1-i\kappa/\omega_{c})},$
(19)

where the direction of the spin is given by:

$\displaystyle\cos{\theta_{\rm
SP}}=\frac{4\left(1+\kappa^{2}/\omega_{c}^{2}\right)}{\tilde{\lambda}_{-}^{2}+\tilde{\lambda}_{+}^{2}+\sqrt{4\tilde{\lambda}_{+}^{2}\tilde{\lambda}_{-}^{2}-\frac{\kappa^{2}}{\omega^{2}_{c}}\left(\tilde{\lambda}_{-}^{2}-\tilde{\lambda}_{+}^{2}\right)^{2}}},$
$\displaystyle\sin{2\phi_{\rm
SP}}=-\frac{\kappa}{\omega_{c}}\frac{\tilde{\lambda}_{-}^{2}-\lambda_{+}^{2}}{2\tilde{\lambda}_{-}\tilde{\lambda}_{+}}.$
(20)

The stability analysis is rather involved and we must resort to numerical
linearization around the steady-state. The final result is that SP↓ exists in
the purple and yellow regions of Fig. 1, where it is always stable. The phase
boundary can be divided into a second-order critical line (marked as a solid
curve in Fig. 1) and first-order transitions (dashed lines). The second-order
critical line can be parameterized by $r$ as
$\tilde{\lambda}_{-}=\tilde{\lambda}_{c}^{-}$ and
$\tilde{\lambda}_{+}=r\tilde{\lambda}_{c}^{-}$, where
$\tilde{\lambda}_{c}^{-}$ is given in Eq. (17). This curve ends at $r=r_{\pm}$
with two tricritical points [49, 50]. The one marked by an open circle in Fig.
1 is at $\tilde{\lambda}_{-}=\tilde{\lambda}_{c}^{\rm T}$ and
$\tilde{\lambda}_{+}=r_{-}\tilde{\lambda}_{c}^{\rm T}$, where $r_{-}$ is given
in Eq. (18) and:

$\tilde{\lambda}_{c}^{\rm
T}=\sqrt{2\left(1+\kappa^{2}/\omega_{c}^{2}+\sqrt{1+\kappa^{2}/\omega_{c}^{2}}\right)}.$
(21)

Beyond the tricritical points, the first-order phase boundaries (dashed lines
in Fig. 1) are simply given by
$\tilde{\lambda}_{+}=r_{\pm}\tilde{\lambda}_{-}$. Within the mean-field
approximation, when entering the normal phase across these boundaries the
order parameter $|c|^{2}$ of the superradiant state drops discontinuously to
zero. We also mention that another superradiant solution exists in the yellow
regions of Fig. 1. However, this fixed point is always unstable and does not
appear in the phase diagram. Furthermore, as stated already, Eq. (13) does not
give any chaotic or oscillatory phase.

In closing this section about the mean-field analysis, we expand upon our
earlier discussion about the regime of applicability. While mean-field theory
becomes accurate in the limit $\eta\to\infty$, we are interested here in the
stationary state of the system, i.e., we also impose $t\to\infty$. It is then
crucial to specify in which order the two limits are applied. From the
derivation of Eq. (13), it is clear that the phase diagram of Fig. 1 is
obtained by taking $\eta\to\infty$ before $t\to\infty$. More precisely, as we
will establish explicitly later on, the mean-field treatment is only valid on
a timescale $t_{\rm mf}$ which is proportional to $\eta$:

$t\ll t_{\rm mf}\propto\eta.$ (22)

If we take $\eta\to\infty$ from the onset, the above restriction does not play
any role. The mean-field theory predicts correctly that, in this limit, the
states appearing in Fig. 1 are stable for an arbitrarily long time.

Figure 2:  Mean-field phase diagram obtained by considering the unique steady
state of Eq. (3) in the limit $\eta\to\infty$. In the NP phase the system is
in a statistical mixture of NP↑ and NP↓, while the SP phase is formed by
statistical mixtures of NP↑ and SP↓. The dotted and solid curves mark first-
and second-order critical lines, respectively. We used the same parameters of
Fig. 2.

Alternatively, one might be interested in the properties of the steady state
at finite but arbitrarily large $\eta$. This is a different question which
assumes the opposite order of limits, i.e., $t\to\infty$ before
$\eta\to\infty$. Since the condition Eq. (22) is violated, the phase diagram
of Fig. 1 gets modified. In particular, multi-stable phases are not allowed in
this limit, since the quantum master equation has always a unique steady-
state. The phase diagram of the steady state at large $\eta$ (i.e.,
considering the regime $t\gg t_{\rm mf}$) is shown in Fig. 2, which summarizes
some of the main conclusions of the following sections. In this phase diagram
there are only a normal (NP) and superradiant (SP) phase. In each phase, the
system is in a well-defined statistical mixture of the mean-field states: NP↓
and NP↑ (SP↓ and NP↑) in the former (latter) case. Notably, the C phase of
Fig. 1 has disappeared, and the dotted line has become a first-order phase
boundary between the normal and superradiant phases.

##  IV Finite-$\eta$ effects

We now start our analysis of quantum fluctuations, which determine the fate of
the general multi-stable behavior found in Fig. 1. We first consider the
simpler yet instructive treatment of the normal phase. The cumulant expansion
is a rigorous approach to include quantum effects, and we show that it indeed
allows us to derive in a systematic manner the correct stationary properties
of the system (Sec. IV.1). However, the cumulant expansion is not physically
very transparent, and it is cumbersome to apply it to the superradiant state.
This is why we switch to an effective master-equation approach, recovering the
correct stationary properties of the normal phase (Sec. IV.2). We then
generalize this treatment to the superradiant phase (Sec. IV.3) and to the
tristable phase (Sec. IV.4), to establish the complete phase diagram shown in
Fig. 2.

###  IV.1 Cumulant expansion in the normal phase

In mean-field theory, we have used the factorization assumption
$\langle\hat{A}\hat{B}\rangle=\langle\hat{A}\rangle\langle\hat{B}\rangle$ to
decouple qubit-bosonic field correlations and obtained a set of dynamical
equations which only depend on $s_{z},s_{+}$ and the (properly rescaled)
bosonic field amplitude, $\langle\hat{a}\rangle$. Beyond the mean-field
approach, qubit-bosonic field correlations should be taken into account. For
example, from Eq. (II) we immediately obtain:

$\frac{d\langle\hat{\sigma}_{z}\rangle}{dt}=2i[\lambda_{-}(\langle\hat{a}\hat{\sigma}_{+}\rangle-\langle\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle)-\lambda_{+}(\langle\hat{a}\hat{\sigma}_{-}\rangle-\langle\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle)],$
(23)

where the evolution of the second-order correlators can be described through
the cumulant expansion [59, 60]. To leading-order, we assume that third-order
moments are zero, i.e., set
$\langle\hat{A}\hat{B}\hat{C}\rangle=\langle\hat{A}\hat{B}\rangle\langle\hat{C}\rangle+\langle\hat{A}\hat{C}\rangle\langle\hat{B}\rangle+\langle\hat{B}\hat{C}\rangle\langle\hat{A}\rangle-2\langle\hat{A}\rangle\langle\hat{B}\rangle\langle\hat{C}\rangle$.
By doing so, we obtain a set of non-linear equations such as:

$\displaystyle\frac{d\langle\hat{a}\hat{\sigma}_{+}\rangle}{dt}=$
$\displaystyle-(\kappa+i\omega_{c}-i\Omega)\langle\hat{a}\hat{\sigma}_{+}\rangle+\frac{i\lambda_{-}}{2}\left(1+\langle\hat{\sigma}_{z}\rangle\right)$
$\displaystyle+i(\lambda_{-}\langle\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{z}\rangle+\lambda_{+}\langle\hat{a}^{2}\hat{\sigma}_{z}\rangle),$
(24)

where the two last terms are approximated by:

$\displaystyle\langle\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{z}\rangle=$
$\displaystyle\langle\hat{a}^{\dagger}\hat{a}\rangle\langle\hat{\sigma}_{z}\rangle+\langle\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle\langle\hat{a}\rangle+\langle\hat{a}\hat{\sigma}_{z}\rangle\langle\hat{a}^{\dagger}\rangle$
$\displaystyle-2\langle\hat{a}\rangle\langle\hat{a}^{\dagger}\rangle\langle\hat{\sigma}_{z}\rangle,$
$\displaystyle\langle\hat{a}^{2}\hat{\sigma}_{z}\rangle=$
$\displaystyle\langle\hat{a}^{2}\rangle\langle\hat{\sigma}_{z}\rangle+2\langle\hat{a}\hat{\sigma}_{z}\rangle\langle\hat{a}\rangle-2\langle\hat{a}\rangle^{2}\langle\hat{\sigma}_{z}\rangle.$
(25)

The complete set of equations for the second-order correlators is presented in
Appendix A.

In the normal phase (i.e., the white region of Fig. 1) the order parameter
$\langle\hat{a}\rangle$ is zero, which allows us to simplify the non-linear
equations and obtain analytic expressions for the stationary expectation
values. The qubit population is given by:

$\displaystyle\langle\hat{\sigma}_{z}\rangle\simeq\frac{\lambda_{+}^{2}-\lambda_{-}^{2}}{\lambda_{-}^{2}+\lambda_{+}^{2}}+\Delta,$
(26)

where $\Delta$ is a small ($\sim 1/\eta$) correction, whose explicit
expression is given in Appendix A. Higher-order terms are omitted, as they are
not significant (to derive them, it would be necessary to include higher-order
correlators in the cumulant expansion).

The leading term of Eq. (26) already captures the effects of quantum
fluctuations, which are beyond mean-field theory. In the normal phase, mean-
field theory predicts that both NP↑ and NP↓ are stable fixed points with
distinct basins of attraction, thus the stationary expectation value at
$t\to\infty$ is $s_{z}=\pm 1$, depending on the initial condition. Instead,
Eq. (26) gives a very different result. For example, in the balanced case
$\lambda_{-}=\lambda_{+}$ we obtain a small expectation value
$\langle\hat{\sigma}_{z}\rangle\simeq-2/\eta$. This difference is due to
fluctuations between the two mean-field solutions, induced by the quantum
evolution. At any finite $\eta\gg 1$, the exact dynamics described by Eq. (3)
leads to a stationary state which is a well-defined statistical mixture of
$\rm{NP}_{\uparrow}$ and $\rm{NP}_{\downarrow}$. From Eq. (26) we can infer
that, to leading order in $\eta$, the populations of these two states are
respectively given by $P_{+}$ and $P_{-}$, where:

$P_{\pm}=\frac{\lambda_{\pm}^{2}}{\lambda_{+}^{2}+\lambda_{-}^{2}}.\quad{\rm(normal~{}phase)}$
(27)

Figure 3:  Time evolution of $\langle\hat{\sigma}_{z}\rangle$ for the
isotropic Rabi model ($\tilde{\lambda}_{+}=\tilde{\lambda}_{-}$) at a finite
value $\eta=100$. The purple squares are calculated with the master equation
(3), considering both $|\uparrow\rangle$ and $|\downarrow\rangle$ as initial
states (the bosonic mode is initially in the vacuum). The green dashed-dotted
lines are computed numerically from the cumulant expansion, see Eq. (IV.1) and
Appendix A, and show excellent agreement with the exact dynamics. The red line
is the approximate stationary value of Eq. (26). The blue dotted lines are
calculated based on the mean-field equations (8–10). In the main panel, the
mean-field evolution does not display any significant time dependence. For
completeness, we also show the mean-field evolution for much longer times
(inset). We used the parameters $\kappa=0.5\omega_{c}$ and
$\tilde{\lambda}_{+}=\tilde{\lambda}_{-}=0.618$.

To verify the validity of the cumulant expansion, we have performed an exact
calculation based on Eq. (3), see the purple squares in Fig. 3. In this
numerical example we consider the isotropic model ($\lambda_{-}=\lambda_{+}$)
and two initial conditions, corresponding to either NP↑ or NP↓. As seen from
Fig. 3, mean-field theory (dotted curves) predicts that the two initial values
$\langle\hat{\sigma}_{z}\rangle=\pm 1$ are constant on a timescale which is
much longer than the true relaxation time. Instead, the numerical solutions
obtained from the cumulant expansion (green dashed-dotted curves) are in
excellent agreement with the master equations. The stationary value of
$\langle\hat{\sigma}_{z}\rangle$ is also very close to the one predicted by
Eq. (26).

Figure 3 allows us to discuss more concretely the previous arguments about the
$\eta\to\infty$ and $t\to\infty$ limits. The semiclassical treatment becomes
formally exact when $\eta\to\infty$, but Fig. 3 shows that this is only true
if the evolution is shorter than the typical relaxation time associated with
the quantum evolution. As seen, mean-field theory fails over the long
evolution time considered in Fig. 3 (where $\tau\gtrsim 10^{3}\gg\eta$). If in
Fig. 3 we restrict ourselves to shorter times (e.g., $\tau\ll\eta$), the
discrepancy between exact evolution and semiclassical dynamics remains small.
It is also important to note that the relaxation time of the exact quantum
dynamics grows with $\eta$, so mean-field theory becomes accurate for longer
time periods when taking the $\eta\to\infty$ limit.

The cumulant expansion can be also applied to the other multi-stable phases of
Fig. 1. However, in this case we cannot easily obtain analytic results. For
this reason, we prefer the alternative approach described in the following
sections. Working to leading order in $\eta$, we will derive effective master
equations which describe in a transparent manner the transitions between
different fixed points. In the normal phase, we can recover the first term of
Eq. (26) and the treatment can be extended to the rest of the phase diagram.

###  IV.2 Effective master equation in the normal phase

The derivation of the effective master equation in the normal phase is
relatively straightforward. As a first step, we apply a unitary transformation
$\hat{U}_{np}$, which approximately decouples the qubit and cavity field [29,
31]:

$\displaystyle\hat{U}_{np}=e^{\frac{1}{2\sqrt{\eta}}\left[\left(\tilde{\lambda}_{-}\hat{a}+\tilde{\lambda}_{+}\hat{a}^{\dagger}\right)\hat{\sigma}_{+}-\left(\tilde{\lambda}_{+}\hat{a}+\tilde{\lambda}_{-}\hat{a}^{\dagger}\right)\hat{\sigma}_{-}\right]}.$
(28)

By neglecting terms of order $\omega_{c}/\sqrt{\eta}$ in the effective
Hamiltonian, $H_{np}=\hat{U}_{np}^{\dagger}\hat{H}\hat{U}_{np}$, we obtain:

$\displaystyle\hat{H}_{np}$
$\displaystyle\simeq\,\frac{\Omega}{2}\left(1+\frac{\tilde{\lambda}_{-}^{2}+\tilde{\lambda}_{+}^{2}}{4\eta}\right)\hat{\sigma}_{z}+\omega_{c}\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\frac{\omega_{c}}{4}\left[\left(\tilde{\lambda}_{-}^{2}+\tilde{\lambda}_{+}^{2}\right)\hat{a}^{\dagger}\hat{a}+\tilde{\lambda}_{-}\tilde{\lambda}_{+}\left(\hat{a}^{2}+\hat{a}^{\dagger
2}\right)\right]\hat{\sigma}_{z}.$ (29)

Since $\hat{H}_{np}$ is diagonal in $\hat{\sigma}_{z}$ and, to leading order
in $\eta$, only cavity decay takes place, the $\hat{\sigma}_{z}=\pm 1$
subspaces are decoupled also for the dissipative time evolution. This leading-
order approximation is consistent with the mean-field treatment, when there
are two stable fixed points NP↑ and NP↓ with no transition taking place
between them. However, to describe the evolution of
$\hat{\widetilde{\rho}}=\hat{U}_{np}^{\dagger}\hat{\rho}\hat{U}_{np}$ more
accurately, it is necessary to consider how the Lindblad term is affected by
the unitary transformation. The hybridization with the cavity field has the
well-known consequence of inducing decay of the two-level atom (i.e., Purcell
relaxation). For completeness, the transformation of the Lindlad term is
discussed in detail in Appendix B, while we include below the most important
effect:

$\displaystyle\dot{\hat{\widetilde{\rho}}}=$
$\displaystyle-i\left[\hat{H}_{np},\hat{\widetilde{\rho}}\right]+\kappa\mathcal{D}[\hat{a}]\hat{\widetilde{\rho}}+\frac{\kappa}{4\eta}\sum_{\mu=\pm}\tilde{\lambda}^{2}_{\mu}\mathcal{D}[\hat{\sigma}_{\mu}]\hat{\widetilde{\rho}}.$
(30)

The last term of Eq. (30) describes spin flips beyond the mean-field
approximation. As mentioned in Eq. (22), these spin-flips occur on a long
timescale $\sim\eta/\kappa$, but are essential to establish the unique
stationary state at any finite $\eta$. Since the interaction contains co- and
counter-rotating terms, both relaxation and excitation processes appear. The
corresponding rates:

$\displaystyle\Gamma_{{\rm NP}_{\uparrow}\to{\rm
NP}_{\downarrow}}=\frac{\tilde{\lambda}_{-}^{2}}{2\eta}\kappa,$ (31)
$\displaystyle\Gamma_{{\rm NP}_{\downarrow}\to{\rm
NP}_{\uparrow}}=\frac{\tilde{\lambda}_{+}^{2}}{2\eta}\kappa,$ (32)

give the populations in Eq. (27) and the leading term of Eq. (26). By
interpreting the rates in terms of an effective temperature $T_{\rm eff}$, we
immediately have (with $k_{\rm B}=1$):

$T_{\rm eff}=\frac{\Omega}{2\ln(\lambda_{-}/\lambda_{+})}.$ (33)

Interestingly, the incoherent pumping induced by the counter-rotating terms
leads to an infinite-temperature spin state in the isotropic model
($\lambda_{+}=\lambda_{-}$), while population inversion is realized when
$\lambda_{+}>\lambda_{-}$.

###  IV.3 Second-order NP-SP transition

After having discussed the normal phase, we now consider the superradiant
states. As we will see, the discussion of transition rates becomes more subtle
than in the normal phase. In this subsection we characterize the stationary
state in the purple region of Fig. 1, where the stable fixed points of mean-
field theory are NP↑ and SP↓. The coexistence region C, where also NP↓ is
stable, will be discussed in the next subsection.

To derive an effective master equation which is valid close to the
superradiant fixed points, we should take into account the finite expectation
value of the cavity field. We consider one of the two mean-field solutions for
$c$, while the other broken-symmetry state can be treated in a completely
analogous manner. It is convenient to first perform a displacement
transformation
$\hat{U}_{d}=e^{\sqrt{\eta}\left(c\hat{a}^{\dagger}-c^{*}\hat{a}\right)}$
[29], giving the effective Hamiltonian
$\hat{H}_{d}=\hat{U}_{d}^{\dagger}\hat{H}\hat{U}_{d}+i\sqrt{\eta}\kappa(c^{*}\hat{a}-c\hat{a}^{\dagger})$.
Furthermore, since the natural spin quantization axis is tilted away from the
$z$ direction, we define:

$\displaystyle\vec{\hat{\tau}}=Y(\theta_{\rm SP})Z(-\phi_{\rm
SP})\vec{\hat{\sigma}},$ (34)

where $Z(\pm\phi_{\rm SP})$ and $Y(\theta_{\rm SP})$ are rotation matrices
around the $z$ and $y$ axis, respectively, with $\theta_{\rm SP},\phi_{\rm
SP}$ given by Eq. (III). After these preliminary transformations, we obtain:

$\displaystyle\hat{H}_{d}=\,$
$\displaystyle\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2\cos\theta_{\rm
SP}}\hat{\tau}_{z}-\frac{\sqrt{\eta}\omega_{c}}{2}(1+\hat{\tau}_{z})\left(\tilde{\chi}_{z}\hat{a}+\tilde{\chi}^{*}_{z}\hat{a}^{\dagger}\right)$
$\displaystyle~{}-\frac{\sqrt{\eta}\omega_{c}}{2}\left[\tilde{\chi}_{-}\hat{\tau}_{+}\hat{a}+\tilde{\chi}_{+}\hat{\tau}_{+}\hat{a}^{\dagger}+{\rm
H.c.}\right],$ (35)

where the dimensionless couplings are:

$\displaystyle\tilde{\chi}_{\pm}=\tilde{\lambda}_{\pm}e^{i\phi_{\rm
SP}}\cos^{2}\frac{\theta_{\rm SP}}{2}+\tilde{\lambda}_{\mp}e^{-i\phi_{\rm
SP}}\sin^{2}\frac{\theta_{\rm SP}}{2},$
$\displaystyle\tilde{\chi}_{z}=-\frac{1}{2}\sin\theta_{\rm
SP}\left(\tilde{\lambda}_{-}e^{i\phi_{\rm
SP}}+\tilde{\lambda}_{+}e^{-i\phi_{\rm SP}}\right).$ (36)

The second line of Eq. (IV.3) contains co- and counter-rotating interactions
with modified strength, given by $\tilde{\chi}_{\pm}$. To leading order, we
can eliminate these off-diagonal terms by applying a unitary transformation
$\hat{U}_{sp}$ analogous to Eq. (28). Details are given in Appendix B, where
we derive the following master equation:

$\displaystyle\dot{\hat{\tilde{\rho}}}=$
$\displaystyle-i\left[\hat{H}_{sp},\hat{\bar{\rho}}\right]+\kappa\mathcal{D}[\hat{a}]\hat{\bar{\rho}}$
$\displaystyle+\frac{\kappa}{4\eta}\cos^{2}{\theta_{\rm
SP}}\left(|\tilde{\chi}_{-}|^{2}\mathcal{D}[\hat{\tau}_{-}]\hat{\bar{\rho}}+|\tilde{\chi}_{+}|^{2}\mathcal{D}[\hat{\tau}_{+}]\hat{\bar{\rho}}\right).$
(37)

In the above Eq. (IV.3), the effective Hamiltonian reads:

$\displaystyle\hat{H}_{sp}\simeq\frac{\Omega}{2}\left(\frac{1}{\cos\theta_{\rm
SP}}+\frac{|\tilde{\chi}_{-}|^{2}+|\tilde{\chi}_{+}|^{2}}{4\eta}\cos\theta_{\rm
SP}\right)\hat{\tau}_{z}+\omega_{c}\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\frac{\omega_{c}}{4}\left[(|\tilde{\chi}_{-}|^{2}+|\tilde{\chi}_{+}|^{2})\hat{a}^{\dagger}\hat{a}+(\tilde{\chi}_{-}\tilde{\chi}_{+}^{*}\hat{a}^{2}+\tilde{\chi}_{-}^{*}\tilde{\chi}_{+}\hat{a}^{\dagger
2})\right]\hat{\tau}_{z},$
$\displaystyle-\frac{\sqrt{\eta}\omega_{c}}{2}(\tilde{\chi}_{z}\hat{a}+\tilde{\chi}_{z}^{*}\hat{a}^{{\dagger}})(1+\hat{\tau}_{z}).$
(38)

Despite the more complex transformation (here the density matrix is
$\hat{\bar{\rho}}=\hat{U}_{sp}^{\dagger}\hat{U}_{d}^{\dagger}\hat{\rho}\hat{U}_{d}\hat{U}_{sp}$),
the type of spin-flip processes appearing in the second line of Eq. (IV.3) are
completely analogous to the last term of Eq. (30). In the superradiant case,
$\tilde{\chi}_{\pm}$ and $\hat{\tau}_{\alpha}$ play the role of
$\tilde{\lambda}_{\pm}$ and $\hat{\sigma}_{\alpha}$, respectively. The spin-
flip rates are now given by:

$\displaystyle\Gamma_{{\rm SP}_{\uparrow}\to{\rm
SP}_{\downarrow}}=\frac{|\tilde{\chi}_{-}|^{2}}{2\eta}\kappa\cos^{2}\theta_{\rm
SP},$ (39) $\displaystyle\Gamma_{{\rm SP}_{\downarrow}\to{\rm
SP}_{\uparrow}}=\frac{|\tilde{\chi}_{+}|^{2}}{2\eta}\kappa\cos^{2}\theta_{\rm
SP},$ (40)

where the additional factors $\cos^{2}\theta_{\rm SP}$ are related to the fact
that the $\tau_{z}=\pm 1$ eigenstates have spin polarization
$\langle\hat{\sigma}_{z}\rangle=\pm\cos\theta_{\rm SP}$. On the other hand, if
we consider the Hamiltonian in Eq. (IV.3), we find an important difference
from the normal phase. The first three lines have the same structure of the
normal-phase Hamiltonian $\hat{H}_{np}$, see Eq. (IV.2), but the last term of
$H_{sp}$ has no analogy to the treatment of the normal state and leads to
important physical consequences described below.

The third line of Eq. (IV.3) is a spin-dependent drive of the cavity, which
vanishes when $\hat{\tau}_{z}=-1$. Therefore, if we restrict ourselves to the
subspace $\hat{\tau}_{z}=-1$, the drive term is absent and the system will
relax to a stationary state around $\langle\hat{a}^{\dagger}\hat{a}\rangle=0$
(in the displaced frame). Instead, in the subspace $\hat{\tau}_{z}=1$ the
drive term is large and leads to a rapid increase of the average photon
number, bringing the system beyond the regime of applicability of Eqs. (IV.3)
and (IV.3). The interpretation of this behavior is straightforward from the
point of view of the mean-field treatment: SP↓ is a stable fixed point, thus
the evolution with $\hat{\tau}_{z}=-1$ leads to a well-defined stationary
state. However, the state SP↑ (obtained from SP↓ by flipping the spin
direction) is not a stationary mean-field solution. SP↑ undergoes a fast
evolution which brings the photon field away from the expectation value
$\langle\hat{a}\rangle\simeq c\sqrt{\eta}$ (in the original frame), and
eventually leads to the other stable fixed point, NP↑.

Figure 4:  Relevant transitions taking place in the superradiant phase. The
filled dots (NP↑ and SP↓) are stable fixed points within the mean-field
approximation, while the empty dots (NP↓ and SP↑) are unstable states. Fast
transitions are indicated by thick arrows. The gray dashed arrows indicate
slow spin flips described by the effective master equations. The corresponding
rates (all $\propto\kappa/\eta$) are given in Eqs. (31), (32), (39), and (40).

We now take into account all slow spin-flip processes, by including both
normal and superradiant states. For convenience, the various transitions
discussed are schematically represented in Fig. 4. First, we consider SP↓ as
the initial state. On a long time scale corresponding to Eq. (40) a spin flip
takes place, bringing the system to SP↑. Afterward, SP↑ evolves to NP↑ within
a short timescale. As described in the previous paragraph, this
SP${}_{\uparrow}\to$ NP↑ process is governed by the mean-field evolution, thus
the corresponding transition rate is not suppressed by $\eta$. Since ${\rm
SP}_{\uparrow}\to{\rm NP}_{\uparrow}$ can be considered as essentially
instantaneous, the small spin-flip rate $\Gamma_{{\rm SP}_{\uparrow}\to{\rm
SP}_{\downarrow}}$ in Eq. (39) does not play any significant role and we
obtain a dominant chain of transitions SP${}_{\downarrow}\to$
SP${}_{\uparrow}\to$ NP↑, with the effective rate given by Eq. (40):

$\Gamma_{{\rm SP}_{\downarrow}\to{\rm
NP}_{\uparrow}}\simeq\frac{|\tilde{\chi}_{+}|^{2}}{2\eta}\kappa\cos^{2}\theta_{\rm
SP}.$ (41)

Similarly, if the initial state is NP↑ a slow spin-flip process to NP↓ takes
place, with a rate given by Eq. (31). Since NP↓ is unstable in the
superradiant phase, the system quickly reaches ${\rm SP}_{\downarrow}$. The
effective rate for the whole process, ${\rm NP}_{\uparrow}\to{\rm
NP}_{\downarrow}\to{\rm SP}_{\downarrow}$, is given by:

$\Gamma_{{\rm NP}_{\uparrow}\to{\rm
SP}_{\downarrow}}\simeq\frac{\tilde{\lambda}_{-}^{2}}{2\eta}\kappa.$ (42)

From Eqs. (41) and (42) we find the probability of being in the superradiant
state:

$P_{{\rm
SP}_{\downarrow}}=\frac{\tilde{\lambda}_{-}^{2}}{\tilde{\lambda}_{-}^{2}+|\tilde{\chi}_{+}|^{2}\cos^{2}\theta_{\rm
SP}},$ (43)

while the population of NP↑ is simply given by $P_{{\rm
NP}_{\uparrow}}=1-P_{{\rm SP}_{\downarrow}}$. At the second-order phase
boundary (i.e., the solid curve in Fig. 1), these populations join smoothly
with the $P_{\pm}$ of Eq. (27), since $\theta_{\rm SP}\to 0$ and
$|\tilde{\chi}_{+}|^{2}\to\tilde{\lambda}_{+}^{2}$.

Figure 5:  Dependence of $\langle\hat{a}^{\dagger}\hat{a}\rangle/\eta$ on the
coupling strength $\tilde{\lambda}_{-}$ across the second-order NP-SP
transition (the solid critical line of Fig. 2). Squares are exact numerical
data, obtained by solving Eq. (3). The solid curves are the analytical
treatment, see Eq. (44). The dashed curves are the order parameter $|c|^{2}$
of the pure ${\rm SP}_{\downarrow}$ state, see Eq. (III). In panel (a) we have
used $\tilde{\lambda}_{+}/\tilde{\lambda}_{-}=1$ and $\eta=400$ while in panel
(b) $\tilde{\lambda}_{+}/\tilde{\lambda}_{-}=3$ and $\eta=500$. Other
parameters are the same as Fig. 1.

The value of $P_{{\rm SP}_{\downarrow}}$ has a direct effect on the order
parameter of the second-order phase transition. The average photon number in
the superradiant phase is:

$\langle\hat{a}^{\dagger}\hat{a}\rangle=P_{{\rm SP}_{\downarrow}}|c|^{2}\eta,$
(44)

where $c$, given by Eq. (III), is for the pure ${\rm SP}_{\downarrow}$ state.
In Fig. 5 we show a comparison of Eq. (44) to the numerical results obtained
by solving the master equation at large $\eta$, finding excellent agreement.
The admixture with ${\rm NP}_{\uparrow}$ can cause a significant reduction
from the photon number of the pure superradiant state (dashed curves). The
discrepancy is more obvious close to the phase transition, where $\theta_{\rm
SP}\simeq 0$. Instead, deep in the superradiant phase the spin direction moves
away from the $z$ axis, which causes a reduction of $\cos\theta_{\rm SP}$ and
a growth of $P_{{\rm SP}_{\downarrow}}$ [see Eq. (43)]. As a result, we see in
Fig. 5(a) that the exact photon number tends to approach the dashed curve by
increasing $\tilde{\lambda}_{-}$. Furthermore, close to the critical point we
have $|\tilde{\chi}_{+}|\simeq\tilde{\lambda}_{+}$, thus we can approximate
$P_{{\rm SP}_{\downarrow}}\simeq(1+\lambda^{2}_{+}/\lambda^{2}_{-})^{-1}$.
This expression indicates that a larger value of $\lambda_{+}/\lambda_{-}$
leads to a more dramatic reduction of
$\langle\hat{a}^{\dagger}\hat{a}\rangle$, in agreement with panel (b). The
physical interpretation is quite clear: Strong counter-rotating terms can
excite simultaneously the qubit and the cavity. This, in conjunction with
cavity decay, drives a large population to the NP↑ high-energy state,
suppressing superradiance.

###  IV.4 First-order phase transition

We now examine the coexistence region C of Fig. 1, where the stationary state
is expected to be a statistical mixture of the ${\rm NP}_{\uparrow}$, ${\rm
NP}_{\downarrow}$, and ${\rm SP}_{\downarrow}$. Remarkably, instead, we find
that in region C the system returns to the same normal state described before,
i.e., the populations of ${\rm NP}_{\uparrow/\downarrow}$ are as in Eq. (27)
and $P_{{\rm SP}_{\downarrow}}=0$. Therefore, the phase diagram for the
stationary states (i.e., taking $t\to\infty$ before $\eta\to\infty$) becomes
as in Fig. 2, where the mean-field phase boundary between the superradiant
(SP) and C regions has now become a first-order transition between SP and NP.

Figure 6:  Dependence of the photon number on the coupling strength. In panel
(a) we plot $\langle\hat{a}^{\dagger}\hat{a}\rangle/\eta$ as a function of
$\tilde{\lambda}_{-}$. The three curves are for $\eta=200$ (solid), 400
(dashed), and 600 (dot-dashed). In panels (b-e) we show the detailed photon
distribution with $\eta=200$ and four representative values of
$\tilde{\lambda}_{-}$, marked by symbols in panel (a). Panel (b) and (c) are
respectively in the normal and superradiant phases. Panels (d) and (e) are
both in the ‘coexistence’ phase C of Fig. 1. However, a significant mixture of
normal and superradiant states only occurs for panel (d). At large
$\tilde{\lambda}_{-}$ the system returns to the normal state, see panel (e).
The vertical dashed lines in panels (c-e) indicate the mean-field value of
$\langle\hat{a}^{\dagger}\hat{a}\rangle$ for the superradiant state, given by
Eq. (III). In all panels $\tilde{\lambda}_{+}=0.3\tilde{\lambda}_{-}$ and
$\kappa=0.5\omega_{0}$.

In Fig. 6, we show evidence of this behavior from numerical simulations based
on the master equation. We consider relatively large values of $\eta$ and
compute the average photon number as a function of $\tilde{\lambda}_{-}$ (at a
fixed ratio $\tilde{\lambda}_{+}/\tilde{\lambda}_{-}=0.3$), finding that the
monotonic increase of $\langle\hat{a}^{\dagger}\hat{a}\rangle$ breaks down
after entering the C region. As seen in Fig. 6(a), there is a sudden drop at
sufficiently large $\tilde{\lambda}_{-}$, bringing the system close to
$\langle\hat{a}^{\dagger}\hat{a}\rangle/\eta\approx 0$. By increasing $\eta$,
the drop occurs at progressively smaller values of $\tilde{\lambda}_{-}$ and
becomes sharper, suggesting the occurrence of a first-order transition when
$\eta\to\infty$.

The various regimes can also be studied through the photon number distribution
in Fock space, shown in panels (b-e) of Fig. 6 at four representative values
of $\tilde{\lambda}_{-}$. The behavior in panels (b) and (c) is as expected:
in the former (latter) case a prominent peak around
$\langle\hat{a}^{\dagger}\hat{a}\rangle=0$
($\langle\hat{a}^{\dagger}\hat{a}\rangle=\eta|c|^{2}$) is found, corresponding
to the NP (SP↓) mean-field solution. Note that in panel (c), which should
generically show a mixture of NP↑ and SP↓ states, there is no visible peak
around zero photons. The small population of ${\rm NP}_{\uparrow}$ is due to
large value of $\tilde{\lambda}_{-}$, which is far from the second-order
critical point, and the small ratio
$\tilde{\lambda}_{+}/\tilde{\lambda}_{-}=0.3$ (see the discussions at the end
of Sec. IV.3). In panel (d) of Fig. 6 we show a state which is clearly a
nontrivial mixture of superradiant and normal states. The photon distribution
has two prominent peaks, but this behavior is only found in the narrow
transition region where $\langle\hat{a}^{\dagger}\hat{a}\rangle$ drops. By
increasing $\tilde{\lambda}_{-}$, the superradiant peak progressively
disappears and the system eventually returns to the normal state, see panel
(e).

Figure 7:  Dependence of
$\Delta\tilde{\lambda}_{-}=\tilde{\lambda}_{-}^{*}-\tilde{\lambda}_{c}^{+}$ on
$\eta$. Here $\tilde{\lambda}_{-}^{*}$ is the position at which there is a
rapid drop in photon number, see Fig. 6(a). It is found from the inflection
points, i.e., by imposing
$d^{2}\langle\hat{a}^{\dagger}\hat{a}\rangle/d\tilde{\lambda}_{-}^{2}=0$. The
solid line shows a power-law dependence of type $1/\sqrt{\eta}$. The
parameters are as in Fig. 6.

Unfortunately, the limit $\eta\to\infty$ is difficult to approach through
exact numerical simulations. Due to the large photon number in the SP phase
(where $\langle\hat{a}^{\dagger}\hat{a}\rangle\propto\eta$), computing
accurately the drop of $\langle\hat{a}^{\dagger}\hat{a}\rangle$ requires a
large Hilbert space. This point can be also appreciated from panels (b-e) of
Fig. 6. At the largest value of $\eta$ shown in Fig. 6(a), convergence to the
$\eta\to\infty$ limit has not yet been achieved. To infer the position of the
critical point we plot in Fig. 7 the dependence on $\eta$ of
$\Delta\tilde{\lambda}_{-}=\tilde{\lambda}_{-}^{*}-\tilde{\lambda}_{c}^{+}$,
where $\tilde{\lambda}_{-}^{*}$ marks the inflection point of
$\langle\hat{a}^{\dagger}\hat{a}\rangle$ and $\tilde{\lambda}_{c}^{+}$ is the
SP/C boundary, given in Eq. (17). We find a dependence of the type
$\Delta\tilde{\lambda}_{-}\sim\eta^{-1/2}$, indicating that the critical point
is indeed at $\tilde{\lambda}_{c}^{+}$.

Figure 8:  Relevant transitions taking place in the C region. The notation is
the same of Fig. 4 but now also ${\rm NP}_{\downarrow}$ is a stable fixed
point, with its own basin of attraction.

The above numerical analysis is consistent with a physical argument based on
effective master equations. Similarly to Fig. 4, we represent in Fig. 8 the
transitions between mean-field solutions in the C region. As discussed before,
the slow process ${\rm SP}_{\downarrow}\to{\rm SP}_{\uparrow}$ allows the
transition pathway ${\rm SP}_{\downarrow}\to{\rm SP}_{\uparrow}\to{\rm
NP}_{\uparrow}$ and, like in the superradiant phase, the metastable state
${\rm NP}_{\uparrow}$ undergoes the slow transition ${\rm
NP}_{\uparrow}\to{\rm NP}_{\downarrow}$. However, differently from Fig. 4,
${\rm NP}_{\downarrow}$ is now a stable fixed point, thus the process ${\rm
NP}_{\downarrow}\to{\rm SP}_{\downarrow}$ is not available. Once the system
has reached one of the two normal states, it will perform ${\rm
NP}_{\uparrow}\leftrightarrow{\rm NP}_{\downarrow}$ transitions with the same
rates of Eqs. (31) and (32) and will never return to the superradiant state.
Therefore, we find the same statistical mixture of the normal phase, with no
discontinuity at the dashed boundaries of Fig. 1. The main predictions of Fig.
8, i.e., the occurrence of a first-order phase transition at
$\tilde{\lambda}_{-}=\tilde{\lambda}_{c}^{+}$ (dotted lines in Fig. 1) and the
absence of a superradiant contribution in the C region, are in agreement with
our numerical findings.

##  V Scaling analysis

In this section we return to the second-order phase transition between the NP
and SP and perform a finite-$\eta$ scaling analysis based on an effective
equilibrium theory. This approach was originally proposed for the dissipative
Dicke model [55], where it shows good agreement with the more rigorous Keldysh
path-integral approach [55, 61, 62]. Here, since the phase transition is due
to the $\downarrow$ subspace, we initially assume $\hat{\sigma}_{z}=-1$. We
need, however, to go beyond the quadratic approximation of the effective
Hamiltonian adopted so far and include the leading-order nonlinear terms. With
the quadratures $\hat{x}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ and
$\hat{p}=i(\hat{a}^{\dagger}-\hat{a})/\sqrt{2}$, we find (see Appendix B):

$\displaystyle\hat{H}_{eff}\simeq\omega_{c}\left[\frac{\alpha_{y}}{2}\hat{p}^{2}+\frac{\alpha_{x}}{2}\hat{x}^{2}+\frac{(\tilde{\lambda}_{x}^{2}\hat{x}^{2}+\tilde{\lambda}_{y}^{2}\hat{p}^{2})^{2}}{16\eta}\right],$
(45)

where the rescaled coupling constants
$\tilde{\lambda}_{x,y}=(\tilde{\lambda}_{+}\pm\tilde{\lambda}_{-})/2$
correspond to Eq. (2). Furthermore:

$\alpha_{x,y}=1-\tilde{\lambda}_{x,y}^{2}.$ (46)

The quadratic part of Eq. (45) is in agreement with $\hat{H}_{np}$ of Eq.
(IV.2) and the additional terms are suppressed by a factor of $1/\eta$, thus
they can be neglected in the normal phase. These quartic terms, however,
become crucial in the critical region since
$\langle\hat{x}^{2}\rangle\sim\eta$ in the superradiant phase. At variance
with the closed system [31], both $\langle\hat{x}^{2}\rangle$ and
$\langle\hat{p}^{2}\rangle$ scale linearly with $\eta$ in the superradiant
phase. Therefore, we cannot neglect $\tilde{\lambda}_{y}^{2}p^{2}$ in the
quartic term.

To compute the critical properties, we consider the following equation of
motion:

$\frac{d^{2}\hat{x}}{dt^{2}}=-\omega_{c}^{2}\left(\alpha^{2}\hat{x}+\frac{\beta^{3}}{\eta}\hat{x}^{3}\right)-2\kappa\frac{d\hat{x}}{dt}+\hat{\mathcal{G}},$
(47)

which can be derived from the Heisenberg-Langevin equations associated to Eq.
(45). Details are presented in Appendix C, where we obtain:

$\displaystyle\alpha^{2}=\alpha_{x}\alpha_{y}+\frac{\kappa^{2}}{\omega_{c}^{2}},\quad\beta^{3}=\alpha_{y}\tilde{\lambda}^{4}_{x}-\frac{1}{\alpha_{y}^{3}}\left(\frac{\kappa\tilde{\lambda}_{y}}{\omega_{c}}\right)^{4},$
(48)

and the expression of the noise operator $\hat{\mathcal{G}}(t)$, given in Eq.
(82). While Eq. (47) is an operator equation, we follow the idea that the
critical behavior is determined by classical fluctuations [55]. Then, an
accurate description of the critical behavior can be obtained from the
classical limit of Eq. (47), describing a Brownian particle. The potential is:

$V(x)=\omega_{c}\left(\frac{1}{2}\alpha^{2}x^{2}+\frac{1}{4\eta}\beta^{3}x^{4}\right),$
(49)

which is in agreement with the mean-field treatment, since the second-order
phase boundary is recovered by setting $\alpha^{2}=0$ and the equilibrium
position $x=-\alpha^{2}/\beta^{3}$ is consistent with the mean-field value of
$c$ in Eq. (III). In Appendix C we also derive the effective temperature
associated with the noise operator $\hat{\mathcal{G}}$:

$T=\frac{\kappa^{2}+\alpha_{y}^{2}\omega_{c}^{2}}{4\omega_{c}}.$ (50)

To compute the critical properties, we should finally recall that the above
classical limit only refers to the $\hat{\sigma}_{z}=-1$ subspace. The
fluctuations in the $\hat{\sigma}_{z}=+1$ subspace do not diverge with $\eta$,
so the contribution of these states is negligible. This implies that a
prefactor $P_{-}$ (i.e., the population of $\hat{\sigma}_{z}=-1$) should be
included when computing the expectation values of $\langle\hat{x}^{2}\rangle$
by a classical average:

$\displaystyle\langle\hat{x}^{2}\rangle=$ $\displaystyle P_{-}\frac{\int
x^{2}e^{-V(x)/T}dx}{\int e^{-V(x)/T}dx}$ $\displaystyle=2P_{-}\sqrt{\frac{\eta
T}{\beta^{3}\omega_{c}}}Q\left(-\frac{\alpha^{2}}{2}\sqrt{\frac{\eta\omega_{c}}{\beta^{3}T}}\right).$
(51)

The scaling form in the second line of Eq. (V) is obtained after a simple
change of variable in the integrals, giving:

$Q(x)=\frac{\int z^{2}e^{-(x-z^{2})^{2}}dz}{\int e^{-(x-z^{2})^{2}}dz}.$ (52)

It is also useful to linearize $\alpha^{2}$ close to the critical point. By
considering the line $\tilde{\lambda}_{+}=r\tilde{\lambda}_{-}$, we expand
$\alpha^{2}$ in terms of
$\delta\tilde{\lambda}_{-}=\tilde{\lambda}_{-}-\tilde{\lambda}_{c}^{-}$, where
the critical point $\lambda_{c}^{-}$ is given by Eq. (17) and:

$\displaystyle\alpha^{2}\simeq-\tilde{\lambda}_{c}^{-}\left[1+r^{2}-\frac{(1-r^{2})^{2}}{4}(\tilde{\lambda}_{c}^{-})^{2}\right]\delta\tilde{\lambda}_{-}\equiv-\Lambda\delta\tilde{\lambda}_{-}.$
(53)

From Eq. (V), we can also easily obtain the expression for the average photon
number $\langle\hat{a}^{\dagger}\hat{a}\rangle$, using the approximate
relation $\hat{p}\simeq\kappa\hat{x}/\alpha_{y}\omega_{c}$ (see Appendix C).

###  V.1 Critical fluctuations

We now evaluate Eq. (V) on the second-order phase boundary, which can be
parametrized with $r=\tilde{\lambda}_{+}/\tilde{\lambda}_{-}$. The
fluctuations at the critical point are:

$\langle\hat{x}^{2}\rangle=\left(\frac{2Q(0)}{1+r^{2}}\sqrt{\frac{T}{\beta^{3}\omega_{c}}}\right)\sqrt{\eta},$
(54)

where $Q(0)=\Gamma(3/4)/\Gamma(1/4)\simeq 0.338$ and the values of $\beta^{3}$
and $T$ can be computed using:

$\tilde{\lambda}_{x}=\frac{r+1}{2}\tilde{\lambda}_{c}^{-},\qquad\tilde{\lambda}_{y}=\frac{r-1}{2}\tilde{\lambda}_{c}^{-},$
(55)

where $\tilde{\lambda}_{c}^{-}$ is given in Eq. (17). Since both $T$ and
$\beta^{3}$ are evaluated using Eqs. (55) and (17), the prefactor in Eq. (54)
becomes a function of $r$ and $\kappa/\omega_{c}$. For the isotropic Rabi
model ($r=1$) we find the simpler formula:

$\langle\hat{x}^{2}\rangle=\frac{Q(0)}{2}\sqrt{\frac{\eta}{1+\kappa^{2}/\omega_{c}^{2}}}.\qquad(r=1)$
(56)

In Fig. 9, we show a comparison of the approximate Eq. (54) to exact numerical
results at three different points on the phase boundary, finding excellent
agreement at sufficiently large $\eta$.

Figure 9:  Dependence of $\langle\hat{x}^{2}\rangle$ on $\eta$ at different
points on the second-order phase boundary. The critical points are
parameterized by $r$, see Eq. (55), and we use $r=1,0.3,2$ in panels (a,b,c),
respectively. The blue circles are obtained from the master equation while the
solid curves are the analytical prediction of Eq. (54). The dashed curves omit
the $P_{-}$ prefactor. As explained in the main text, these curves may be
referred to the thermodynamic limit of the Dicke model (where
$\eta=N\Omega/\omega_{c}$ [31]). We used $\kappa=0.5\omega_{c}$.

The above results give a critical scaling:

$\displaystyle\langle x^{2}\rangle\propto\eta^{\zeta},$ (57)

with the same exponent $\zeta=1/2$ of the isotropic limit [32]. This
conclusion is hardly surprising, as a finite anisotropy is not expected to
modify the universality class of the model. Instead, it is more interesting to
perform a detailed comparison to the critical scaling of the Dicke model. The
comparison was done numerically in Ref. [32] for the isotropic case and is
immediate within the above treatment. In fact, the effective equilibrium
theory of the anisotropic Dicke model leads to the same potential and
temperature of Eqs. (49) and (50), with the $\eta$ parameter generalized as
follows [31]:

$\eta\to\frac{N\Omega}{\omega_{c}}.$ (58)

Besides this replacement, an important difference is due to spin fluctuations,
which introduce $P_{-}$ in Eq. (V). Within the effective equilibrium theory,
this is the only non-universal factor relating the two models. In particular,
our Eq. (56) differs by $P_{-}=1/2$ from the corresponding result for the
isotropic Dicke model [55]. The non-universal factor found in Ref. [32] is
0.507, close to the exact value predicted here. In the presence of a finite
anisotropy ($r\neq 1$), the non-universal factor is simply $P_{-}=1/(1+r^{2})$
and determines the difference between dashed and solid curves of Fig. 9.
Therefore, as shown in panel (b), the difference in the scaling behavior of
the Dicke and Rabi models should be strongly reduced at small $r$. In
conclusion, the above arguments give a physical interpretation of the non-
universal factor, which has been identified with the $\downarrow$ spin
population.

###  V.2 Scaling functions

Figure 10:  Finite-$\eta$ scaling function obtained from the evaluation of
$\langle\hat{x}^{2}\rangle$. The numerical data (symbols) are computed at
different values of $\eta$ and $\tilde{\lambda}_{\pm}$ and collapse to a
universal function when properly rescaled. The scaling coefficients $c_{1,2}$
are given in Eq. (61) and the solid curve is $xQ(x)$, see Eq. (52). As in
previous figures, $\kappa=0.5\omega_{c}$.

The effective equilibrium theory also allows us to describe accurately the
scaling functions which, as clear from Eq. (V), can be related to $Q(x)$.
Specifically, to make contact with previous literature [32], we first compute
from Eq. (V) the exponent $\nu$, giving the divergence of the fluctuations
when approaching the critical point in the limit $\eta\to\infty$:

$\langle\hat{x}^{2}\rangle\simeq
P_{-}\frac{T}{\omega_{c}\alpha^{2}}\propto\frac{1}{|\delta\tilde{\lambda}_{-}|^{\nu}},$
(59)

where we assumed $\tilde{\lambda}_{-}<\tilde{\lambda}_{c}^{-}$ and used
$Q(x)\simeq-\frac{1}{4x}$ when $x\to-\infty$. Expanding $\alpha^{2}$ as in Eq.
(53) we immediately find $\nu=1$. Furthermore, we derive the following scaling
relation

$(\delta\tilde{\lambda}_{-})^{\nu}\langle\hat{x}^{2}\rangle=F\left(\eta(\delta\tilde{\lambda}_{-})^{\nu/\zeta}\right),$
(60)

where $F(x)=(c_{1}/c_{2})\sqrt{x}Q(c_{1}\sqrt{x})$ can be obtained from simple
manipulations of Eq. (V). Here:

$\displaystyle
c_{1}=\frac{\Lambda}{2}\sqrt{\frac{\omega_{c}}{\beta^{3}T}},\qquad
c_{2}=\frac{\Lambda\omega_{c}}{4P_{-}T},$ (61)

where $\Lambda$ is defined by Eq. (53). The above expressions of $c_{1,2}$
allow us to collapse the scaling functions obtained at different values of $r$
along the phase boundary. In Fig. 10 we show that the collapse of numerical
data is excellent, giving a scaling function which agrees with the analytical
prediction of Eq. (V). Other scaling functions, e.g., describing
$\langle\hat{a}^{\dagger}\hat{a}\rangle$, can be obtained in a similar manner.

##  VI Discussion

We have presented a detailed characterization of various dissipative phase
transitions of the anisotropic Rabi model with cavity dissipation. A main
point brought out by our analysis is the crucial role played by spin
fluctuations in determining the properties of the system (e.g., the order
parameter), and in reaching a proper understanding of the critical behavior.
This finding is unexpected in the context of superradiant phase transitions:
In the Dicke model, the transition is realized for $N\to\infty$, i.e., when
the system approaches its classical limit. Also the ground-state phase
transition of the Rabi model is well understood in terms of small quantum
corrections around a ‘classical oscillator limit’ [27]. Instead, here large
quantum fluctuations persist in the stationary state. They lead to a general
suppression of superradiance, due to effective incoherent pumping of the
qubit, and completely determine the non-universal scaling factors of the
second-order transition. Furthermore, superradiance can disappear in certain
regions of parameters where, in principle, a mean-field superradiant state
would be stable.

We expect the type of processes described here to be generic for dissipative
transitions involving few-body systems. Qualitatively, quantum fluctations
persists due to the small Hilbert space of the qubit, which makes the high-
energy states accessible via a single spin-flip. Therefore, a similar
effective spin pumping process would likely play an important role in phase
transitions of the Dicke model with fixed $N\geq 2$ (still considering the
$\Omega/\omega_{c}\to\infty$ limit). A wider variety of related few-body
transitions, including few-site lattices and central-spin systems, has been
discussed in the recent literature [63, 64, 33, 35, 36, 37, 38, 39, 40], and
several aspects of the present analysis might be applicable to those systems
as well.

It is also interesting to make contact with other recently discussed non-
equilibrium effects. The physics of this model is strictly related to concept
of effective spin temperature, which is also important for the dissipative
Dicke model (see, e.g., Ref. [55]). In the Rabi model, however, the non-
equilibrium state is more extreme since, as mentioned, the effective spin
temperature is infinite or negative in half of the parameter space. Therefore,
a closer analogy can be perhaps found to non-equilibrium phenomena such as the
‘breakdown of photon blockade’ [63, 64] of a single spin coupled to a cavity,
or the ‘breakdown phase transition’ of the Dicke model [65], where the system
is led to a highly excited state. Especially, in the ‘breakdown phase
transition’, infinite heating of the cavity is induced by a combination of a
strong coherent interaction and qubit decay. The former can simultaneously
excite the cavity and the qubits, while dissipation tends to continuously
reset the spin state [65].

The basic mechanism at play here is the mirror-image of the one described
above. Spin pumping is induced by a combined effect of the counter-rotating
terms (which can simultaneously excite the qubit and the cavity) and cavity
relaxation. However, there are significant differences from the ‘breakdown
phase transition’. Here, for example, a strong collective interaction is not
necessary. The effective temperature only depends on the ratio of co- and
counter-rotating couplings, thus an infinite or negative spin temperature is
realized already at weak coupling. In general, we expect that dissipative
phase transitions of few-body systems may lead more easily to this type of
states under extreme non-equilibrium conditions.

Another natural implication of our treatment is that qubit relaxation should
be favorable to the formation of the superradiant state. In fact, since the
spin-flip rates due to the cavity are suppressed by $\eta$, natural
decoherence mechanisms of the qubit eventually dominate when $\eta\to\infty$.
This is actually an interesting point, as in more realistic systems the
perfect correspondence between the Dicke and Rabi model might be restored
(i.e., the non-universal factor $P_{-}$ could be absent, like in the zero-
temperature QPT [31]). On the other hand, the cavity-induced spin-flip rates
only decay as $\sim\kappa/\eta$. If the qubit relaxation rate $\gamma$ is much
smaller than $\kappa$, a sufficiency large suppression factor might be
difficult to achieve (a maximum value $\eta=25$ has been realized in Ref.
[43]). In that case, ‘finite-size’ contributions to the spin-flip transition
rates are likely to play an important role in the whole range of accessible
parameters.

Future studies may address more precisely this interplay between qubit and
cavity-induced dissipation. It would be also interesting to determine if an
enhancement of superrdiance by qubit dissipation could be helpful for recently
proposed quantum sensing protocols, based on the critical behavior of the Rabi
model [44, 45, 46, 47]

_Note added:_ During the last stages of preparing this manuscript, we became
aware of a related work studying the dissipative phase transitions of the
anisotropic quantum Rabi model [66].

##  VII Acknowledgments

This work has been supported by the ERC under grant agreement No. 101053159
(RAVE), and by the PNRR MUR project PE0000023-NQSTI. Y.D.W. acknowledges
support from the NSFC (Grant No. 12275331), the Penghuanwu Innovative Research
Center (Grant No. 12047503), and the Innovation Program for Quantum Science
and Technology (Grant No. 2021ZD031602). S.C. acknowledges support from the
Innovation Program for Quantum Science and Technology (Grant No.
2021ZD0301602), the National Science Association Funds (Grant No. U2230402),
and the National Natural Science Foundation of China (Grant Nos. 11974040 and
12150610464). The numerical solutions of the master equation were obtained
using QuTiP [67, 68].

##  Appendix A Cumulant expansion

As described in Sec. IV.1 of the main text, we have obtained a set of non-
linear equations for the second-order correlation functions. Besides Eq.
(IV.1), we have:

$\displaystyle\frac{d\langle\hat{a}^{\dagger}\hat{a}\rangle}{dt}=$
$\displaystyle-2\kappa\langle\hat{a}^{\dagger}\hat{a}\rangle-i\lambda_{-}(\langle\hat{a}\hat{\sigma}_{+}\rangle-\langle\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle)$
$\displaystyle-i\lambda_{+}(\langle\hat{a}\hat{\sigma}_{-}\rangle-\langle\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle),$
(62) $\displaystyle\frac{d\langle\hat{a}^{2}\rangle}{dt}=$
$\displaystyle-2\left(\kappa+i\omega_{c}\right)\langle\hat{a}^{2}\rangle+2i\left(\lambda_{-}\langle\hat{a}\hat{\sigma}_{-}\rangle+\lambda_{+}\langle\hat{a}\hat{\sigma}_{+}\rangle\right),$
(63) $\displaystyle\frac{d\langle\hat{a}^{\dagger 2}\rangle}{dt}=$
$\displaystyle-2\left(\kappa-i\omega_{c}\right)\langle\hat{a}^{\dagger
2}\rangle-2i(\lambda_{-}\langle\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle+\lambda_{+}\langle\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle),$
(64) $\displaystyle\frac{d\langle\hat{a}\hat{\sigma}_{+}\rangle}{dt}=$
$\displaystyle-\left(\kappa+i\omega_{c}-i\Omega\right)\langle\hat{a}\hat{\sigma}_{+}\rangle+i\lambda_{-}\langle\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{z}\rangle$
$\displaystyle+i\lambda_{+}\langle\hat{a}^{2}\hat{\sigma}_{z}\rangle+i\lambda_{-}\langle\hat{\sigma}_{-}\hat{\sigma}_{+}\rangle,$
(65)
$\displaystyle\frac{d\langle\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle}{dt}\\!=$
$\displaystyle\\!-\\!(\kappa-i\omega_{c}+i\Omega)\langle\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle-i\lambda_{-}\langle\hat{a}\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle$
$\displaystyle-i\lambda_{+}\langle\hat{a}^{\dagger
2}\hat{\sigma}_{z}\rangle-i\lambda_{-}\langle\hat{\sigma}_{+}\hat{\sigma}_{-}\rangle,$
(66) $\displaystyle\frac{d\langle\hat{a}\hat{\sigma}_{-}\rangle}{dt}=$
$\displaystyle-(\kappa+i\omega_{c}+i\Omega)\langle\hat{a}\hat{\sigma}_{-}\rangle-i\lambda_{-}\langle\hat{a}^{2}\hat{\sigma}_{z}\rangle$
$\displaystyle-i\lambda_{+}\langle(\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{z}\rangle-\langle\hat{\sigma}_{-}\hat{\sigma}_{+}\rangle),$
(67)
$\displaystyle\frac{d\langle\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle}{dt}=$
$\displaystyle-\left(\kappa-i\omega_{c}-i\Omega\right)\langle\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle+i\lambda_{-}\langle\hat{a}^{\dagger
2}\hat{\sigma}_{z}\rangle$
$\displaystyle+i\lambda_{+}(\langle\hat{a}\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle-\langle\hat{\sigma}_{+}\hat{\sigma}_{-}\rangle),$
(68) $\displaystyle\frac{d\langle\hat{a}\hat{\sigma}_{z}\rangle}{dt}=$
$\displaystyle\\!-\\!\left(\kappa+i\omega_{c}\right)\langle\hat{a}\hat{\sigma}_{z}\rangle+2i\left[\lambda_{-}(\langle\hat{a}^{2}\hat{\sigma}_{+}\rangle-\langle\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{-}\rangle)\right.$
$\displaystyle-\left.\lambda_{+}(\langle\hat{a}^{2}\hat{\sigma}_{-}\rangle-\langle\hat{a}^{\dagger}\hat{a}\hat{\sigma}_{+}\rangle)\right]-i\lambda_{-}\langle\hat{\sigma}_{-}\rangle+i\lambda_{+}\langle\hat{\sigma}_{+}\rangle,$
(69)
$\displaystyle\frac{d\langle\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle}{dt}=$
$\displaystyle\\!-\\!\left(\kappa+i\omega_{c}\right)\langle\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle-2i\left[\lambda_{-}(\langle\hat{a}^{\dagger
2}\hat{\sigma}_{-}\rangle-\langle\hat{a}\hat{a}^{\dagger}\hat{\sigma}_{+}\rangle)\right.$
$\displaystyle\left.-\lambda_{+}(\langle\hat{a}^{\dagger
2}\hat{\sigma}_{+}\rangle-\langle\hat{a}\hat{a}^{\dagger}\hat{\sigma}_{-}\rangle)\right]+i\lambda_{-}\langle\hat{\sigma}_{+}\rangle-i\lambda_{+}\langle\hat{\sigma}_{-}\rangle.$
(70)

In the above equations the third-order correlators appear explicitly. As
explained in the main text, we treat them by neglecting the third-order
moments and obtain a closed set of differential equations. Considerable
simplifications occur in the normal phase, where we can set
$\langle\hat{a}\rangle=\langle\hat{a}^{\dagger}\rangle=\langle\hat{\sigma_{\pm}}\rangle=\langle\hat{a}\hat{\sigma}_{z}\rangle=\langle\hat{a}^{\dagger}\hat{\sigma}_{z}\rangle=0$.
Using this condition, and discarding higher-order terms of the large-$\eta$
expansion, we obtain Eq. (26) of the main text, where $\Delta$ is:

$\displaystyle\Delta=\frac{4\lambda_{-}^{2}\lambda_{+}^{2}\left[4\omega_{c}^{2}\lambda_{\rm
diff}^{2}-2\omega_{c}\Omega\left(\kappa^{2}+\omega_{c}^{2}\right)\right]}{\lambda_{\rm
diff}^{8}-(\lambda_{-}^{2}+\lambda_{+}^{2})^{2}[2\omega_{c}\Omega\lambda_{\rm
diff}^{2}-\Omega^{2}(\kappa^{2}+\omega_{c}^{2})]}.$ (71)

Above, we defined $\lambda_{\rm diff}^{2}=\lambda_{-}^{2}-\lambda_{+}^{2}$.
One can see that this expression is of order $1/\eta$ by rescaling the
couplings as $\lambda_{\pm}=\tilde{\lambda}_{\pm}\lambda_{c}$, where
$\tilde{\lambda}_{\pm}$ are of order unity and
$\lambda_{c}=\sqrt{\Omega\omega_{c}}/2$. In Fig. 11 we perform a comparison of
Eq. (71) with the subleading correction to $\langle\hat{\sigma}_{z}\rangle$,
obtained from a direct solution of the master equation. As seen, there is good
agreement between numerical and analytical results.

Figure 11:  Comparison of analytical and numerical results for the stationary
value of $\langle\hat{\sigma}_{z}\rangle$. The circles are values of
$\langle\hat{\sigma}_{z}\rangle$ obtained from the exact master equation, see
Eq. (3). The dashed line is Eq. (71). We used the parameters
$\kappa=0.5\omega_{c}$ and $\tilde{\lambda}_{+}=\tilde{\lambda}_{-}=0.618$.

##  Appendix B Effective master equations

We give here some additional details on the effective master equation
approach. First we consider how the Hamiltonian is transformed. In the normal
phase, the appropriate unitary transformation $\hat{U}_{np}$ is given by Eq.
(28) which, to second order in the coupling parameters, leads to Eq. (IV.2) of
the main text. That effective Hamiltonian, besides being diagonal with respect
to $\hat{\sigma}_{z}$, is quadratic in the bosonic operators. To include
nonlinear effects due to the spin-cavity interaction, higher-order terms of
$\hat{H}_{np}$ should be computed, which are rather cumbersome to obtain
directly using $\hat{U}_{np}$. It is more convenient to rely on the method
developed in Ref. [31], from which the following effective Hamiltonian is
obtained, in terms of the quadratures
$\hat{x}=(\hat{a}^{\dagger}+\hat{a})/\sqrt{2}$ and
$\hat{p}=i(\hat{a}^{\dagger}-\hat{a})/\sqrt{2}$:

$\displaystyle\hat{H}_{eff}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\frac{\Omega}{2}\hat{\sigma}_{z}\sqrt{1+\frac{2}{\eta}\left[\tilde{\lambda}_{y}^{2}\hat{p}^{2}+\tilde{\lambda}_{x}^{2}\hat{x}^{2}-\tilde{\lambda}_{x}\tilde{\lambda}_{y}\hat{\sigma}_{z}\right]}.$
(72)

Here, to simplify the notation, we define the couplings
$\tilde{\lambda}_{x}=(\tilde{\lambda}_{+}+\tilde{\lambda}_{-})/2$ and
$\tilde{\lambda}_{y}=(\tilde{\lambda}_{+}-\tilde{\lambda}_{-})/2$. Expanding
the square root of Eq. (72) in a Taylor series, and setting
$\hat{\sigma}_{z}=-1$, we immediately obtain Eq. (45) of the main text. In Eq.
(45) we have omitted a constant, as well as negligible quadratic terms
$\propto 1/\eta$.

For the superradiant states, the transformation of the Hamiltonian is more
involved. In the following, we assume the mean-field solution
$\langle\hat{a}\rangle=\sqrt{\eta}c$, while the other broken-symmetry state
(with $c\to-c$) can be treated in a completely analogous manner. As described
in the main text, first a displacement transformation $\hat{U}_{d}$ is applied
to the photon field, to yield $\langle\hat{a}^{\dagger}\hat{a}\rangle\simeq
0$. This step leads to the effective Hamiltonian (IV.3), given in the main
text. Next, the following $\hat{U}_{sp}$ is applied:

$\displaystyle\hat{U}_{sp}=\exp\left[\frac{\cos\theta_{\rm
SP}}{2\sqrt{\eta}}\left(\tilde{\chi}_{-}\hat{a}\hat{\tau}_{+}+\tilde{\chi}_{+}\hat{a}^{\dagger}\hat{\tau}_{+}-{\rm
H.c.}\right)\right],$ (73)

which is analogous to the $\hat{U}_{np}$ transformation given in Eq. (28).
Here $\tilde{\chi}_{\pm}$ corresponds to $\tilde{\lambda}_{\pm}$ and the
appearance of $\cos\theta_{\rm SP}$ can be attributed to the modified spin-
splitting in Eq. (IV.3). A direct evaluation of the transformed Hamiltonian
gives:

$\displaystyle\hat{U}_{sp}^{\dagger}\hat{H}_{d}$
$\displaystyle\hat{U}_{sp}\simeq\hat{H}_{sp}-\frac{\omega_{c}}{2}\cos\theta_{\rm
SP}\Big{[}\tilde{\chi}_{z}\hat{a}^{\dagger}\hat{a}(\tilde{\chi}_{+}\hat{\tau}_{+}+\tilde{\chi}_{-}^{*}\hat{\tau}_{-})$
$\displaystyle+\tilde{\chi}_{z}\hat{a}^{2}(\tilde{\chi}_{+}^{*}\hat{\tau}_{-}+\tilde{\chi}_{-}\hat{\tau}_{+})+\tilde{\chi}_{z}\tilde{\chi}_{+}\hat{\tau}_{+}+{\rm
H.c.}\Big{]}.$ (74)

where we neglected terms of order $\omega_{c}/\sqrt{\eta}$. The main
contribution to Eq. (B) is $\hat{H}_{sp}$, given in Eq. (IV.3) of the main
text. Here we also include a correction to $\hat{H}_{sp}$ which is not
diagonal in $\hat{\tau}_{z}$ and originates from the spin-dependent drive of
$H_{d}$ [i.e., the last term in the first line of Eq. (IV.3)]. Such type of
terms does not appear in the normal-phase Hamiltonian $\hat{H}_{np}$ which, to
the same order of approximation, is diagonal in $\hat{\tau}_{z}$.

Despite this difference between the normal and superradiant states, the off-
diagonal terms of Eq. (B) are negligible for our discussions. First we note
that they are much smaller than the original off-diagonal terms (which in
$\hat{H}_{d}$ were of order $\sqrt{\eta}\omega_{c}$). At the same time, the
gap of $\hat{H}_{sp}$ between the $\hat{\tau}_{z}=\pm 1$ states is large (it
is approximately $\Omega=\eta\omega_{c}$). Therefore, if desired, we could
explicitly eliminate the off-diagonal corrections by applying a second unitary
transformation $\hat{U}^{\prime}_{sp}=\exp[\hat{S}^{\prime}]$, where
$\hat{S}^{\prime}\sim 1/\eta$. It is then not difficult to check that the
dominant part of the Hamiltonian, $\hat{H}_{sp}$, would remain unchanged and
all the higher-order corrections would be at most of order
$\omega_{c}/\sqrt{\eta}$. This argument justifies omitting in Eq. (IV.3) the
off-diagonal term of Eq. (B).

After having obtained Eqs. (IV.2) and (IV.3) as the appropriate effective
Hamiltonians, we discuss the effect of dissipation. Before applying
$\hat{U}_{np}$ and $\hat{U}_{sp}$, the master equations have a similar form,
i.e., Eq. (3) in the original frame and

$\displaystyle\dot{\hat{\rho}}_{d}=-i\left[\hat{H}_{d},\hat{\rho}_{d}\right]+\kappa\mathcal{D}[\hat{a}]\hat{\rho}_{d},$
(75)

after the displacement transformation, where
$\hat{\rho}_{d}=\hat{U}_{d}^{\dagger}\hat{\rho}\hat{U}_{d}$. We first discuss
the master equation for the normal states where, after applying $\hat{U}_{np}$
to Eq. (3) and neglecting corrections smaller than $\kappa/\eta$, we get:

$\displaystyle\dot{\hat{\widetilde{\rho}}}\simeq$
$\displaystyle-i\left[\hat{H}_{np},\hat{\widetilde{\rho}}\right]+\kappa\mathcal{D}[\hat{a}]\hat{\widetilde{\rho}}$
$\displaystyle-\frac{\kappa}{2\sqrt{\eta}}\left(\mathcal{D}_{C}[\hat{a},\hat{O}]\hat{\widetilde{\rho}}+\mathcal{D}_{C}[\hat{O}^{\dagger},\hat{a}^{\dagger}]\hat{\widetilde{\rho}}\right)$
$\displaystyle-\frac{\kappa}{4\eta}\tilde{\lambda}_{-}\tilde{\lambda}_{+}\left(\mathcal{D}_{C}[\hat{\sigma}_{+},\hat{\sigma}_{+}]\hat{\widetilde{\rho}}+\mathcal{D}_{C}[\hat{\sigma}_{-},\hat{\sigma}_{-}]\hat{\widetilde{\rho}}\right)$
$\displaystyle+\frac{\kappa}{4\eta}\left(\tilde{\lambda}_{-}^{2}-\tilde{\lambda}_{+}^{2}\right)\left(\mathcal{D}_{C}[\hat{a},\hat{a}]\hat{\widetilde{\rho}}+\mathcal{D}_{C}[\hat{a}^{\dagger},\hat{a}^{\dagger}]\hat{\widetilde{\rho}}\right)$
$\displaystyle+\frac{\kappa}{4\eta}\left(\tilde{\lambda}_{-}^{2}\mathcal{D}[\hat{\sigma}_{-}]\hat{\widetilde{\rho}}+\tilde{\lambda}_{+}^{2}\mathcal{D}[\hat{\sigma}_{+}]\hat{\widetilde{\rho}}\right).$
(76)

In Eq. (B),
$\mathcal{D}_{C}[\hat{O}_{1},\hat{O}_{2}]\hat{\rho}=2\hat{O}_{1}\hat{\rho}\hat{O}_{2}-\hat{O}_{1}\hat{O}_{2}\hat{\rho}-\hat{\rho}\hat{O}_{1}\hat{O}_{2}$
and (in the second line)
$\hat{O}=\tilde{\lambda}_{-}\hat{\sigma}_{-}-\tilde{\lambda}_{+}\hat{\sigma}_{+}$.
Since the first line of Eq. (B) is dominant in the limit of large $\eta$ and
$\hat{H}_{np}$ is diagonal in $\hat{\sigma}_{z}$, we consider the effect of
the remaining terms on a density matrix with the following form:

$\hat{\widetilde{\rho}}=P_{+}\hat{\rho}_{c,\uparrow}|\uparrow\rangle\langle\uparrow|+P_{-}\hat{\rho}_{c,\downarrow}|\downarrow\rangle\langle\downarrow|,$
(77)

where $\hat{\rho}_{c,\sigma}$ are spin-resolved density matrices for the
cavity mode.

Formally, the second line of Eq. (B) is the dominant correction to the
dissipator, as it is only suppressed by $\sqrt{\eta}$. However, the effect of
applying this term to Eq. (77) is to induce coherence between the $\uparrow$
and $\downarrow$ subspaces. Due to the large gap in $\hat{H}_{np}$, these off-
diagonal terms are highly oscillating, thus they can be neglected. As for the
third line of Eq. (B), it is easily seen that it is zero on a diagonal
$\hat{\widetilde{\rho}}$ like Eq. (77). Finally, the fourth line of Eq. (B)
has an influence on the detailed form of the cavity density matrices
$\hat{\rho}_{c,\sigma}$, but it does not introduce qualitatively new physical
effects. At large $\eta$, it is negligible compared to the regular photon
decay, included in the first line. We conclude that the non-Lindblad terms in
the second, third and fourth lines of Eq. (B) have little effect on our
results, thus are omitted in the effective master equation. On the other hand,
the last term of Eq. (B) is crucial to establish the populations $P_{\pm}$ of
the spin states, thus is included in Eq. (30).

Finally, we briefly consider the superradiant fixed point. Due to the
similarity of $\hat{U}_{np}$ and $\hat{U}_{sp}$, the transformation of the
dissipator of Eq. (75) can be performed in a way completely analogous to the
normal phase, which leads to the effective master equation for the
superradiant states Eq. (IV.3).

##  Appendix C Effective potential and noise at the second-order critical
point

Close to the second-order transition between the normal and superradiant
phases, we base our treatment on the effective Hamiltonian in the
$\hat{\sigma}_{z}=-1$ subspace, given by Eq. (45). The associated Heisenberg-
Langevin equations read:

$\displaystyle\frac{d\hat{x}}{dt}$
$\displaystyle=\omega_{c}\left[\alpha_{y}\hat{p}+\frac{\tilde{\lambda}_{y}^{4}\hat{p}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{p},\hat{x}^{2}\\}}{\eta}\right]-\kappa\hat{x}+\hat{\mathcal{F}}_{+},$
(78) $\displaystyle\frac{d\hat{p}}{dt}$
$\displaystyle=-\omega_{c}\left[\alpha_{x}\hat{x}+\frac{\tilde{\lambda}_{x}^{4}\hat{x}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{x},\hat{p}^{2}\\}}{\eta}\right]-\kappa\hat{p}-i\hat{\mathcal{F}}_{-},$
(79)

where we defined
$\\{\hat{A},\hat{B}\\}=\frac{1}{2}(\hat{A}\hat{B}+\hat{B}\hat{A})$ and
$\hat{\mathcal{F}}_{\pm}=(\hat{\mathcal{F}}\pm\hat{\mathcal{F}}^{\dagger})/\sqrt{2}$,
with $\hat{\mathcal{F}}(t)$ the same noise operator entering Eq. (II). It is
also useful to rewrite Eq. (78) as follows:

$\hat{p}=\frac{\frac{d\hat{x}}{dt}+\kappa\hat{x}-\hat{\mathcal{F}}_{+}}{\omega_{c}\alpha_{y}}-\frac{\tilde{\lambda}_{y}^{4}\hat{p}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{p},\hat{x}^{2}\\}}{\eta\alpha_{y}}.$
(80)

Next, we compute $d^{2}\hat{x}/dt^{2}$ from Eq. (78), leading to:

$\displaystyle\frac{d^{2}\hat{x}}{dt^{2}}=$
$\displaystyle-\omega_{c}^{2}\alpha^{2}\hat{x}-2\kappa\frac{d\hat{x}}{dt}+\mathcal{G}$
$\displaystyle-\frac{\omega_{c}^{2}}{\eta}\alpha_{y}\left(\tilde{\lambda}_{x}^{4}\hat{x}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{x},\hat{p}^{2}\\}\right)$
$\displaystyle+\frac{\kappa\omega_{c}}{\eta}\left(\tilde{\lambda}_{y}^{4}\hat{p}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{p},\hat{x}^{2}\\}\right)$
$\displaystyle+\frac{\omega_{c}}{\eta}\frac{d}{dt}\left(\tilde{\lambda}_{y}^{4}\hat{p}^{3}+\tilde{\lambda}_{x}^{2}\tilde{\lambda}_{y}^{2}\\{\hat{p},\hat{x}^{2}\\}\right),$
(81)

where in the first line $\alpha^{2}$ is given by Eq. (48) and the noise
operator is:

$\hat{\mathcal{G}}(t)=\frac{1}{\sqrt{2}}\left[(\kappa-i\alpha_{-}\omega_{c})\hat{\mathcal{F}}(t)+\frac{d\hat{\mathcal{F}}(t)}{dt}+{\rm
H.c.}\right].$ (82)

Note that, to obtain Eq. (C) we have substituted $d\hat{p}/dt$ with Eq. (79),
where $\hat{p}$ is given by Eq. (80).

We should now simplify the nonlinear terms of Eq. (C). In the second and third
lines, we could in principle substitute the full expression of $\hat{p}$,
given by Eq. (80). However, doing this leads to a number of uninteresting
corrections. Firstly, we can safely neglect in Eq. (C) the higher-order terms,
$\propto 1/\eta^{2}$, thus the $1/\eta$ term of Eq. (80) can be omitted.
Furthermore, Eq. (80) would induce $1/\eta$ corrections to the damping and
noise terms, appearing in the first line of Eq. (C). These effects are not
important when $\eta\to\infty$. We conclude that, to deal with the nonlinear
terms, we can effectively simplify Eq. (80) to
$\hat{p}\simeq\kappa\hat{x}/(\omega_{c}\alpha_{y})$. After this substitution,
and omitting the last line of Eq. (C), we obtain Eq. (47) of the main text,
which for convenience we repeat here:

$\frac{d^{2}\hat{x}}{dt^{2}}=-\omega_{c}^{2}\left(\alpha^{2}\hat{x}+\frac{\beta^{3}}{\eta}\hat{x}^{3}\right)-2\kappa\frac{d\hat{x}}{dt}+\hat{\mathcal{G}}.$
(83)

The coefficient $\beta^{2}$, obtained from the second and third lines of Eq.
(C), is given in Eq. (48). Finally, we justify why the fourth line Eq. (C) was
omitted. Due to the time derivative, this contribution will give two types of
terms. The ones proportional to $d\hat{x}/dt$, following our previous
arguments, can be neglected in comparison to the dominant damping term. The
second type of terms is $\propto d{\hat{p}}/dt$, and we should express them
using Eq. (79). However, to obtain the effective potential we can neglect the
$1/\eta$ and noise contributions, thus Eq. (79) simplifies to
$d{\hat{p}}/dt\simeq\omega_{c}\alpha_{x}\hat{x}-\kappa\hat{p}$ and,
eventually:

$\frac{d\hat{p}}{dt}\simeq\omega_{c}\alpha_{x}\hat{x}-\kappa\left(\frac{\kappa\hat{x}}{\omega_{c}\alpha_{y}}\right)=\frac{\omega_{c}}{\alpha_{y}}\alpha^{2}\hat{x}.$
(84)

In conclusion, we find that, within these approximations, $d{\hat{p}}/dt$
vanishes at the critical point (where $\alpha^{2}=0$), thus the corresponding
contribution to the effective potential can be omitted.

We conclude this section by deriving the effective temperature of the bath,
obtained by a comparison to the following equation of motion of a Brownian
particle (we set $k_{\rm B}=1$):

$m\frac{d^{2}y}{dt^{2}}=-\frac{dV(y)}{dy}-\gamma\frac{d^{2}y}{dt^{2}}+\sqrt{2\gamma
T}\xi,$ (85)

where $m$ is the mass, $V(y)$ is the potential, and $\gamma$ is the damping
coefficient. The noise term satisfies
$\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$ and
$\langle\xi(t)\rangle=0$. Since here $y$ is a position and $\hat{x}$ in Eq.
(83) is dimensionless we define $x=y/y_{0}$, where $y_{0}$ is an arbitrary
conversion length. Equivalently, we choose an energy scale $\omega_{0}$ giving
$y_{0}^{-1}=\sqrt{\omega_{0}m}$. First, comparing Eq. (85) to the classical
limit of Eq. (83), we find the potential in the $x$ coordinate:

$\displaystyle
V(x)=\frac{\omega_{c}^{2}}{\omega_{0}}\left(\frac{1}{2}\alpha^{2}x^{2}+\frac{1}{4\eta}\beta^{3}x^{4}\right).$
(86)

Furthermore, proceeding as in Ref. [55], we compute the symmetrized correlator
of the noise in Fourier space:

$\frac{1}{2}\langle\hat{\mathcal{G}}_{\omega}\hat{\mathcal{G}}_{\omega^{\prime}}+\hat{\mathcal{G}}_{\omega^{\prime}}\hat{\mathcal{G}}_{\omega}\rangle=\kappa(\kappa^{2}+\alpha^{2}_{y}\omega_{c}^{2}+\omega^{2})\delta(\omega+\omega^{\prime}),$
(87)

where
$\hat{\mathcal{G}}_{\omega}=\int\frac{dt}{\sqrt{2\pi}}\hat{\mathcal{G}}(t)e^{i\omega
t}$. The analogous correlator of the classical noise is $\frac{\sqrt{2\gamma
T}}{my_{0}}\delta(\omega+\omega^{\prime})$. Considering the low-frequency
limit and using $\gamma/m=2\kappa$, we finally identify:

$T=\frac{\kappa^{2}+\alpha^{2}_{y}\omega_{c}^{2}}{4\omega_{0}}.$ (88)

We see that the choice of the energy scale is arbitrary, as $\omega_{0}$
simplifies when computing $V(x)/T$. The effective temperature of Ref. [55] is
recovered in the isotropic limit (giving $\alpha_{y}=1$) and with the natural
choice $\omega_{0}=\omega_{c}$, which we assume in the main text.

## References

  * Raimond  _et al._ [2001] J. M. Raimond, M. Brune, and S. Haroche, Rev. Mod. Phys. 73, 565 (2001).
  * Ritsch  _et al._ [2013] H. Ritsch, P. Domokos, F. Brennecke, and T. Esslinger, Rev. Mod. Phys. 85, 553 (2013).
  * Fitzpatrick  _et al._ [2017] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch, and A. A. Houck, Phys. Rev. X 7, 011016 (2017).
  * Collodo  _et al._ [2019] M. C. Collodo, A. Potočnik, S. Gasparinetti, J.-C. Besse, M. Pechal, M. Sameti, M. J. Hartmann, A. Wallraff, and C. Eichler, Phys. Rev. Lett. 122, 183601 (2019).
  * Smith  _et al._ [2016] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nat. Phys. 12, 907 (2016).
  * Guo  _et al._ [2021] Q. Guo, C. Cheng, Z.-H. Sun, Z. Song, H. Li, Z. Wang, W. Ren, H. Dong, D. Zheng, Y.-R. Zhang, R. Mondaini, H. Fan, and H. Wang, Nat. Phys. 17, 234 (2021).
  * Zhang  _et al._ [2017a] J. Zhang, P. W. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I. D. Potirniche, A. C. Potter, A. Vishwanath, N. Y. Yao, and C. Monroe, Nature 543, 217 (2017a).
  * Jurcevic  _et al._ [2017] P. Jurcevic, H. Shen, P. Hauke, C. Maier, T. Brydges, C. Hempel, B. P. Lanyon, M. Heyl, R. Blatt, and C. F. Roos, Phys. Rev. Lett. 119, 080501 (2017).
  * Zhang  _et al._ [2017b] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z. X. Gong, and C. Monroe, Nature 551, 601 (2017b).
  * Monroe  _et al._ [2021] C. Monroe, W. C. Campbell, L.-M. Duan, Z.-X. Gong, A. V. Gorshkov, P. W. Hess, R. Islam, K. Kim, N. M. Linke, G. Pagano, P. Richerme, C. Senko, and N. Y. Yao, Rev. Mod. Phys. 93, 025001 (2021).
  * Blais  _et al._ [2021] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Rev. Mod. Phys. 93, 025005 (2021).
  * Borjans  _et al._ [2020] F. Borjans, X. G. Croot, X. Mi, M. J. Gullans, and J. R. Petta, Nature 577, 195 (2020).
  * Harvey-Collard  _et al._ [2022] P. Harvey-Collard, J. Dijkema, G. Zheng, A. Sammak, G. Scappucci, and L. M. K. Vandersypen, Phys. Rev. X 12, 021026 (2022).
  * Dicke [1954] R. H. Dicke, Phys. Rev. 93, 99 (1954).
  * Kirton  _et al._ [2019] P. Kirton, M. M. Roses, J. Keeling, and E. G. Dalla Torre, Adv. Quantum Technol. 2, 1800043 (2019).
  * Hepp and Lieb [1973] K. Hepp and E. H. Lieb, Ann. Phys. (N.Y.) 76, 360 (1973).
  * Wang and Hioe [1973] Y. K. Wang and F. T. Hioe, Phys. Rev. A 7, 831 (1973).
  * Dimer  _et al._ [2007] F. Dimer, B. Estienne, A. S. Parkins, and H. J. Carmichael, Phys. Rev. A 75, 013804 (2007).
  * Kirton and Keeling [2017] P. Kirton and J. Keeling, Phys. Rev. Lett. 118, 123602 (2017).
  * Baumann  _et al._ [2010] K. Baumann, C. Guerlin, F. Brennecke, and T. Esslinger, Nature 464, 1301 (2010).
  * Klinder  _et al._ [2015] J. Klinder, H. Keßler, M. Wolke, L. Mathey, and A. Hemmerich, Proc. Natl. Acad. Sci. U.S.A. 112, 3290 (2015).
  * Zhiqiang  _et al._ [2017] Z. Zhiqiang, C. H. Lee, R. Kumar, K. J. Arnold, S. J. Masson, A. S. Parkins, and M. D. Barrett, Optica 4, 424 (2017).
  * Zhang  _et al._ [2018] Z. Zhang, C. H. Lee, R. Kumar, K. J. Arnold, S. J. Masson, A. L. Grimsmo, A. S. Parkins, and M. D. Barrett, Phys. Rev. A 97, 043858 (2018).
  * Safavi-Naini  _et al._ [2018] A. Safavi-Naini, R. J. Lewis-Swan, J. G. Bohnet, M. Gärttner, K. A. Gilmore, J. E. Jordan, J. Cohn, J. K. Freericks, A. M. Rey, and J. J. Bollinger, Phys. Rev. Lett. 121, 040503 (2018).
  * Levine and Muthukumar [2004] G. Levine and V. Muthukumar, Phys. Rev. B 69, 113203 (2004).
  * Hines  _et al._ [2004] A. P. Hines, C. M. Dawson, R. H. McKenzie, and G. J. Milburn, Phys. Rev. A 70, 022303 (2004).
  * Bakemeier  _et al._ [2012] L. Bakemeier, A. Alvermann, and H. Fehske, Phys. Rev. A 85, 043821 (2012).
  * Ashhab [2013] S. Ashhab, Phys. Rev. A 87, 013826 (2013).
  * Hwang  _et al._ [2015] M.-J. Hwang, R. Puebla, and M. B. Plenio, Phys. Rev. Lett. 115, 180404 (2015).
  * Puebla  _et al._ [2017] R. Puebla, M.-J. Hwang, J. Casanova, and M. B. Plenio, Phys. Rev. Lett. 118, 073001 (2017).
  * Liu  _et al._ [2017] M. Liu, S. Chesi, Z.-J. Ying, X. Chen, H.-G. Luo, and H.-Q. Lin, Phys. Rev. Lett. 119, 220601 (2017).
  * Hwang  _et al._ [2018] M.-J. Hwang, P. Rabl, and M. B. Plenio, Phys. Rev. A 97, 013825 (2018).
  * Hwang and Plenio [2016] M.-J. Hwang and M. B. Plenio, Phys. Rev. Lett. 117, 123602 (2016).
  * Felicetti and Le Boité [2020] S. Felicetti and A. Le Boité, Phys. Rev. Lett. 124, 040404 (2020).
  * Wang  _et al._ [2020] Y. Wang, M. Liu, W.-L. You, S. Chesi, H.-G. Luo, and H.-Q. Lin, Phys. Rev. A 101, 063843 (2020).
  * Shao  _et al._ [2023] L. Shao, R. Zhang, W. Lu, Z. Zhang, and X. Wang, Phys. Rev. A 107, 013714 (2023).
  * Zhang  _et al._ [2021] Y.-Y. Zhang, Z.-X. Hu, L. Fu, H.-G. Luo, H. Pu, and X.-F. Zhang, Phys. Rev. Lett. 127, 063602 (2021).
  * Fallas Padilla  _et al._ [2022] D. Fallas Padilla, H. Pu, G.-J. Cheng, and Y.-Y. Zhang, Phys. Rev. Lett. 129, 183602 (2022).
  * Grimaudo  _et al._ [2023] R. Grimaudo, A. S. M. a. de Castro, A. Messina, E. Solano, and D. Valenti, Phys. Rev. Lett. 130, 043602 (2023).
  * Huang  _et al._ [2023] W. Huang, Y. Wu, and X.-Y. Lü, Phys. Rev. A 107, 033702 (2023).
  * De Filippis  _et al._ [2023] G. De Filippis, A. de Candia, G. Di Bello, C. A. Perroni, L. M. Cangemi, A. Nocera, M. Sassetti, R. Fazio, and V. Cataudella, Phys. Rev. Lett. 130, 210404 (2023).
  * Lv  _et al._ [2018] D. Lv, S. An, Z. Liu, J.-N. Zhang, J. S. Pedernales, L. Lamata, E. Solano, and K. Kim, Phys. Rev. X 8, 021027 (2018).
  * Cai  _et al._ [2021] M. L. Cai, Z. D. Liu, W. D. Zhao, Y. K. Wu, Q. X. Mei, Y. Jiang, L. He, X. Zhang, Z. C. Zhou, and L. M. Duan, Nat. Comm. 12, 1126 (2021).
  * Garbe  _et al._ [2020] L. Garbe, M. Bina, A. Keller, M. G. A. Paris, and S. Felicetti, Phys. Rev. Lett. 124, 120504 (2020).
  * Chu  _et al._ [2021] Y. Chu, S. Zhang, B. Yu, and J. Cai, Phys. Rev. Lett. 126, 010502 (2021).
  * Ilias  _et al._ [2022] T. Ilias, D. Yang, S. F. Huelga, and M. B. Plenio, PRX Quantum 3, 010354 (2022).
  * Zhu  _et al._ [2023] X. Zhu, J.-H. Lü, W. Ning, F. Wu, L.-T. Shen, Z.-B. Yang, and S.-B. Zheng, Sci. China Phys. Mech. Astron. 66, 250313 (2023), arXiv:2302.13676 [quant-ph] .
  * Keeling  _et al._ [2010] J. Keeling, M. J. Bhaseen, and B. D. Simons, Phys. Rev. Lett. 105, 043001 (2010).
  * Soriente  _et al._ [2018] M. Soriente, T. Donner, R. Chitra, and O. Zilberberg, Phys. Rev. Lett. 120, 183603 (2018).
  * Stitely  _et al._ [2020] K. C. Stitely, A. Giraldo, B. Krauskopf, and S. Parkins, Phys. Rev. Research 2, 033131 (2020).
  * Baksic and Ciuti [2014] A. Baksic and C. Ciuti, Phys. Rev. Lett. 112, 173601 (2014).
  * Xie  _et al._ [2014] Q.-T. Xie, S. Cui, J.-P. Cao, L. Amico, and H. Fan, Phys. Rev. X 4, 021046 (2014).
  * Yang and Wang [2017] W.-J. Yang and X.-B. Wang, Phys. Rev. A 95, 043823 (2017).
  * Shen  _et al._ [2017] L.-T. Shen, Z.-B. Yang, H.-Z. Wu, and S.-B. Zheng, Phys. Rev. A 95, 013819 (2017).
  * Dalla Torre  _et al._ [2013] E. G. Dalla Torre, S. Diehl, M. D. Lukin, S. Sachdev, and P. Strack, Phys. Rev. A 87, 023831 (2013).
  * Wang  _et al._ [2018] Y. Wang, W.-L. You, M. Liu, Y.-L. Dong, H.-G. Luo, G. Romero, and J. You, New J. Phys. 20, 053061 (2018).
  * Cheng  _et al._ [2018] X.-H. Cheng, I. Arrazola, J. S. Pedernales, L. Lamata, X. Chen, and E. Solano, Phys. Rev. A 97, 023624 (2018).
  * Zhong  _et al._ [2019] P.-G. Zhong, C. Li, Y. Wang, J. Song, S.-T. Liu, Y.-Y. Jiang, and Y. Xia, Phys. Rev. A 99, 043829 (2019).
  * Öztop  _et al._ [2012] B. Öztop, M. Bordyuh, Ö. E. Müstecaplıoğlu, and H. E. Türeci, New J. Phys. 14, 085011 (2012).
  * Kirton and Keeling [2018] P. Kirton and J. Keeling, New J. Phys. 20, 015009 (2018).
  * Nagy and Domokos [2015] D. Nagy and P. Domokos, Phys. Rev. Lett. 115, 043601 (2015).
  * Sieberer  _et al._ [2016] L. M. Sieberer, M. Buchhold, and S. Diehl, Rep. Prog. Phys. 79, 096001 (2016).
  * Carmichael [2015] H. Carmichael, Phys. Rev. X 5, 031028 (2015).
  * Fink  _et al._ [2017] J. M. Fink, A. Dombi, A. Vukics, A. Wallraff, and P. Domokos, Phys. Rev. X 7, 011012 (2017).
  * Reiter  _et al._ [2020] F. Reiter, T. L. Nguyen, J. P. Home, and S. F. Yelin, Phys. Rev. Lett. 125, 233602 (2020).
  * Lyu  _et al._ [2023] G. Lyu, K. Kottmann, M. B. Plenio, and M.-J. Hwang, (2023), arXiv:2311.11346 .
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, Comput. Phys. Commun. 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, Comput. Phys. Commun. 184, 1234 (2013).