# Unconventional Dicke model: Multistabilities and nonequilibrium dynamics

Farokh Mivehvar<EMAIL_ADDRESS>Institut für Theoretische Physik,
Universität Innsbruck, A-6020 Innsbruck, Austria

###### Abstract

The Dicke model describes the collective behavior of a sub-wavelength–size
ensemble of two-level atoms (i.e., spin-1/2) interacting identically with a
single quantized radiation field of a cavity. Across a critical coupling
strength it exhibits a zero-temperature phase transition from the normal state
to the superradian phase where the field is populated and and the collective
spin acquires a nonzero $x$-component, which can be imagined as
_ferromagnetic_ ordering of the atomic spins along $x$. Here we introduce a
variant of this model where two sub-wavelength–size ensembles of spins
interact with a single quantized radiation field with opposite strengths. In
addition to the common $x$-ferromagnetic superradiance, we analytically find
an exotic superradiant state with $x$-_ferrimagnetic_ ordering, coexisting
with $x$-ferromagnetically ordered superradiant state in large parameter
regimes. The stability and dynamics of the system in the thermodynamic limit
are then examined using a semiclassical approach, which predicts non-
stationary behaviors due to the multistabilities. At the end, we also perform
small-scale full quantum-mechanical calculations, with results consistent with
the semiclassical ones.

_Introduction._ —The nonequilibrium Dicke model is one of the most celebrated
models in quantum optics Hepp and Lieb (1973a, b); Wang and Hioe (1973). It
describes the collective behavior of $N$ two-level atoms (i.e., spin-1/2)
cooperatively interacting with a single quantized radiation field. The model
has a fairly simple Hamiltonian (we set $\hbar=1$ throughout the paper)
Garraway (2011),

$\displaystyle\hat{H}_{\rm
D}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{a}\hat{S}_{z}+\lambda(\hat{a}^{\dagger}+\hat{a})\hat{S}_{x},$
(1)

where $\omega_{c}$ is the frequency of the cavity mode, $\omega_{a}$ the
transition frequency of one of the atoms, and $\lambda$ the single-atom–field
coupling strength. Here, $\hat{a}$ is the bosonic annihilation operator of the
cavity radiation field and
$\hat{\mathbf{S}}=(\hat{S}_{x},\hat{S}_{y},\hat{S}_{z})$ is the collective
atomic spin operator of length $N/2$ with components
$\hat{S}_{\alpha}=\sum_{j=1}^{N}\hat{s}_{j,\alpha}$, where
$\alpha=\\{x,y,z\\}$, defined in terms of the common single-atom spin
operators $\hat{s}_{j,\alpha}$. Despite its simple form, the Dicke Hamiltonian
(1) is expected to exhibit a variety of interesting phenomena, most notably
the zero-temperature phase transition from the normal (N) state with the field
being in the vacuum state to the superradiant (SR) phase where the field
acquires a nonzero photon number $\langle\hat{a}^{\dagger}\hat{a}\rangle\neq
0$ across the critical coupling strength
$\sqrt{N}\lambda_{c}=\sqrt{\omega_{c}\omega_{a}}$ Kirton  _et al._ (2018).
Correspondingly, the collective atomic spin completely in the spin-down state
$\langle\hat{S}_{z}\rangle=N/2$ initially—i.e., _ferromagnetically_ (Fo)
ordered along the negative $z$ direction, referred to as “$-$zFo-N” in this
work—obtains a nonzero $x$-component $\langle\hat{S}_{x}\rangle\neq 0$. The
collective atomic spin orients almost completely in the positive or negative
$x$ direction $\langle\hat{S}_{x}\rangle\to\pm N/2$, depending on the broken
$\mathbb{Z}_{2}$ parity symmetry of the Dicke Hamiltonian on the onset of the
superradiant phase transition, in the deep superradiant phase
$\lambda\to\infty$. Therefore, we refer to this superradiant phase with its
asymptotic _ferromagnetic_ ordering along the positive or negative $x$
direction as “$+$xFo-SR” and “$-$xFo-SR”, respectively, designated
collectively by “$\pm$xFo-SR” (the “$+$” and “$-$” signs specify the direction
of the total spin along a desired axis).

The Dicke model has been successfully implemented using cavity-assisted two-
photon Raman transitions between atomic momentum states for both bosonic
Baumann  _et al._ (2010); Klinder  _et al._ (2015); Kollár  _et al._ (2017)
and fermionic Zhang  _et al._ (2021); Helson  _et al._ (2023) atoms as well as
using cavity-assisted two-photon Raman transitions between atomic hyperfine
ground states Zhiqiang  _et al._ (2017), bypassing the no-go theorems
Bialynicki-Birula and Rza¸żewski (1979); Nataf and Ciuti (2010). Generalized
implementations of it including both momentum and hyperfine internal states
have also been proposed Safaei  _et al._ (2013); Mivehvar  _et al._ (2017,
2019); Masalaeva  _et al._ (2021) and realized experimentally Landini  _et
al._ (2018); Kroeze  _et al._ (2018, 2019) in quantum-gas–cavity platforms.
There is also great interest in realizing Dicke-type models and superradiance
in waveguide-QED setups Chang  _et al._ (2018); Sheremet  _et al._ (2023) and
cavity quantum materials Schlawin  _et al._ (2022).

Motivated by the recent progress in quantum-gas–cavity QED Mivehvar  _et al._
(2021), in this Letter we introduce a variant of the Dicke model, which we
coin the name “ _non-standard_ Dicke model” for it, where two independent
ensembles of spin-1/2 atoms are coupled to a single mode of a cavity, but with
opposite coupling-strength signs. The system exhibits a wealth of intriguing
phenomena. In particular, the semi-classical approach reveals the existence of
multistable steady-state phases as shown in Figs. 1 and 2; especially, a
bistable superradiant regime. This bistable region contains $\pm$xFo-SR states
where the total spins in the both ensembles orient in the same $x$ direction
(i.e., both in either positive or negative $x$ direction, respectively) in the
strong coupling limit. However, it also includes more exotic superradiant
phases where the total spins in the two ensembles point in opposite $x$
directions (i.e., one in the positive and another in the negative $x$
directions) in the strong coupling limi, thus forming _ferrimagnetic_ (Fi)
order “$\pm$xFi-SR”. This order evolves into _antiferromagnetic_ (aF) order,
“$\pm$xaF-SR”, when two ensembles have an equal number of spins.

The linear-stability analysis as well as the nonequilibrium dynamics of the
system assert the stability of the steady states of the system. Furthermore,
we find initial states in the multistable regimes where the ensuing
nonequilibrium dynamics from them do not lead to any steady state of the
system, rather give rise to non-stationary oscillating trajectories due to
competing fixed points in these regimes; see Fig. 3. Finally, the full
quantum-mechanical calculations also confirm the coexistence of $\pm$xFo-SR
and $\pm$xFi-SR as shown in Fig. 4.

_Model and Hamiltonian._ —Consider two independent ensembles of $N_{1,2}$
spins-1/2 located at negative and positive maxima of a cavity field and
coupled to it with strengths $\pm\lambda$, respectively. The system is
described by a non-standard Dicke Hamiltonian,

$\displaystyle\hat{H}_{\rm
nsD}=\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{a}\hat{S}_{z}+\lambda(\hat{a}^{\dagger}+\hat{a})\Delta\hat{S}_{x},$
(2)

where
$\hat{S}_{\alpha}\equiv\hat{S}_{1,\alpha}+\hat{S}_{2,\alpha}=\sum_{j=1}^{N_{1}}\hat{s}_{j,\alpha}^{(1)}+\sum_{j=1}^{N_{2}}\hat{s}_{j,\alpha}^{(2)}$
is the total collective spin of the two ensembles and
$\Delta\hat{S}_{\alpha}\equiv\hat{S}_{1,\alpha}-\hat{S}_{2,\alpha}$ is a
staggered collective spin. Unlike the Dicke model (1), the total collective
spin here is not a constant of motion as
$\hat{\mathbf{S}}^{2}=\hat{S}_{x}^{2}+\hat{S}_{y}^{2}+\hat{S}_{z}^{2}$ does
not commute with $\hat{H}_{\rm nsD}$. This means that the non-standard Dicke
Hamiltonian $\hat{H}_{\rm nsD}$ mixes manifolds with different total spins
Hotter  _et al._ (2023). That said, the total spin for each ensemble is
conserved. As we show in supplemental material (SM) Miv , remarkably the non-
standard Dicke model also describes nonequilibrium dynamics of the combined
density-spin self-ordering of transversely-driven two-component quantum gases
inside single-mode cavities deep in the superradiant phase Safaei  _et al._
(2013); Mivehvar  _et al._ (2017, 2019); Masalaeva  _et al._ (2021); Landini
_et al._ (2018); Kroeze  _et al._ (2018, 2019). Due to the symmetry, in the
following we restrict ourselves solely to $N_{2}/N_{1}\leq 1$.

_Mean-field approach: steady states, stability, and dynamics._ —In order to
obtain insight into the system, we start with the mean-field approach that
omits quantum fluctuations and replaces quantum operators with classical
variables, namely, $\hat{a}\to a=\langle\hat{a}\rangle$ and
$\hat{\mathbf{S}}_{l}\to\mathbf{S}_{l}=\langle\hat{\mathbf{S}}_{l}\rangle$
with $l=1,2$, justified for large ensembles of spins Keeling  _et al._ (2010).
The system is then described by a set of seven coupled differential equations
obtained from Heisenberg equations of motion,

$\displaystyle\dot{a}$
$\displaystyle=-(i\omega_{c}+\kappa)a-i\lambda(S_{1,x}-S_{2,x}),$
$\displaystyle\dot{S}_{l,x}$ $\displaystyle=-\omega_{a}S_{l,y},$
$\displaystyle\dot{S}_{l,y}$
$\displaystyle=\omega_{a}S_{l,x}+(-1)^{l}\lambda(a^{*}+a)S_{l,z},$
$\displaystyle\dot{S}_{l,z}$
$\displaystyle=(-1)^{l+1}\lambda(a^{*}+a)S_{l,y},$ (3)

and endowed by two spin-conservation constraints
$S_{l}=|{\mathbf{S}}_{l}|=\sqrt{S_{l,x}^{2}+S_{l,y}^{2}+S_{l,z}^{2}}=N_{l}/2$.
Here, $\kappa$ is the cavity-field decay rate.

The fixed-point solutions are obtained by setting $\dot{a}=0$ and
$\dot{S}_{l,\alpha}=0$ in Eq. (Unconventional Dicke model: Multistabilities
and nonequilibrium dynamics). This immediately implies that in a steady state,
$S_{l,y}^{\rm(ss)}=0$ for any $\omega_{a}\neq 0$. After some algebraic
manipulation Miv , one obtains two distinct classes of nontrivial solutions:

$\displaystyle
S_{1,x}^{\rm(ss)}=\frac{N_{1}}{N_{2}}S_{2,x}^{\rm(ss)}=\pm\frac{N_{1}}{2}\sqrt{1-\bigg{(}\frac{\lambda_{c}^{(\rm
xFo)}}{\lambda}\bigg{)}^{4}},$ (4a) with corresponding
$S_{l,z}^{\rm(ss)}=(-1)^{l}\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$
determined through the normalization constraints, and $\displaystyle
S_{1,x}^{\rm(ss)}=-\frac{N_{1}}{N_{2}}S_{2,x}^{\rm(ss)}=\pm\frac{N_{1}}{2}\sqrt{1-\bigg{(}\frac{\lambda_{c}^{(\rm
xFi)}}{\lambda}\bigg{)}^{4}},$ (4b)

with corresponding
$S_{l,z}^{\rm(ss)}=-\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$. Here, we
have introduced the critical couplings as

$\displaystyle\lambda_{c}^{(\rm
xFo/xFi)}=\sqrt{\frac{\omega_{a}(\omega_{c}^{2}+\kappa^{2})}{(N_{1}\mp
N_{2})\omega_{c}}}.$ (5)

The choice of the sign for $S_{l,x}^{\rm(ss)}$, i.e., the “$\pm$” sings
outside of the square roots in the right-hand sides of Eq. (4), indicates two
possible $\mathbb{Z}_{2}$ solutions originating from the invariance of the
Hamiltonian (2) under the parity transformation $\hat{a}\to-\hat{a}$ and
$\hat{S}_{l,x}\to-\hat{S}_{l,x}$.

Figure 1: Steady-state phase diagram of the system and schematic
representations of all phases. The two white curves obtained by rearranging
Eq. (5) in the phase diagram indicate the thresholds for the emergence of the
$\pm$xFi-SR and $\pm$xFo-SR states. The numbers inside squares show the total
number of fixed points in each regime. However, only the stable solutions are
stated explicitly. All the phases are illustrated schematically for finite
$\lambda$ in the table in the right. The other parameters are set to
$\omega_{c}=\omega_{a}=\kappa$.

The two nontrivial classes of the solutions, Eqs. (4a) and (4b), have
fundamentally different properties. The first class corresponds to
superradiant states with $x$-ferromagnetic ordering where both
$S_{1,x}^{\rm(ss)}$ and $S_{2,x}^{\rm(ss)}$ point in the same direction, while
the second class corresponds to exotic superradiant states with
$x$-ferrimagnetic ordering where $S_{1,x}^{\rm(ss)}$ and $S_{2,x}^{\rm(ss)}$
point in the opposite directions. We refer to these superradiant states,
respectively, as “$\pm$xFo-SR” and “$\pm$xFi-SR”, corresponding to their
asymptotic spin-ordering behaviors at $\lambda\to\infty$; see the table in the
right side of Fig. 1 for the schematic representation of these states at
finite $\lambda$. To re-express it plainly, in the $+$xFo-SR ($-$xFo-SR) state
both $S_{1,x}^{\rm(ss)}$ and $S_{2,x}^{\rm(ss)}$ have the same plus (minus)
$\mathbb{Z}_{2}$ parity-symmetry sign, while in the +xFi-SR ($-$xFi-SR) state
$S_{1,x}^{\rm(ss)}$ has the plus (minus) sign and $S_{2,x}^{\rm(ss)}$ the
minus (plus) sign. The $\pm$xFi-SR states cross over into
$x$-antiferromagnetic superradiant states, “$\pm$xaF-SR”, for $N_{1}=N_{2}$.

In addition to the four superradiant fixed points discussed above, there are
four trivial fixed points corresponding to $a_{\rm ss}=S_{l,x}^{\rm(ss)}=0$
and either $S_{z}^{\rm(ss)}=\pm(N_{1}+N_{2})/2$ referred to as “$\pm$zFo-N” or
$S_{z}^{\rm(ss)}=\pm(N_{1}-N_{2})/2$ designated by “$\pm$zFi-N”; see the table
in the right side of Fig. 1. Once again, the $\pm$zFi-N states cross over into
$z$-antiferromagnetic normal states, “$\pm$zaF-N”, for $N_{1}=N_{2}$.

The superradiant thresholds are obtained by setting $S_{l,x}^{\rm(ss)}=0$ in
Eq. (4), yielding $\lambda=\lambda_{c}^{(\rm xFo/xFi)}$. The fixed points
$\pm$xFo-SR ($\pm$xFi-SR) only emerge beyond the threshold $\lambda_{c}^{(\rm
xFo)}$ [$\lambda_{c}^{(\rm xFi)}$]. Note, however, that the xFo-SR threshold
diverges for $N_{2}\to N_{1}$ as $\alpha_{\rm
ss}\propto\Delta{S}_{x}^{\rm(ss)}\to 0$.

We now turn our attention to the linear stability of the fixed points of the
system obtained above. To this end, we write $a=a_{\rm ss}+\delta a$ and
$S_{l,\alpha}=S_{l,\alpha}^{\rm(ss)}+\delta S_{l,\alpha}$ in the mean-field
equations of motion (Unconventional Dicke model: Multistabilities and
nonequilibrium dynamics) and subsequently linearize them to obtain
$\partial_{t}\delta\mathbf{X}=\mathbf{J}\delta\mathbf{X}$, where
$\delta\mathbf{X}$ is a vector of fluctuations $\delta a$ and $\delta
S_{l,\alpha}$, and $\mathbf{J}$ the Jacobian matrix given explicitly in SM Miv
. A fixed point is stable provided all eigenvalues of the Jacobian matrix for
that given fixed point have negative real parts Roussel (2019). We analyze
eigenvalues of the Jacobian matrix for all the fixed points. In particular, we
find that all $\pm$xFo-SR and $\pm$xFi-SR states are stable in their entire
corresponding parameter regimes, i.e., beyond the $\lambda_{c}^{(\rm
xFo/xFi)}$ thresholds, respectively. On the other hand, the trivial fixed
points $-$zFo-N and $-$zFi-N are solely stable below the $\lambda_{c}^{(\rm
xFi)}$ and $\lambda_{c}^{(\rm xFo)}$ thresholds, respectively, and lose their
stability beyond these thresholds. The other two trivial fixed points +zFo-N
and +zFi-N are always unstable. This means the system undergoes a
supercritical pitchfork bifurcation from $-$zFo-N into $\pm$xFi-SR at the
threshold $\lambda_{c}^{(\rm xFi)}$, and from $-$zFi-N into $\pm$xFo-SR at the
threshold $\lambda_{c}^{(\rm xFo)}$.

The steady-state phase diagram of the system in the parameter plane of
$\\{N_{2}/N_{1},\lambda/\kappa\\}$ is shown in Fig. 1. The two white curves
are the analytical boundaries for the xFo-SR and xFi-SR transitions, obtained
from Eq. (5). The total number of fixed points in each parameter regime is
indicated in the phase diagram. However, only the stable fixed points in each
parameter regime are indicated explicitly. The phase diagram features regions
of multistability Zhou  _et al._ (2009); Ali  _et al._ (2022), in particular,
a region with multiple coexistent superradiant phases.

Figure 2: Steady-state behavior of the system across the xFi-SR and xFo-SR
thresholds: the $\\{x,z\\}$-components $S_{x,z}^{(\rm ss)}$ of the total spin
(a,b), the $x$-component $\Delta S_{x}^{(\rm ss)}$ of the staggered spin (c),
and the energy $E_{0}$ of the system (d) as a function of the atom-field
coupling $\lambda$ at a fixed $N_{2}/N_{1}=0.3$ for the fixed points $-$zFi-N
(orange), $-$zFo-N (green), $\pm$xFo-SR (red), and $\pm$xFi-SR (blue). The $x$
component of the total and staggered spin both acquire non-zero values in the
$\pm$xFi-SR and $\pm$xFo-SR states, thus implying that the system can undergo
two independent supercritical pitchfork bifurcations. The $z$-component
$S_{z}^{(\rm ss)}$ of the total spin and the energy $E_{0}$ change
continuously in the transition from the $-$zFo-N to $\pm$xFi-SR states and
from the $-$zFi-N to $\pm$xFo-SR states. The dashed green and orange lines
indicate that the corresponding fixed points are unstable for the given
coupling $\lambda$. The other parameters are the same as Fig. 1.

The steady-state behavior of the $\\{x,z\\}$-components of the total spin
$S_{x,z}^{(\rm ss)}$, the $x$-component of the staggered spin
$\Delta{S}_{x}^{(\rm ss)}$, and the energy of the system
$E_{0}=\langle\hat{H}_{\rm nsD}\rangle$ as a function of the atom-field
coupling $\lambda$ at a fixed $N_{2}/N_{1}=0.3$ are depicted in Fig. 2 for the
fixed points $-$zFo-N, $-$zFi-N, $\pm$xFo-SR, and $\pm$xFi-SR. As expected,
$\Delta{S}_{x}^{(\rm ss)}$ acquires nonzero values across the
$\lambda_{c}^{(\rm xFi/xFo)}$ thresholds in the $\pm$xFi-SR (blue curves) and
$\pm$xFo-SR (red curves) states, respectively; see Fig. 2(c). The two branches
in each case correspond to the two possible $\mathbb{Z}_{2}$ solutions.
Accordingly, the field amplitude $\alpha_{\rm ss}\propto\Delta{S}_{x}^{(\rm
ss)}$ also grows from zero and the system enters superradiant phases, implying
that $\Delta{S}_{x}^{(\rm ss)}$ can be identified as the order parameter of
the system. Note that $\Delta{S}_{x}^{(\rm ss)}$ is zero in both $-$zFo-N and
$-$zFi-N states [green and orange lines, respectively; obstructed somewhat by
each other in Fig. 2(c)]. The $x$ component of the total spin, $S_{x}^{(\rm
ss)}$, exhibits a similar behavior, though it grows much faster in the
$\pm$xFo-SR states compared to the $\pm$xFi-SR states; see Fig. 2(a). From
Fig. 2(b) and (d), one sees that the $z$-component $S_{z}^{(\rm ss)}$ of the
total spin and the energy $E_{0}$ change continuously from the $-$zFo-N to
$\pm$xFi-SR states and from the $-$zFi-N into $\pm$xFo-SR states, signaling
second-order phase transitions. This is consistent with Fig. 2(c) and the
linear-stability analysis predicting supercritical pitchfork bifurcations. The
behavior of $S_{x,z}^{(\rm ss)}$, $\Delta{S}_{x}^{(\rm ss)}$, and $E_{0}$ in
the full parameter space of $\\{N_{2}/N_{1},\lambda/\kappa\\}$ are given in SM
Miv .

From the energy point of view, the fixed points $-$zFo-N and $\pm$xFi-SR in
their corresponding stability regimes are the lowest-energy statets, implying
that the other equilibrium points are metastable states. However, in open
dynamical systems, dynamical stability is more important than the minimal
energy condition Bhaseen  _et al._ (2012). Therefore, there is an essential
distinction between the dynamics of an open system with $\kappa\neq 0$ and its
equilibrium counterpart with $\kappa=0$.

We now examine the semiclassical dynamics of the system by numerically
integrating the mean-field equations of motion (Unconventional Dicke model:
Multistabilities and nonequilibrium dynamics). As expected, the stable steady
states are the attractors of the long-time dynamics, while a small
perturbation in the unstable steady states leads to phase-space trajectories
being repelled from these points into one of the stable fixed points. We study
these by adding a small perturbation $\delta a$ to the steady-state field
amplitude $a_{\rm ss}$ in each fixed point and looking at the ensuing
dynamics. In particular, this reveals that both $\pm$zFo-N ($\pm$zFi-N) are
unstable towards one of the parity-symmetric pair $\pm$xFi-SR ($\pm$xFo-SR)
above the threshold $\lambda_{c}^{(\rm xFi)}$ [$\lambda_{c}^{(\rm xFo)}$]; see
SM for these types of phase-space dynamics. While below the threshold, +zFo-N
(+zFi-N) evolves to the corresponding low-energy state $-$zFo-N ($-$zFi-N).

Note that due to the multistability, the long-time dynamics depend crucially
on the initial state and can exhibit intriguing features. In particular, in
some parameter regimes we find initial states where the following phase-space
dynamics from them do not lead to any of the stable fixed points, rather
exhibit oscillatory behaviors owing to competing fixed points Thompson (1984).
This is especially interesting in the parameter regime where all $\pm$xFo-SR
and $\pm$xFi-SR are the stable fixed points of the system. Starting from the
initial state $\mathbf{S}_{1}=(1/\sqrt{2},1/\sqrt{2},0)N_{1}/2$ and
$\mathbf{S}_{2}=(1/\sqrt{2},0,-1/\sqrt{2})N_{2}/2$ with a small field
$a=\delta a$, the ensuing trajectories in the Bloch spheres of
$\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ encircle the two fixed points $-$xFo-SR
and $-$xFi-SR and exhibit limit cycles as shown in Fig. 3(b) Kosior  _et al._
(2022). That said, the corresponding trajectory in the phase space of the
cavity-field amplitude $a$ as shown in Fig. 3(a) exhibits an inward spiral
behavior toward an emergent focus lying between the two fixed points $-$xFo-SR
and $-$xFi-SR. The dynamics of the system in the Bloch spheres of the total
$\mathbf{S}$ and staggered $\Delta\mathbf{S}$ spins with radius
$(N_{1}+N_{2})/2$ are displayed in the insets of Fig. 3(a). Note that although
the dynamics in the Bloch spheres $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ are
strictly restricted to the Bloch surfaces, in the $\mathbf{S}$ and
$\Delta\mathbf{S}$ phase spaces the dynamics move inside the Bloch spheres,
confirming that the total spin $\mathbf{S}$ is not conserved.

Figure 3: The nonequilibrium dynamics of the system. The system is prepared in
the initial state $\mathbf{S}_{1}=(1/\sqrt{2},1/\sqrt{2},0)N_{1}/2$,
$\mathbf{S}_{2}=(1/\sqrt{2},0,-1/\sqrt{2})N_{2}/2$, and a small $a=\delta a$.
The ensuing dynamics of the system are displayed in the phase space of the the
cavity-field amplitude $a$ (a) and the $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$
spin Bloch spheres (b). The color gradient from light blue to darker blue in
trajectories indicates schematically the arrow of the time evolution. The
trajectories in the spin Bloch spheres exhibit limit-cycle behaviors,
enclosing the two fixed points $-$xFo-SR and $-$xFi-SR. On the other hand, the
corresponding trajectory in the field-amplitude space is attracted to an
emergent focus lying between the $-$xFo-SR and $-$xFi-SR fixed points. The
insets in panel (a) display the dynamics in the the Bloch spheres of the total
$\mathbf{S}$ and staggered $\Delta\mathbf{S}$ spins. The radii of the Bloch
spheres are indicated in each case, where $\mathbf{S}$ and $\Delta\mathbf{S}$
have the same radius $(N_{1}+N_{2})/2$. Here $\lambda/\kappa=2$ and
$N_{2}/N_{1}=0.3$. The other parameters are the same as Fig. 1.

_Quantum description._ —Finally, we address briefly the full quantum
description of the system via the master equation for the density-matrix
operator, $\hat{\dot{\rho}}=i[\hat{\rho},\hat{H}_{\rm
nsD}]+\hat{\mathcal{L}}[\hat{\rho}]$. In the Born-Markov approximation, the
Liouvillian can be expressed in the Lindblad form as
$\hat{\mathcal{L}}[\hat{\rho}]=\kappa(2\hat{a}\hat{\rho}\hat{a}^{\dagger}-\hat{a}^{\dagger}\hat{a}\hat{\rho}-\hat{\rho}\hat{a}^{\dagger}\hat{a})$,
where we have ignored the decay of the atomic excited states.

We perform small-scale, full-quantum calculations with $N_{1}=4$ and $N_{2}=3$
and find that indeed depending on the initial state, $\pm$xFo-SR or $\pm$xFi-
SR can be the final state of the long-time quantum dynamics of the system in
appropriate parameter regimes. Remarkably, we also find initial states and
parameter regimes where the ensuing quantum dynamics from them lead to a state
which is a superposition of all $\pm$xFo-SR and $\pm$xFi-SR. This can be best
seen in the $Q$-representation of the cavity field Sandner  _et al._ (2015);
Halati  _et al._ (2020). The $Q$-function is presented in Fig. 4 and comprised
of four, partially disjoint lobes, an indication of the coexistence of
multiple superradiant states, i.e., $\pm$xFo-SR and $\pm$xFi-SR [cf. Fig.
3(a)].

Figure 4: The Husimi $Q$ representation of the cavity field after long-time,
quantum dynamics. It consists of four bubbles, an indication of the
superposition of four $\pm$xFo-SR and $\pm$xFi-SR superradiant states. The
parameters are $N_{1}=4$, $N_{2}=3$, and $\lambda=1.01\kappa$, with the rest
being the same as Fig. 1.

_Conclusions._ —We have introduced a variant of the Dicke model, which
exhibits a wealth of steady-state and non-stationary phenomena. Our model can
readily be implemented in state-of-the-art experiments Gupta  _et al._ (2007);
Muniz  _et al._ (2020); Bohr  _et al._ (2023) and opens a new avenue for
studying various nonequilibrium magnetic ordering Davis  _et al._ (2019) and
dynamical phenomena Dogra  _et al._ (2019); Kongkhambut  _et al._ (2022) in
cavity-QED experimental setups.

###### Acknowledgements.

We acknowledge inspiring discussions with Natalia Masalaeva, Karol Gietka,
Arkadiusz Kosior, and Helmut Ritsch. F. M. is supported financially by the
Stand-alone project P 35891-N of the Austrian Science Fund (FWF), Tyrolean
Science Promotion Fund (TWF), and the ESQ-Discovery Grant of the Austrian
Academy of Sciences (ÖAW).

## References

  * Hepp and Lieb (1973a) Klaus Hepp and Elliott H Lieb, “On the superradiant phase transition for molecules in a quantized radiation field: the dicke maser model,” Annals of Physics 76, 360–404 (1973a).
  * Hepp and Lieb (1973b) Klaus Hepp and Elliott H. Lieb, “Equilibrium statistical mechanics of matter interacting with the quantized radiation field,” Phys. Rev. A 8, 2517–2525 (1973b).
  * Wang and Hioe (1973) Y. K. Wang and F. T. Hioe, “Phase transition in the dicke model of superradiance,” Phys. Rev. A 7, 831–836 (1973).
  * Garraway (2011) Barry M. Garraway, “The dicke model in quantum optics: Dicke model revisited,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 369, 1137–1155 (2011).
  * Kirton  _et al._ (2018) Peter Kirton, Mor M. Roses, Jonathan Keeling, and Emanuele G. Dalla Torre, “Introduction to the dicke model: From equilibrium to nonequilibrium, and vice versa,” Advanced Quantum Technologies 2, 1800043 (2018).
  * Baumann  _et al._ (2010) Kristian Baumann, Christine Guerlin, Ferdinand Brennecke, and Tilman Esslinger, “Dicke quantum phase transition with a superfluid gas in an optical cavity,” Nature 464, 1301–1306 (2010).
  * Klinder  _et al._ (2015) Jens Klinder, Hans Keßler, Matthias Wolke, Ludwig Mathey, and Andreas Hemmerich, “Dynamical phase transition in the open dicke model,” Proceedings of the National Academy of Sciences 112, 3290–3295 (2015).
  * Kollár  _et al._ (2017) Alicia J. Kollár, Alexander T. Papageorge, Varun D. Vaidya, Yudan Guo, Jonathan Keeling, and Benjamin L. Lev, “Supermode-density-wave-polariton condensation with a bose–einstein condensate in a multimode cavity,” Nature Communications 8 (2017), 10.1038/ncomms14386.
  * Zhang  _et al._ (2021) Xiaotian Zhang, Yu Chen, Zemao Wu, Juan Wang, Jijie Fan, Shujin Deng, and Haibin Wu, “Observation of a superradiant quantum phase transition in an intracavity degenerate fermi gas,” Science 373, 1359–1362 (2021).
  * Helson  _et al._ (2023) Victor Helson, Timo Zwettler, Farokh Mivehvar, Elvia Colella, Kevin Roux, Hideki Konishi, Helmut Ritsch, and Jean-Philippe Brantut, “Density-wave ordering in a unitary fermi gas with photon-mediated interactions,” Nature 618, 716–720 (2023).
  * Zhiqiang  _et al._ (2017) Zhang Zhiqiang, Chern Hui Lee, Ravi Kumar, K. J. Arnold, Stuart J. Masson, A. S. Parkins, and M. D. Barrett, “Nonequilibrium phase transition in a spin-1 dicke model,” Optica 4, 424 (2017).
  * Bialynicki-Birula and Rza¸żewski (1979) Iwo Bialynicki-Birula and Kazimierz Rza¸żewski, “No-go theorem concerning the superradiant phase transition in atomic systems,” Phys. Rev. A 19, 301–303 (1979).
  * Nataf and Ciuti (2010) Pierre Nataf and Cristiano Ciuti, “No-go theorem for superradiant quantum phase transitions in cavity QED and counter-example in circuit QED,” Nature Communications 1 (2010), 10.1038/ncomms1069.
  * Safaei  _et al._ (2013) S Safaei, Ö E Müstecaplıoğlu, and B Tanatar, “Raman superradiance and spin lattice of ultracold atoms in optical cavities,” New Journal of Physics 15, 083037 (2013).
  * Mivehvar  _et al._ (2017) Farokh Mivehvar, Francesco Piazza, and Helmut Ritsch, “Disorder-driven density and spin self-ordering of a bose-einstein condensate in a cavity,” Phys. Rev. Lett. 119, 063602 (2017).
  * Mivehvar  _et al._ (2019) Farokh Mivehvar, Helmut Ritsch, and Francesco Piazza, “Cavity-quantum-electrodynamical toolbox for quantum magnetism,” Phys. Rev. Lett. 122, 113603 (2019).
  * Masalaeva  _et al._ (2021) Natalia Masalaeva, Wolfgang Niedenzu, Farokh Mivehvar, and Helmut Ritsch, “Spin and density self-ordering in dynamic polarization gradients fields,” Phys. Rev. Res. 3, 013173 (2021).
  * Landini  _et al._ (2018) M. Landini, N. Dogra, K. Kroeger, L. Hruby, T. Donner, and T. Esslinger, “Formation of a spin texture in a quantum gas coupled to a cavity,” Phys. Rev. Lett. 120, 223602 (2018).
  * Kroeze  _et al._ (2018) Ronen M. Kroeze, Yudan Guo, Varun D. Vaidya, Jonathan Keeling, and Benjamin L. Lev, “Spinor self-ordering of a quantum gas in a cavity,” Phys. Rev. Lett. 121, 163601 (2018).
  * Kroeze  _et al._ (2019) Ronen M. Kroeze, Yudan Guo, and Benjamin L. Lev, “Dynamical spin-orbit coupling of a quantum gas,” Phys. Rev. Lett. 123, 160404 (2019).
  * Chang  _et al._ (2018) D. E. Chang, J. S. Douglas, A. González-Tudela, C.-L. Hung, and H. J. Kimble, “Colloquium: Quantum matter built from nanoscopic lattices of atoms and photons,” Rev. Mod. Phys. 90, 031002 (2018).
  * Sheremet  _et al._ (2023) Alexandra S. Sheremet, Mihail I. Petrov, Ivan V. Iorsh, Alexander V. Poshakinskiy, and Alexander N. Poddubny, “Waveguide quantum electrodynamics: Collective radiance and photon-photon correlations,” Rev. Mod. Phys. 95, 015002 (2023).
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef, “Cavity quantum materials,” Applied Physics Reviews 9, 011312 (2022).
  * Mivehvar  _et al._ (2021) Farokh Mivehvar, Francesco Piazza, Tobias Donner, and Helmut Ritsch, “Cavity QED with quantum gases: new paradigms in many-body physics,” Advances in Physics 70, 1–153 (2021).
  * Hotter  _et al._ (2023) Christoph Hotter, Laurin Ostermann, and Helmut Ritsch, “Cavity sub- and superradiance for transversely driven atomic ensembles,” Phys. Rev. Res. 5, 013056 (2023).
  * (26) See supplemental material for the details of mapping the combined density-spin self-ordering of two-component quantum gases into the non-standard Dicke model, derivation of the steady-state solutions and their linear-stability analyses, and the dependence of the long-time dynamics on the initial state .
  * Keeling  _et al._ (2010) J. Keeling, M. J. Bhaseen, and B. D. Simons, “Collective dynamics of bose-einstein condensates in optical cavities,” Phys. Rev. Lett. 105, 043001 (2010).
  * Roussel (2019) Marc R Roussel,  _Nonlinear Dynamics_ (Morgan & Claypool Publishers, 2019).
  * Zhou  _et al._ (2009) Lu Zhou, Han Pu, Hong Y. Ling, and Weiping Zhang, “Cavity-mediated strong matter wave bistability in a spin-1 condensate,” Phys. Rev. Lett. 103, 160403 (2009).
  * Ali  _et al._ (2022) Abid Ali, Farhan Saif, and Hiroki Saito, “Phase separation and multistability of a two-component bose-einstein condensate in an optical cavity,” Phys. Rev. A 105, 063318 (2022).
  * Bhaseen  _et al._ (2012) M. J. Bhaseen, J. Mayoh, B. D. Simons, and J. Keeling, “Dynamics of nonequilibrium dicke models,” Phys. Rev. A 85, 013817 (2012).
  * Thompson (1984) J.M.T. Thompson, “An introduction to nonlinear dynamics,” Applied Mathematical Modelling 8, 157–168 (1984).
  * Kosior  _et al._ (2022) Arkadiusz Kosior, Helmut Ritsch, and Farokh Mivehvar, “Nonequilibrium phases of ultracold bosons with cavity-induced dynamic gauge fields,” Pre-print: arXiv:2208.04602 (2022).
  * Sandner  _et al._ (2015) R. M. Sandner, W. Niedenzu, F. Piazza, and H. Ritsch, “Self-ordered stationary states of driven quantum degenerate gases in optical resonators,” EPL (Europhysics Letters) 111, 53001 (2015).
  * Halati  _et al._ (2020) Catalin-Mihai Halati, Ameneh Sheikhan, Helmut Ritsch, and Corinna Kollath, “Numerically exact treatment of many-body self-organization in a cavity,” Phys. Rev. Lett. 125, 093604 (2020).
  * Gupta  _et al._ (2007) Subhadeep Gupta, Kevin L. Moore, Kater W. Murch, and Dan M. Stamper-Kurn, “Cavity nonlinear optics at low photon numbers from collective atomic motion,” Phys. Rev. Lett. 99, 213601 (2007).
  * Muniz  _et al._ (2020) Juan A. Muniz, Diego Barberena, Robert J. Lewis-Swan, Dylan J. Young, Julia R. K. Cline, Ana Maria Rey, and James K. Thompson, “Exploring dynamical phase transitions with cold atoms in an optical cavity,” Nature 580, 602–607 (2020).
  * Bohr  _et al._ (2023) Eliot Bohr, Sofus L. Kristensen, Christoph Hotter, Stefan Alaric Schäffer, Julian Robinson-Tait, Jan W. Thomsen, Tanya Zelevinsky, Helmut Ritsch, and Jörg Helge Müller, “Collectively enhanced ramsey readout by cavity sub- to superradiant transition,” pre-print: arXiv:2306.12544 (2023).
  * Davis  _et al._ (2019) Emily J. Davis, Gregory Bentsen, Lukas Homeier, Tracy Li, and Monika H. Schleier-Smith, “Photon-mediated spin-exchange dynamics of spin-1 atoms,” Phys. Rev. Lett. 122, 010405 (2019).
  * Dogra  _et al._ (2019) Nishant Dogra, Manuele Landini, Katrin Kroeger, Lorenz Hruby, Tobias Donner, and Tilman Esslinger, “Dissipation-induced structural instability and chiral dynamics in a quantum gas,” Science 366, 1496–1499 (2019).
  * Kongkhambut  _et al._ (2022) Phatthamon Kongkhambut, Jim Skulte, Ludwig Mathey, Jayson G. Cosme, Andreas Hemmerich, and Hans Keßler, “Observation of a continuous time crystal,” Science 377, 670–673 (2022).

##  I Supplemental Material

In the supplemental material, we present the details of the connection of the
non-standard Dicke model to the combined density-spin self-ordering in
transversely-driven two-component quantum gases inside cavities, the
analytical derivation of the steady states, their linear-stability analysis,
and the dependence of the semi-classical dynamics of the system on the initial
state.

###  I.1 The combined density-spin self-ordering in quantum gases: The non-
standard Dicke model

An effective two-component Bose-Einstein condensate (BEC) coupled to a single
mode of a standing-wave cavity and driven in a transverse direction by two
pump lasers as depicted in Fig. S1(a) is described by the Hamiltonian ($\hbar$
is set to one)
$\hat{H}=-\Delta_{c}\hat{a}^{\dagger}\hat{a}+\int\hat{\Psi}^{\dagger}(x)\hat{\mathcal{H}}_{0}\hat{\Psi}(x)dx$,
where $\Delta_{c}$ is the pump-cavity detuning, and $\hat{a}$ and
$\hat{\Psi}=(\hat{\psi}_{\uparrow},\hat{\psi}_{\downarrow})^{\top}$ are the
photonic and two-component atomic bosonic annihilation field operators,
respectively; see also Refs. Mivehvar  _et al._ (2017, 2019) for more details.
The single-particle Hamiltonian density reads,

$\displaystyle\hat{\mathcal{H}}_{0}=\left[\frac{p_{x}^{2}}{2M}+\frac{1}{2}\hat{U}_{\uparrow}(x)+\frac{1}{2}\hat{U}_{\downarrow}(x)\right]\mathbb{I}_{2\times
2}+\frac{1}{2}\begin{pmatrix}\delta+\hat{U}_{\uparrow}(x)-\hat{U}_{\downarrow}(x)&2\hat{\Omega}_{\rm
R}(x)\\\ 2\hat{\Omega}^{\dagger}_{\rm
R}(x)&-\delta-\hat{U}_{\uparrow}(x)+\hat{U}_{\downarrow}(x)\end{pmatrix},$
(S1)

with the potential operators
$\hat{U}_{\uparrow}(x)=U_{0\uparrow}\hat{a}^{\dagger}\hat{a}\cos^{2}(k_{c}x)$
and
$\hat{U}_{\downarrow}(x)=U_{0\downarrow}\hat{a}^{\dagger}\hat{a}\cos^{2}(k_{c}x)$,
and the two-photon Raman-Rabi coupling operator $\hat{\Omega}_{\rm
R}(x)=\lambda(\hat{a}^{\dagger}+\hat{a})\cos(k_{c}x)$. Here, $M$ is the atomic
mass, $p_{x}$ the atomic center-of-mass momentum operator, $\delta$ the
effective detuning between the atomic ground states, and $\mathbb{I}_{2\times
2}$ the identity matrix in the atomic ground-state space. The parameters
$U_{0\uparrow}$, $U_{0\downarrow}$, and $\lambda$ are proportional to atom-
cavity and atom-pump coupling strengths as well as the pump-atom detuning
Mivehvar  _et al._ (2017, 2019).

Figure S1: A quantum gas inside a cavity. (a) A BEC with two low-lying
internal states $\\{\left|\uparrow\right>,\left|\downarrow\right>\\}$ is
coupled to a single mode of a linear cavity and further driven transversely by
two pump lasers in a double $\Lambda$ scheme. The inset shows the details of
atomic-field couplings. The atomic excited states
$\\{\left|1\right>,\left|2\right>\\}$ are very far detuned and do not enter
into the low-energy description of the system considered here. (b) The mean-
field atomic spin texture along the cavity axis:
$s_{x}(x)=2\text{Re}[\psi_{\uparrow}^{*}(x)\psi_{\downarrow}(x)]$,
$s_{y}(x)=2\text{Im}[\psi_{\uparrow}^{*}(x)\psi_{\downarrow}(x)]$, and
$s_{z}(x)=|\psi_{\uparrow}(x)|^{2}-|\psi_{\downarrow}(x)|^{2}$. The
$x$-component of the spin $s_{x}(x)$ points in opposite directions at even
$k_{c}x_{m}=2m\pi$ and odd $k_{c}x_{m}=(2m+1)\pi$ lattice sites, where
$m=0,\pm 1,\pm 2,\cdots$. However, $s_{z}(x)$ points in the negative direction
in both sub-lattices and $s_{y}(x)$ is zero everywhere. The parameters are set
to
$(\Delta_{c},\kappa,\delta,NU_{0\downarrow},NU_{0\uparrow},\sqrt{N}\lambda)=(-8,4,28,-4,-4,-12)\omega_{r}$,
with $\omega_{r}=k_{c}^{2}/2M$ being the recoil frequency.

By defining the local pseudospin operator
$\hat{\mathbf{s}}(x)=\hat{\Psi}^{\dagger}(x)\boldsymbol{\sigma}\hat{\Psi}(x)$
($\boldsymbol{\sigma}$ is the vector of Pauli matrices), the Hamiltonian of
the system can be recast as

$\displaystyle\hat{H}=-\Delta_{c}\hat{a}^{\dagger}\hat{a}$
$\displaystyle+\int\hat{\Psi}^{\dagger}(x)\left[\frac{p_{x}^{2}}{2M}+\frac{1}{2}\hat{U}_{\uparrow}(x)+\frac{1}{2}\hat{U}_{\downarrow}(x)\right]\mathbb{I}_{2\times
2}\hat{\Psi}(x)dx$
$\displaystyle+\int\left\\{\frac{1}{2}\left[\delta+\hat{U}_{\uparrow}(x)-\hat{U}_{\downarrow}(x)\right]\hat{s}_{z}(x)+\hat{\Omega}_{\rm
R}(x)\hat{s}_{x}(x)\right\\}dx.$ (S2)

The first term in Eq. (I.1) is the bare cavity Hamiltonian, the second term
describes the self-ordering of the total atomic density, and the last term is
associated with the spin self-ordering.

We numerically solve Eq. (I.1) in the mean-field regime $\hat{a}\to
a=\langle\hat{a}\rangle$ and
$\hat{\Psi}(x)\to\Psi(x)=\langle\hat{\Psi}(x)\rangle$ [and consequently
$\hat{\mathbf{s}}(x)\to\mathbf{s}(x)=\langle\hat{\mathbf{s}}(x)\rangle$]. For
small $\lambda$ the BEC is uniform and the cavity field is in the vacuum
state, while for strong enough $\lambda$ the system enters the superradiant
phase exhibiting density and spin self-ordering; see Fig. S1(b). The self-
ordered density and spin texture are localized at $x_{j}=j\pi/k_{c}$ with
$2\pi/k_{c}$ being the wavelength of the cavity mode and $j=0,\pm 1,\pm
2,\cdots$. Therefore, it is justified to discretize the system into a lattice
with lattice sites $x_{j}$, specially deep in the self-ordered phase. One can
see that $s_{j,z}\equiv s_{z}(x_{j})$ always point down (for the chosen
positive $\delta$), while $s_{j,x}\equiv s_{x}(x_{j})$ alternates its
direction. Note also that $s_{j,y}\equiv s_{y}(x_{j})$ is zero everywhere.

In order to explain these results, we consider the spin part of the
Hamiltonian (I.1) (i.e., the last term) and re-write it in the discrete space,

$\displaystyle\hat{H}_{\rm spin}=\sum_{j=0,\pm 1,\pm
2,\cdots}\left(\hat{h}_{j,z}\hat{s}_{j,z}+\hat{h}_{j,x}\hat{s}_{j,x}\right)=\hat{h}_{0,z}\sum_{j=0,\pm
1,\pm 2,\cdots}\hat{s}_{j,z}+\hat{h}_{0,x}\left(\sum_{j=0,\pm
2,\cdots}\hat{s}_{j,x}-\sum_{j=\pm 1,\pm 3,\cdots}\hat{s}_{j,x}\right),$ (S3)

where
$\hat{h}_{j,z}=[\delta+(U_{0\uparrow}-U_{0\downarrow})\hat{a}^{\dagger}\hat{a}\cos^{2}(k_{c}x_{j})]/2$
and $\hat{h}_{j,x}=\lambda(\hat{a}^{\dagger}+\hat{a})\cos(k_{c}x_{j})$. In the
last equality, we have made use of the fact that $\hat{h}_{j,z}$ is site
independent, while $\hat{h}_{j,x}$ has opposite values in even and odd sites.

Defining two sub-lattices consisting only of even and odd sites, respectively,
and assuming that $U_{0\uparrow}=U_{0\downarrow}$, we obtain

$\displaystyle\hat{H}_{\rm
spin}=\frac{\delta}{2}\hat{S}_{z}+\lambda(\hat{a}^{\dagger}+\hat{a})\left(\hat{S}_{1,x}-\hat{S}_{2,x}\right),$
(S4)

where $\hat{S}_{z}=\sum_{j=0,\pm 1,\pm 2,\cdots}\hat{s}_{j,z}$ is the $z$
component of the total spin in both sub-lattices, and
$\hat{S}_{1,x}=\sum_{j=0,\pm 2,\cdots}\hat{s}_{j,x}$ and
$\hat{S}_{2,x}=\sum_{j=\pm 1,\pm 3,\cdots}\hat{s}_{j,x}$ are the $x$
components of the spin in each sub-lattices. Adding the bare cavity
Hamiltonian into $\hat{H}_{\rm spin}$ (S4) and replacing
$\delta/2\to\omega_{a}$ and $-\Delta_{c}\to\omega_{c}$, we arrive at the non-
standard Dicke Hamiltonian $\hat{H}_{\rm nsD}$, given in Eq. (2) in the main
text. We note that in the deep self-ordered regime the atomic density is
almost frozen and its Hamiltonian [i.e., the second term in Eq. (I.1)] has
been, therefore, omitted. The non-standard Dicke Hamiltonian hence describes
the physics of the combined density-spin self-ordering of transversely-driven
two-component Bose gases inside single-mode cavities deep in the superradiant
phase. Note that here $|\mathbf{S}_{1}|=|\mathbf{S}_{2}|$. However, this
restriction can be relaxed by adding further elements into the model.

###  I.2 Steady-state solutions

In this section we derive the steady-state solutions in more details. By
setting $\dot{a}=0$ and $\dot{S}_{l,\alpha}=0$ in Eq. (Unconventional Dicke
model: Multistabilities and nonequilibrium dynamics) in the main text, one
obtains

$\displaystyle a_{\rm ss}$
$\displaystyle=\frac{\lambda}{-\omega_{c}+i\kappa}\left[S_{1,x}^{\rm(ss)}-S_{2,x}^{\rm(ss)}\right],$
(S5a) $\displaystyle S_{l,x}^{\rm(ss)}$
$\displaystyle=(-1)^{l+1}\frac{\lambda}{\omega_{a}}(a_{\rm ss}+a_{\rm
ss}^{*})S_{l,z}^{\rm(ss)}\quad\text{with}\quad l=1,2.$ (S5b)

Dividing the two equations in Eq. (S5b) leads to a relation between
$S_{1,x/z}^{\rm(ss)}$ and $S_{2,x/z}^{\rm(ss)}$,

$\displaystyle\frac{S_{1,x}^{\rm(ss)}}{S_{2,x}^{\rm(ss)}}=-\frac{S_{1,z}^{\rm(ss)}}{S_{2,z}^{\rm(ss)}}.$
(S6)

Here, $S_{1,z}^{\rm(ss)}$ and $S_{2,z}^{\rm(ss)}$ can have the same or
opposite signs, which lead to two distinct classes of nontrivial solutions. In
particular, when $S_{1,z}^{\rm(ss)}$ and $S_{2,z}^{\rm(ss)}$ have the same
sign (opposite signs), $S_{1,x}^{\rm(ss)}$ and $S_{2,x}^{\rm(ss)}$ must have
opposite signs (the same sign), in order for Eq. (S6) to be satisfied. These
correspond, respectively, to the $\pm$xFi-SR and $\pm$xFo-SR ordering.

By noting that in the steady states
$|S_{l,z}^{\rm(ss)}|=\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$, Eq. (S6)
can be expressed in terms of only $S_{l,x}^{\rm(ss)}$ as,

$\displaystyle\frac{S_{1,x}^{\rm(ss)}}{S_{2,x}^{\rm(ss)}}=\pm\frac{\sqrt{\left(\frac{N_{1}}{2}\right)^{2}-[S_{1,x}^{\rm(ss)}]^{2}}}{\sqrt{\left(\frac{N_{2}}{2}\right)^{2}-[S_{2,x}^{\rm(ss)}]^{2}}}.$
(S7)

Equation (S7) can be simplified to yield the relation
$S_{1,x}^{\rm(ss)}=\pm(N_{1}/N_{2})S_{2,x}^{\rm(ss)}$, where the upper plus
(lower minus) sign corresponds to the xFo-SR (xFi-SR). Using this relation,
the steady-state field amplitude $a_{\rm ss}$ [Eq. (S5a)] can be expressed in
terms of one of the ensemble spins [say, $S_{1,x}^{\rm(ss)}$]. Substituting
this steady-sate field amplitude back in the corresponding equation for the
spin $S_{1,x}^{\rm(ss)}$ [Eq. (S5b)], one obtains

$\displaystyle 1$
$\displaystyle=-\frac{2\omega_{c}\lambda^{2}}{\omega_{a}(\omega_{c}^{2}+\kappa^{2})}\left(1\mp\frac{N_{2}}{N_{1}}\right)\sqrt{\left(\frac{N_{1}}{2}\right)^{2}-\left[S_{1,x}^{\rm(ss)}\right]^{2}}.$
(S8)

This equation can readily be solved for $S_{1,x}^{\rm(ss)}$ and subsequently
$S_{2,x}^{\rm(ss)}$,

$\displaystyle
S_{1,x}^{\rm(ss)}=\pm\sqrt{\left(\frac{N_{1}}{2}\right)^{2}-\left[\frac{\omega_{a}(\omega_{c}^{2}+\kappa^{2})}{2\left(1\mp\frac{N_{2}}{N_{1}}\right)\omega_{c}\lambda^{2}}\right]^{2}},$
$\displaystyle
S_{2,x}^{\rm(ss)}=\pm\sqrt{\left(\frac{N_{2}}{2}\right)^{2}-\left[\frac{\omega_{a}(\omega_{c}^{2}+\kappa^{2})}{2\left(\frac{N_{1}}{N_{2}}\mp
1\right)\omega_{c}\lambda^{2}}\right]^{2}},$ (S9)

yielding, after defining the critical couplings $\lambda_{c}^{(\rm
xFo/xFi)}=\sqrt{\omega_{a}(\omega_{c}^{2}+\kappa^{2})/(N_{1}\mp
N_{2})\omega_{c}}$, Eq. (4) in the main text.

After obtaining $S_{l,x}^{\rm(ss)}$, the steady-state $S_{l,z}^{\rm(ss)}$ can
be obtained from the normalization constraints,
$S_{l,z}^{\rm(ss)}=\pm\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$. As
discussed above, $S_{1,z}^{\rm(ss)}$ and $S_{2,z}^{\rm(ss)}$ have the same
sign (opposite signs) in the $\pm$xFi-SR ($\pm$xFo-SR). In order to determine
the exact signs, however, we numerically solve the steady-state equations. We
find that in the $\pm$xFi-SR states the minus sign is always picked up for
both $S_{l,z}^{\rm(ss)}$, that is,
$S_{l,z}^{\rm(ss)}=-\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$. On the
other hand, for the $\pm$xFo-SR states the minus sign is picked up for
$S_{1,z}^{\rm(ss)}$ and the plus sign for $S_{2,z}^{\rm(ss)}$, i.e.,
$S_{l,z}^{\rm(ss)}=(-1)^{l}\sqrt{(N_{l}/2)^{2}-[S_{l,x}^{\rm(ss)}]^{2}}$.
Indeed, the linear-stability analyses confirm that with these choices, all
$\pm$xFi-SR and $\pm$xFo-SR are stable in their entire corresponding parameter
regimes, while the other choices lead to instabilities in some parameter
regimes. Furthermore, Figs. 2(b) and (d) also show that with these choices,
the $z$-component $S_{z}^{\rm(ss)}$ of the total spin and the energy $E_{0}$
of the system change continuously between the $-$zFo-N and $\pm$xFi-SR states,
and the $-$zFi-N and $\pm$xFo-SR states across the phase transitions.

Figure S2: The steady-state behavior of the $x$ and $z$ components of the
total spin, $S_{x}$ and $S_{z}$, and the $x$ component of the staggered spin,
$\Delta S_{x}$, in the parameter plan of $N_{2}/N_{1}$ vs. $\lambda/\kappa$,
similar to the phase diagram of Fig. 1 in the main text. Each column shows the
behavior of a fixed point which is indicate in the top. The white areas in the
figures for the fixed points $\pm$xFi-SR and $\pm$xFo-SR indicate parameter
regimes where the corresponding fixed points do not exist. The parameters are
the same as Fig. 1 in the main text.

In Fig. S2, we show the steady-state behavior of the $x$ and $z$ components of
the total spin, $S_{x}^{(\rm{ss})}$ and $S_{z}^{(\rm{ss})}$, and the $x$
component of the staggered spin, $\Delta S_{x}^{(\rm{ss})}$, in the parameter
plan of $N_{2}/N_{1}$ vs. $\lambda/\kappa$ in accordance with the phase
diagram of Fig. 1 in the main text. Note that +xFi-SR and -xFi-SR, and +xFo-SR
and -xFo-SR are the parity-symmetric pairs, related to one another by the
parity transformation $a\to-a$ and $S_{l,x}\to-S_{l,x}$ (and equivalently
$S_{x}\to-S_{x}$ and $\Delta S_{x}\to-\Delta S_{x}$).

###  I.3 Linear-stability analysis

We begin by writing the mean-field equations of motion (Unconventional Dicke
model: Multistabilities and nonequilibrium dynamics) in the form,

$\displaystyle\frac{\partial\mathbf{X}}{\partial t}=\mathbf{F}(\mathbf{X}),$
(S10)

where $\mathbf{X}=(a,a^{*},\mathbf{S}_{1},\mathbf{S}_{2})^{\top}$ is a vector
of the mean-field variables and $\mathbf{F}(\mathbf{X})$ an $8\times 8$ matrix
obtained from the right-hand side of Eq. (Unconventional Dicke model:
Multistabilities and nonequilibrium dynamics). The fixed points
$\dot{\mathbf{X}}_{\rm ss}=0$ of the system are obtained by solving the
coupled algebraic equations $\mathbf{F}(\mathbf{X}_{\rm ss})=0$, as discussed
in the previous section. Taking a Taylor expansion (up to the linear term) of
the right-hand side of Eq. (S10) around the fixed point $\mathbf{X}_{\rm ss}$
and noting that $\partial_{t}{\mathbf{X}}_{\rm ss}=\mathbf{F}(\mathbf{X}_{\rm
ss})=0$ yields

$\displaystyle\frac{\partial}{\partial
t}\delta\mathbf{X}=\mathbf{J}\delta\mathbf{X},$ (S11)

where $\delta\mathbf{X}=\mathbf{X}-\mathbf{X}_{\rm ss}=(\delta a,\delta
a^{*},\delta\mathbf{S}_{1},\delta\mathbf{S}_{2})^{\top}$ is a vector of
fluctuations, and $\mathbf{J}$ the Jacobian matrix

$\displaystyle\mathbf{J}=\frac{\partial\mathbf{F}(\mathbf{X})}{\partial\mathbf{X}}\Bigg{|}_{\mathbf{X}_{\rm
ss}}=\begin{pmatrix}-(i\omega_{c}+\kappa)&0&-i\lambda&0&0&i\lambda&0&0\\\
0&-(-i\omega_{c}+\kappa)&i\lambda&0&0&-i\lambda&0&0\\\
0&0&0&-\omega_{a}&0&0&0&0\\\ -\lambda S_{1,z}^{\rm(ss)}&-\lambda
S_{1,z}^{\rm(ss)}&\omega_{a}&0&-\lambda(a_{\rm ss}+a_{\rm ss}^{*})&0&0&0\\\
\lambda S_{1,y}^{\rm(ss)}&\lambda S_{1,y}^{\rm(ss)}&0&\lambda(a_{\rm
ss}+a_{\rm ss}^{*})&0&0&0&0\\\ 0&0&0&0&0&0&-\omega_{a}&0\\\ \lambda
S_{2,z}^{\rm(ss)}&\lambda S_{2,z}^{\rm(ss)}&0&0&0&\omega_{a}&0&\lambda(a_{\rm
ss}+a_{\rm ss}^{*})\\\ -\lambda S_{2,y}^{\rm(ss)}&-\lambda
S_{2,y}^{\rm(ss)}&0&0&0&0&-\lambda(a_{\rm ss}+a_{\rm ss}^{*})&0\end{pmatrix}.$
(S12)

Recall that $S_{l,y}^{\rm(ss)}=0$ in all steady states.

###  I.4 Semi-classical dynamics

Figure S3: The nonequilibrium dynamics of the system. The system is prepared
in the initial state (a,b) $-$zFo-N, i.e., $\mathbf{S}_{1}=(0,0,-1)N_{1}/2$
and $\mathbf{S}_{2}=(0,0,-1)N_{2}/2$, and (c,d) +zFi-N, i.e.,
$\mathbf{S}_{1}=(0,0,1)N_{1}/2$ and $\mathbf{S}_{2}=(0,0,-1)N_{2}/2$, with a
small fluctuation seed for the cavity field $a=\delta a$. The ensuing dynamics
of the system are displayed in the phase space of the cavity-field amplitude
$a$ (a,c) and the $\mathbf{S}_{1}$ and $\mathbf{S}_{2}$ spin Bloch spheres
(b,d). The insets in panels (a,c) exhibit the dynamics in the the Bloch
spheres of the total $\mathbf{S}$ and staggered $\Delta\mathbf{S}$ spins. The
color gradient from light blue to darker blue in trajectories indicates
schematically the arrow of the time evolution. The dynamic evolutions lead to
the stable fixed point -xFi-SR in the first case (a,b), and to the +xFo-SR in
the second case (c,d). However, in the first case trajectories are slowly
inward-moving spirals and the time scale to reach the fixed point is much
longer than the second case. Therefore, $\pm$xFi-SR can be viewed as stable
focuses, while $\pm$xFo-SR as stable nodes in this parameter regime. The radii
of the Bloch spheres are indicated in each case, where $\mathbf{S}$ and
$\Delta\mathbf{S}$ have the same radius $(N_{1}+N_{2})/2$. Here, the
parameters are the same as Fig. 3 in the main text.

Here in Fig. S3 we present the semi-classical dynamics of the system prepared
in the initial state (a,b) $-$zFo-N, i.e., $\mathbf{S}_{1}=(0,0,-1)N_{1}/2$
and $\mathbf{S}_{2}=(0,0,-1)N_{2}/2$, and (c,d) +zFi-N, i.e.,
$\mathbf{S}_{1}=(0,0,1)N_{1}/2$ and $\mathbf{S}_{2}=(0,0,-1)N_{2}/2$, with a
small fluctuation seed for the cavity field $a=\delta a$. The system is
attracted to the -xFi-SR fixed point in the former case, and to +xFo-SR in the
latter case. These indicate that the fixed points $-$zFo-N and +zFi-N are
unstable in this parameter regime (recall, however, that +zFi-N is always
unstable) and the attractor of the long-time dynamics depends crucially on the
initial state due to the multistability and the coexistent fixed points.