††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Noise-Resilient Phase Transitions and Limit-Cycles in Coupled Kerr
Oscillators

H. Alaeian  Department of Physics & Astronomy, Purdue University, West
Lafayette, IN 47907, USA  Department of Electrical & Computer Engineering,
Purdue University, West Lafayette, IN 47907, USA  M. Soriente  Institute for
Theoretical Physics, ETH Zurich, 8093 Zürich, Switzerland  Department of
Physics, Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138,
USA  K. Najafi  Department of Physics, Harvard University, 17 Oxford Street,
Cambridge, Massachusetts 02138, USA  IBM Quantum, IBM T.J. Watson Research
Center, Yorktown Heights, NY 10598 USA  S. F. Yelin  Department of Physics,
Harvard University, 17 Oxford Street, Cambridge, Massachusetts 02138, USA

###### Abstract

Driven-dissipative quantum many-body systems have been the subject of many
studies in recent years. They possess unique, novel classes of dissipation-
stabilized quantum many-body phases including the limit cycle. For a long time
it has been speculated if such a behavior, a recurring phenomenon in non-
linear classical and quantum many-body systems, can be classified as a time
crystal. However, the robustness of these periodic dynamics, against quantum
fluctuations is an open question. In this work we seek the answer to this
question in a canonical yet important system, i.e., a multi-mode cavity with
self and cross-Kerr non-linearity, including the fluctuation effects via
higher order correlations. Employing the Keldysh path integral, we investigate
the Green’s function and correlation of the cavity modes in different regions.
Furthermore, we extend our analysis beyond the mean-field by explicitly
including the effect of two-body correlations via the
$2^{\textrm{nd}}$-cumulant expansion. Our results shed light on the emergence
of dissipative phase transitions in open quantum systems and clearly indicate
the robustness of limit-cycle oscillations in the presence of the quantum
fluctuations.

many-body open quantum system, quantum synchronization and limit-cycle,
Keldysh formalism, Bose-Hubbard model, cumulant expansion

##  I Introduction

Understanding the quantum phase transitions (QPT) of many-body systems and
their time evolution are among the main goals of modern physics. Historically,
phase transitions have been introduced and investigated in the thermodynamics
limit of closed systems where they thermalize in the steady-state. In recent
years, however, phase transitions and critical phenomena in driven-dissipative
many-body quantum systems have emerged as a major field of research. Extensive
studies have been exploring the new classes of quantum phase transitions in
the non-equilibrium steady-state (NESS) of open quantum systems, on account of
experimental realization of dissipative quantum simulators Tindall et al.
(2019); Mur-Petit et al. (2020); Soriente et al. (2021).

Experimental platforms, such as cavity arrays, superconducting circuits, and
exciton-polaritons, have put forward versatile testbeds to examine the
interplay between (in)coherent drive, dissipation, and interaction on NESS
phases. This includes the multi-stability and crystallization in driven-
dissipative nonlinear resonator arrays Rodriguez et al. (2016); Cao et al.
(2016); Biondi et al. (2017); Collodo et al. (2019) and spins Landa et al.
(2020), and synchronized switching in the arrays of coupled Josephson
junctions Leib and Hartmann (2014), to name a few. Limit-cycle (LC), aka
synchronization, is an intriguing phase of open quantum systems, where the
dynamics traverses a closed trajectory in the phase space Pikovsky et al.
(2001); Strogatz (2001); Eneriz et al. (2019).

So far, mean-field (MF) LC-phases have been predicted in a wide variety of
open quantum systems including optomechanical resonators Ludwig et al. (2008);
Rodrigues and Armour (2010); Qian et al. (2012); Nation (2013); Lörch et al.
(2014), Rydberg lattices Lee et al. (2011), Bose–Hubbard arrays with cross-
Kerr interactions Jin et al. (2013, 2014); Heinrich et al. (2010); Buca et al.
(2019); Alaeian et al. (2021), Heisenberg lattices Owen et al. (2017), spin
arrays Chan et al. (2015), Dicke model Bhaseen et al. (2012); Gambetta et al.
(2019), and non-linear photonic crystal cavities Takemura et al. (2020). Due
to its underlying assumptions however, MF treatment does not include the
quantum fluctuation effects, and only a few recent studies have investigated
the notion of LC in the context of super-operator spectrum Tindall et al.
(2020).

Figure 1: (a) The schematics of a multi-mode spherical Fabry-Perot cavity
containing a non-linear medium that leads to an effective photon-photon
interaction at the rate $U_{0}$. (b) The bare cavity modes set by the cavity
geometrical features, e.g. the mirrors radius of curvature and the cavity
length.

In this article, we ask the question of what are the elements that appear
typically in the phase diagram of a dissipative non-equilibrium system. The
focus is on directly experimentally realizable systems that allow in
complexity at least for limit cycles. The theoretical methodology will be
built upon typical tools, such as mean-field theory, and expanded as far as
necessary to show and confirm the complex results. The story we tell in this
article, thus, introduces first the physical system, which consists of a
dissipative multi-mode cavity shown in Fig. 1 which is driven on one mode and
connected to two side modes via a Kerr nonlinearity. Depending on the
parameters, a complex phase diagram with various steady states, bistability,
or limit cycles ensues. As a benchmark, this system is compared to the same
system without the dissipation and a much simpler phase diagram that consists
only of two types of steady-state.

This is a complex interacting system where mean-field theory is not
automatically expected to provide the correct results. We therefore support
and correct the results using two different methods, (i) employing a
thermodynamic-limit Keldysh formalism and (ii) including higher-order
correlations. Both these methods confirm the phase diagram and provide the
quantitative corrections. The results and the methods presented in this work
enable one to study the quantum properties of NESS beyond the semi-classical
approximation and shed light on the emergence of dissipative phase transitions
in open quantum systems.

The paper is structured as follows. In section II we introduce our model of a
driven three-mode cavity with self and mutual Hubbard-type interaction subject
to single-photon (coherent) drive. When the cavity-environment can be
described by a Markovian process, the joint density matrix of the cavity
evolves according to a Lindblad master equation. In section III we describe
the physics of the closed system which we later use as a benchmark to
highlight the profound effects of dissipation. Section IV presents the
numerical results of an attractive interaction and details the MF phase
diagram, explains the uniform and multi-stable regions for the pumped mode,
and the emergence of the LC-phase for un-pumped side modes. Moreover, we
elaborate on the spectral features of various phases obtained from Keldysh
Green’s functions, and showcase dissipation-stabilized cross-overs in the open
system. The beyond-MF results employing the $2^{\textrm{nd}}$-cumulant
approach is presented in sub-section IV.3 where we compare the Gaussian
approximation of the cavity states with their corresponding MF results. The
implications of quantum correlation on LC is detailed in this sub-section. In
section V we briefly summarize our Keldysh formalism and the relevant Green’s
function presented in this work. Finally, section VII concludes our results
and presents an outlook for follow-up works.

##  II The model

We consider the dynamics of a three-mode lossy cavity with self- and cross-
Kerr non-linearity, i.e. a Hubbard interaction, subject to a single-photon
coherent drive and a single-photon loss Alaeian et al. (2021). In the laser
rotated frame, the conserved dynamics of such system is given via the
following Hamiltonian ($\hbar=1$)

$\displaystyle~{}\hat{H}=$
$\displaystyle\sum_{m=1}^{3}\Delta_{m}\hat{a}_{m}^{\dagger}\hat{a}_{m}+\frac{U_{0}}{2}\left(\hat{a}_{m}^{\dagger
2}\hat{a}_{m}^{2}\right)$
$\displaystyle\quad+2U_{0}\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}^{\dagger}\hat{a}_{1}\hat{a}_{2}+\hat{a}_{1}^{\dagger}\hat{a}_{3}^{\dagger}\hat{a}_{1}\hat{a}_{3}+\hat{a}_{3}^{\dagger}\hat{a}_{2}^{\dagger}\hat{a}_{3}\hat{a}_{2}\right)$
$\displaystyle\quad+\left(U_{0}\hat{a}_{2}^{\dagger
2}\hat{a}_{1}\hat{a}_{3}+\Omega_{2}~{}\hat{a}_{2}^{\dagger}+\textrm{H.c.}\right)\,,$
(1)

where the $2^{\textrm{nd}}$-mode is subjected to a single-photon coherent
drive at the rate of $\Omega_{2}$ and the frequency of $\omega_{L}$.
$\hat{a}_{m},\hat{a}_{m}^{\dagger}$ is the annihilation and creation operators
of the $m^{\textrm{th}}$-cavity mode, respectively,
$\Delta_{m}=\omega_{m}-\omega_{L}$ is the detuning of the $m^{\textrm{th}}$
cavity mode from the coherent drive, and $U_{0}$ is the interaction rate.

For $U_{0}\geq 0$ the system energy increases with increasing the particle
number hence, a _repulsive_ interaction. Similarly, when $U_{0}\leq 0$ the
system energy decreases with increasing the particle number and the
interaction is _attractive_.

In the presence of an incoherent single-photon loss with the rate of
2$\gamma_{m}$ from the cavity, the Markov-Born approximation leads to the
following Lindblad dissipator for the cavity density operator

$~{}\mathcal{D}(\hat{\rho})=\sum_{m=1}^{3}\gamma_{m}\left(2\hat{a}_{m}\hat{\rho}\hat{a}_{m}^{\dagger}-\\{\hat{a}_{m}^{\dagger}\hat{a}_{m},\hat{\rho}\\}\right).$
(2)

The time evolution of the multi-mode cavity density operator $\hat{\rho}(t)$
is determined via the following master equation

$~{}\frac{d}{dt}\hat{\rho}(t)=-i[\hat{H},\hat{\rho}]+\mathcal{D}(\hat{\rho}).$
(3)

In this work, we are interested in the long-time solution of this density
matrix, where all the transient dynamics are over. This is obtained by
numerically integrating the coupled equations of motions (EoM) [cf. Appendix
A] implemented by the $4^{\textrm{th}}$-order Runge-Kutta method.

When the driving laser is swept from the red to the blue detuning, the
nonlinear dynamics sets in as a parametric amplification process and photons
are created in the side modes as correlated pairs. As shown in Appendix A, for
an attractive interaction, a multi-stability region exists for the red-detuned
driving fields satisfying $\Delta_{2}\geq\sqrt{3}\gamma_{2}$. Similarly, for
repulsive interactions a multi-stability phase exists for the blue-detuned
coherent drives when $\Delta_{2}\leq-\sqrt{3}\gamma_{2}$. For some pumping
rates and detuning, the amplification moves into a regime of self-sustained
oscillations, _aka_ limit cycle, due to the non-linear self- and cross-Kerr
coupling. In the next section, we detail the phase diagram of the
aforementioned system assuming an attractive interaction.

##  III Closed system

To draw a preliminary understanding of our model we first describe the closed
system. We begin by noting that the Hamiltonian of Eq. (II) is unchanged under
the exchange of the un-pumped modes, $a_{1}\leftrightarrow a_{3}$.
Additionally, there are neither terms that mimic a single- nor a two-photon
pump for these modes taken separately. Therefore, a symmetry-breaking due to
population of modes $a_{1}$ and/or $a_{3}$ is only possible if both of them
are populated at the same time. On the contrary, due to the single-photon pump
of the $2^{\textrm{nd}}$-mode, the system features a phase transition as a
function of the pump strength $\Omega_{2}$. We study the mean-field energy
potential landscape, $\bar{H}_{3}$, for a zero occupation of the
$1^{\textrm{st}}$ and $3^{\textrm{rd}}$-mode 111The study of the full
$6-$dimensional energy functional goes beyond the scope of this work..

$\bar{H}_{3}=\Delta_{2}|\alpha_{2}|^{2}+\frac{U_{0}}{2}|\alpha_{2}|^{4}+\Omega_{2}\alpha_{2}^{*}+\Omega_{2}\alpha_{2}\,,$
(4)

where $\alpha_{2}$ is the order parameter. Without loss of generality, we
assume $\Omega_{2}$ to be real and find the potential extrema to be given by

$\displaystyle
U_{0}{\textrm{Re}{(\alpha_{2})}}^{3}+\Delta_{2}\textrm{Re}{(\alpha_{2})}+\Omega_{2}=0\,,$
(5) $\displaystyle\textrm{Im}{(\alpha_{2})}=0\,.$ (6)

This leaves us with a cubic equation in $\textrm{Re}{(\alpha_{2})}$. We study
the discriminant and find the boundaries between the regions with only one and
three possible solutions

$\Omega_{2}=\frac{2}{3\sqrt{3}}\sqrt{\frac{\Delta_{2}^{3}}{|U_{0}|}}\,,$ (7)

where we explicitly assumed $U_{0}<0$ in our case.

In the parameter regime where three solutions are allowed, we find one _low_
(LP) and one _high_ population (HP) phase accompanied by an un-physical one,
i.e. imaginary eigen-frequencies of excitations.

To understand the nature of the extrema of the potential, we find the
excitation spectrum associated with the uniform case Hamiltonian (4) Soriente
et al. (2021). We obtain the eigen-frequencies

$\displaystyle\omega_{3}^{\pm}=\pm\sqrt{\left(n_{2}U_{0}+\Delta_{2}\right)\left(3n_{2}U_{0}+\Delta_{2}\right)}\,,$
(8)

with eigen-vectors

$\mathbf{v}_{\pm}=\begin{pmatrix}\frac{U\alpha_{2}^{2}}{\omega_{3}^{\pm}-\Delta_{2}-2Un_{2}}\\\
1\end{pmatrix}\,,$ (9)

and their associated symplectic norms

$ds_{\pm}^{2}=\left|\frac{U\alpha^{2}}{\omega_{3}^{\pm}-\Delta_{2}-2Un_{2}}\right|^{2}-1\,,$
(10)

where $n_{2}=|\alpha|^{2}$.

The symplectic norm describes the nature of a state of the system. Whenever
all positive(negative) frequency eigenmodes have a positive (negative)
symplectic norm the state of the is the ground state. On the other hand if a
physical state of the system, i.e., a state with real excitation eigenmodes,
has at least one excitation mode with negative (positive) norm for positive
(negative) excitation frequency, then the state is an excited state of the
closed system. Additionally, a positive symplectic norm is associated with
particle-like processes where an external excitation is absorbed by the system
whereas a negative norm underpins a hole-like process where an excitation in
the system is destroyed Soriente et al. (2021). Finally, the symplectic norm
helps identifying the so called negative-mass instabilities Scarlatella et al.
(2019).

In Fig. 2, we plot the three-mode mean-field energy potential of Eq. (4). Due
to the single-photon pump of the $2^{\textrm{nd}}$-mode, the system features a
phase transition as a function of the pump strength $\Omega_{2}$ between an LP
and an HP phase. The phase diagram comprises of two qualitatively-distinct
regions in parameter space with: I, the energy functional has only one clear
extremum, the HP phase; II, three extrema exists including a “proper” LP phase
ground-state, a saddle point, and a HP phase that reveals itself as a maximum,
i.e. an excited-state.

Figure 2: Closed phase diagram of the three-mode uniform case. Two distinct
regions are indicated by their respective mean-field energy potential
($\bar{H}_{3}$) landscape as a function of the real part of the cavity field,
$\textrm{Re}{(\alpha_{2})}$. The dark-blue region indicates the parameter
regime where the LP phase is the ground state of the system and the HP
represents a physically-allowed excited state. The LP phase disappears in the
light-blue region where the HP is the only physical state. Figure 3: (a),(b)
The order parameter, $|\alpha_{2}|^{2}$, and (c)-(d) the excitation spectrum,
$\omega$, on top of the LP (light hues) and HP (dark hues), along the orange
and brown cut lines of Fig. 2, as a function of the pump strength,
$\Omega_{2}$, and detuning $\Delta_{2}$, respectively. Real (solid) and
imaginary (dashed) values with green [blue] hues encode the particle-
[hole-]excitations, i.e. $ds^{2}>0$ [$ds^{2}<0$]. In (a), (c), at the
$\Omega_{c}$ boundary between the light blue and dark-blue regions in Fig. 2,
the LP ceases to exist and the only possible state is the excited-state HP. In
(b), (d), at the $\Delta_{c}$ boundary, the LP appears as a possible state of
the system. The HP always feature a norm swap with negative norm for the
positive eigenmode. The gray lines indicate the third unphysical solution.

In Fig. 3, we plot the order parameter $|\alpha_{2}|^{2}$ and the excitation
spectra on top of the possible solutions along the vertical orange cut (a,c)
and the horizontal brown cut (b,d) of Fig. 2. In both cases the order
parameter $|\alpha_{2}|^{2}$ has a finite value within the entire range but in
region II, two distinct phases with different populations are possible.
Besides, there is a third phase whose spectrum is fully imaginary and
therefore not a physical state (gray lines in Fig. 3(a),(b)). The HP exists
throughout the entire phase diagram but always presents a negative (positive)
symplectic norm with a positive (negative) eigenfrequency. Therefore, we
confirm that it is indeed an excited-state of the closed system. Due to the
presence of the negative interactions a “true ground-state” of the unbounded
Hamiltonian is only possible in a limited region of the parameter space (see
Fig. 2) and is identified with the LP phase, accompanied by a positive
(negative) symplectic norm with a positive (negative) eigenfrequency Soriente
et al. (2020, 2021).

##  IV Open System

For the open system, we determine the phase diagram by initializing the system
in many randomized initial conditions and examining the dynamical stability of
the end points using the Bogoliubov matrix spectrum [cf. Appendix A]. In this
respect, the phase boundaries can be thought of as dynamical phase transitions
separating distinct long-time asymptotic behaviors Heyl et al. (2013).

In addition to multi-stable phases, where the final state depends on the
initial conditions, it is also possible to find regions of parameter space
where no stable fixed point exists. In such cases the system may be attracted
to some time-dependent solutions such as limit cycles 222Note that in the
rotated frame of these oscillation the limit-cycle phase can be considered as
steady-state stagnation points., as found in other coupled nonlinear systems.
In particular, we search for the complete set of stable attractors of the
long-time dynamics, including the fixed points, the multistable coexistence
phases, and time-dependent trajectories. The possible steady-states for the
evolution of the system under Eq. (3) are phases where either the only
populated mode is the pumped mode, $(\alpha_{n=2}\neq 0,\alpha_{n\neq 2}=0)$,
called the _uniform phase_ , or phases where the side modes, $a_{1},a_{3}$,
are also populated and host a _limit-cycle_.

In Fig. 4 we plot the three-mode cavity open phase diagram. It comprises four
different regions: (I) the white region that presents only one uniform stable
steady-state which we identify as the LP phase, (II) the dark blue region of
tri-stability between two uniform phases namely, the LP and the HP one, and a
non-uniform solution where the un-pumped modes showcase a limit-cycle
behavior, (III) the blue region of bi-stability between the LP and the HP
uniform phases, and finally, (IV) the light blue region that presents a HP
uniform phase. At the boundary between region (I) and (IV), denoted as the
striped region, an exceptional-point appears and the LP and HP are smoothly
connected.

When compared to the phase diagram of the closed system [cf. Fig. 2], it
becomes clear that the dissipation has a marked impact. It lifts the boundary
of the closed phase diagram (red dashed line), stabilizes an excited state of
the closed system, the HP, in region (IV), renders the LP as the only
attractor of the dynamics in region (I), and leads to the emergence of limit-
cycles.

Figure 4:  The multi-mode open cavity phase diagram as a function of the
pumped mode detuning $\Delta_{2}$ and the pumping rate $\Omega_{2}$. Phase
transitions and crossover occur between a uniform low population phase (I,
white) and a high population one (IV, light-blue). The dark-blue (II) and blue
(III) regions indicate the parameter regime, where the low and high population
phases are “co-stable”. Moreover, in the former a limit-cycle appears. The
red-dashed line indicates the closed-system phase diagram boundary [cf. Fig
2]. The white-light-blue striped region, separating region I and IV, delimits
an exceptional-point region and hosts a smooth crossover between the low and
high population phases. Parameters are $\Delta_{1,3}=\Delta_{2}\mp 1$ for the
bare-cavity modes and $U_{0}=-1$ for an attractive interaction, in terms of
the cavity-mode decay.

###  IV.1 Steady-state

In this section, we focus on different points of the open phase diagram in
order to clearly expose the nature of the underlying attractors including
their stability. As discussed in Appendix A, the possible fixed points of the
EoM are determined as $\frac{d}{dt}\braket{\hat{a}_{m}}=0$ in Eq. 18. Among
the possible fixed points in the long-time limit, the semiclassical dynamics
eventually evolves towards some stable fixed points for all initial
conditions, known as the steady-states of the system.

We use the mean-field approximation, i.e. we factorize higher order moments as
$\braket{\hat{a}_{m}^{\dagger}\hat{a}_{n}}\approx\braket{\hat{a}_{m}^{\dagger}}\braket{\hat{a}_{n}}$,
and obtain the equations of motion for the MF order parameters $\alpha_{m}$,
with $m=1,2,3$. When $U_{0}\Delta_{2}\geq 0$ the EoM has one stable solution
only as $(\alpha_{n=2}\neq 0,\alpha_{n\neq 2}=0)$, whereas when
$U_{0}\Delta_{2}\leq 0$, there might exist several stable steady-state
solutions [cf. Appendix A]. We self-consistently determine the steady-state of
the system which, in general, may allow limit-cycle solutions. For these
solutions, the long-time limit of the steady-state has the general form of
$e^{it\omega_{\textrm{LC}}}$. Physically, the LC solutions and their
associated frequencies ($\omega_{\textrm{LC}}$) can be thought as the
frequency of the parametrically generated pair, i.e., the un-pumped modes, via
the parametric Kerr process.

To draw a more comprehensive picture of the steady-state behavior of the
system, in Figure 5(a)-(c) we show the mean-field occupation of each cavity
mode, $n_{m}=|\alpha_{m}|^{2}$, versus the pumping rate $\Omega_{2}$ at
various detunings, $\Delta_{2}$ = -5, 0, +5 (vertical dashed lines in Fig. 4),
for an attractive interaction of $U_{0}=-1$, and bare cavity spacing as
$\Delta_{1,3}=\Delta_{2}\mp 1$. The red lines show the behavior of the
$2^{\textrm{nd}}$-mode while the blue lines correspond to the
$1^{\textrm{st}}$\- and the $3^{\textrm{rd}}$-mode populations. The solid and
dashed lines signify the stable and unstable solutions, respectively. For
$U_{0}\Delta_{2}\geq 0$, panels (a)-(b), there is only one stable steady-
state, i.e. the uniform phase HP, but when $U_{0}\Delta_{2}<\sqrt{3}$ we
observe the regions of multi-stability and transitions to the non-uniform
phases, Fig. 5(c).

Figure 5:  Cavity mode occupation vs. coherent pump rate ($\Omega_{2}$) for
$\Delta_{2}=-5,0,+5$ in (a)-(c), respectively [cf. orange dashed cut-lines in
Fig. 4]. In each panel, the red lines show the pumped mode population,
$n_{2}$, and the blue lines show the un-pumped modes occupation,
$n_{1}=n_{3}$. Solid lines show the stable solutions while the dashed lines
correspond to unstable branches. In (a) and (b), there is only a uniform
stable steady-state (HP) whereas in (c) multiple stable steady-states are
possible (LP, LC, HP). Panels (d)-(f) show the time evolution of the mean-
field equations of motion for the corresponding detuning of (a)-(c) for 100
randomized initial conditions. In each panel the upper row shows the
population of the pumped mode ($n_{2}=|\alpha_{2}|^{2}$) and the lower row
shows $\textrm{Real}(\alpha_{1})$. In (d) and (e) the system always evolves
towards a uniform steady-state with only the pumped mode populated (HP). In
(f) there are three possible steady-states, two uniform phases LP, HP (red and
brown) and a limit-cycle (LC, orange). In the limit-cycle case the side-modes
show an oscillatory behavior at the frequency $\omega_{\textrm{LC}}$, as
depicted with light blue traces in the bottom panel (f). All modes have the
same decay rate $\gamma_{0}$ to which other rate parameters are normalized.
The time is in units of $1/\gamma_{0}$. $U_{0}=-1,\Omega_{2}=3$ in all
temporal calculations. The side modes are equally-spaced around the pumped
mode as $\omega_{1,3}=\omega_{2}\mp 1$.

In Fig. 5(d)-(f) we show three exemplary MF time traces for randomized initial
conditions corresponding to the points $T_{1,2,3}$, i.e. $\Omega_{2}=3$, in
Fig. 4. They are obtained from the direct integration of the EoM and
illustrate the implications of a limit-cycle phase and its difference with the
uniform one. For each detuning, the upper row shows the temporal behavior of
$n_{2}=|\alpha_{2}|^{2}$ while the lower panel shows the real part of the
$1^{\textrm{st}}$-mode order parameter, i.e. $\textrm{Re}(\alpha_{1})$.

For $T_{1}$ and $T_{2}$, corresponding to $\Delta_{2}=-5,0$ respectively,
there is only a uniform phase where $\alpha_{1,3}=0$ and $\alpha_{2}\neq 0$
[cf. Appendix A]. Accordingly, in Fig. 5(d),(e) we see that all time traces
converge to only one non-zero value (HP) for the pumped mode and a zero value
for the side modes, independent of the initial conditions.

On the other hand, for $T_{3}$ corresponding to $\Delta_{2}=+5$, where a
multi-stability is predicted [cf. Appendix A], one can see that the time
traces in the upper panel of Fig. 5(f), converge to three different values for
$n_{2}$ at LP (red traces), LC (orange traces), and HP (brown traces). From
those three phases only the LC corresponds to a non-zero values for the
$1^{\textrm{st}},3^{\textrm{rd}}$ modes (orange line in upper panel of the
Fig. 5(f)), where their order parameter $\textrm{Re}(\alpha_{1})$ shows a
periodic long-time behavior (ligh blue in the lower panel of Fig. 5(f)).

###  IV.2 Uniform phase crossover

To illustrate the dissipation effect in modifying the closed-system phases, in
Fig. 6(a,c) and (b,d), we plot the population of the $2^{\textrm{nd}}$ mode,
$|\alpha_{2}|^{2}$, (top) and the fluctuation eigenvalues (bottom) along the
vertical ($\Delta_{2}=+5$) and horizontal cuts ($\Omega_{2}=0.5$) of Fig. 4,
respectively. In both cases, $|\alpha_{2}|^{2}$ behavior highlights a smooth
cross-over from the HP to the LP between region IV and I. As can be inferred
from the phase diagram, this transition traverses an exceptional point region
where the imaginary part of the eigenvalues coalesce to zero and their real
parts split, hosting over- and under-damped fluctuations.

Figure 6: (a,b) The order parameter, $|\alpha_{2}|^{2}$, and (c,d) real
(dashed) and imaginary (solid) part of the excitation spectrum, $\omega$, on
top of the LP (light green) and HP (dark green), along the vertical
($\Delta_{2}=+5$) and horizontal ($\Omega_{2}=0.5$) cut lines in Fig. 4

Using a Keldysh action approach we readily get access to the Green’s function
of the system and its associated dynamical observable. In Fig. 7, we compare
the spectral functions, $\mathcal{A}(\omega)=-2\text{Im}[G^{R}(\omega)]$, of
the three-modes in the LP (a) and HP (b) phases for different detuning and at
fixed pump strengths $\Omega_{2}=0.5$, points $P_{\textrm{LP,HP}}$ in Fig. 4.
From the left panel we see that the LP phase has the response of a “ground-
state”, i.e., a positive (negative) peak at the positive (negative)
frequencies [cf. Sec. III]. We note that peaks at negative frequencies are
unresolved due to the scale resolution. On the other hand, the HP phase
features a peak swap with a positive (negative) peak at the negative
(positive) frequencies, a hallmark of a stabilized excited-state [cf. Sec. III
and Fig. 10 for more information]. We note that the two uniform phases have
different occupations for the $2^{\textrm{nd}}$-mode, see Fig. 5 LP and HP.
Nevertheless, their spectral functions highlight a more profound difference
rather than just a different mode occupation. Interestingly, the peak swap is
also present in the response of the empty side modes and we can think of it as
being in presence of a normal phase to excited normal phase transition
Soriente et al. (2020).

Even though the LP and HP phases are not described by the same order
parameters, through the study of the dynamical responses we identified a peak
swap between their spectral functions and we traced it back to the LP being
the ground-state of the closed system and the HP an excited-state Soriente et
al. (2021); Ferri et al. (2021). The dissipation stabilizes an excited-state,
the HP, and lifts the boundaries of the closed system phase diagram changing
its topography, see Sec. III. The spectral functions shown in Fig. 7 belong to
two points in the region I and IV of the open phase diagram, i.e. where only
one stable attractor exists. The spectral functions at these two points
present a peak swap, signaling a particle- to hole-like physics transition,
when going from the region I to region IV of the open phase diagram.

Figure 7: Spectral function of the three-mode harmonic cavity uniform phase,
for the points $P_{\textrm{LP}}=(\Delta_{2}=+5,\Omega_{2}=0.5)$, (a), and
$P_{\textrm{HP}}=(\Delta_{2}=-5,\Omega_{2}=0.5)$, (b), [cf. Fig. 4]. The
spectral function of $P_{\textrm{HP}}$ features a peak swap with respect to
the one of $P_{\textrm{LP}}$. This signals the presence of a transition
between particle- to hole-like physics going from region I to IV of the open
phase diagram.

###  IV.3 Gaussian approximation and beyond the MF results

So far, we have studied the MF-phase diagram of the open system and
investigated the disspation-stabilized phases in contrast to the closed
system. Moreover, using Keldysh approach, we have been able to investigate the
spectral signatures of each phase and delimit the onset of a PT, beyond the
MF.

As can be inferred from the Langevin EoM in Eq 17 of Appendix A, the quartic
interaction leads to an infinite hierarchy of moments. Therefore, any semi-
analytic or numerical calculations require a truncation of this hierarchy. The
mean-field treatment ignores the correlation via the factorization
approximation, hence truncating the cumulant expansions to the
$1^{\textrm{st}}$-order. Recent studies however, show that the inclusion of
quantum correlations lead to marked deviations from the mean-field results,
especially close to the phase transition points Reiter et al. (2020). So far,
there have been several approaches to include the effect of higher-order
correlation including exact diagonalization, diagrammatic expansions,
functional renormalization group, numerical or density matrix renormalization
group analysis, and phase space methods. Schleich (2001); Bulla et al. (2003);
Schollwöck (2005); Sieberer et al. (2013); Aoki et al. (2014); Täuber and
Diehl (2014); Mascarenhas et al. (2015); Jin et al. (2016); Mathey and Diehl
(2020); Arndt and Hassler (2021).

In this section, we extend the EoM to include the moment dynamics up to the
$2^{\textrm{nd}}$-order, while assuming a vanishing $3^{\textrm{rd}}$-order,
to find the Gaussian approximation of the NESS with an emphasis on the
robustness of the LC-phase. Further details for this specific problem and the
comparison between the results of this approach with the MF and the full
density matrix calculations for a single-mode Kerr cavity can be found in the
Appendix C.

Figure 8(a)-(c) shows the population of the side modes
($\braket{\hat{a}_{1,3}^{\dagger}\hat{a}_{1,3}}$ in blue), the pumped mode
($\braket{\hat{a}_{2}^{\dagger}\hat{a}_{2}}$ in red), and the correlation
between the generated pairs in the side modes
($|\braket{\hat{a}_{1}\hat{a}_{3}}|$ in brown), for $\Delta_{2}=-5,0,+5$,
respectively.

Figure 8: (a)-(c) Population of the pumped mode
($\braket{\hat{a}_{2}^{\dagger}\hat{a}_{2}}$ in red), the un-pumped modes
($\braket{\hat{a}_{1}^{\dagger}\hat{a}_{1}}$ in blue), and the correlation
between the side-modes ($|\braket{\hat{a}_{1}\hat{a}_{3}}|$ in brown) as a
function of the pumping rate $\Omega_{2}$ at detuning of $\Delta_{2}=-5,0,+5$,
respectively. (d)-(c) The time trace of the pumped mode population for
corresponding detuning and at fixed pump rate of $\Omega_{2}=3$.

For panels (a),(b) corresponding to the MF uniform phase, i.e. one solution
for the pumped mode and no occupation of the side modes, the Gaussian
approximation results are in good agreement with the MF ones depicted in Fig.
5(a),(b). For the side-modes however, the MF predicts zero population while
the $2^{\textrm{nd}}$-cumulant suggests a finite but low value. This can be
understood in terms of the quantum fluctuations that had been ignored in the
MF, and partially resumed in the Gaussian approximation. Besides, the weak
correlation between the side modes, i.e. $|g_{13}|$, reinforces that
interpretation.

For $\Delta_{2}=+5$ however, the MF (Fig. 5(c)) and the Gaussian approximation
(Fig. 8(c)) show marked differences. While the MF shows a multistable
behavior, a $1^{\textrm{st}}$-order phase transitions, and an LC behavior for
the side modes, the Gaussian approximation results are continuous for all
modes. For the side modes and before the LC-phase, the trend is quite similar
to the other detuning in panels (a) and (b), i.e. a finite but low population.
Unlike the MF case however, the transition to the LC-phase is not a
$1^{\textrm{st}}$-order as suggested by the MF and instead it shows a large
but continuous change of the order parameter,
$\braket{\hat{a}_{1}^{\dagger}\hat{a}_{1}}$. This transition is accompanied
with a maximum in the correlation between the side modes,i.e.
$\braket{\hat{a}_{1}\hat{a}_{3}}$, a quantity that signifies the correlated
photon-pair generation within this phase. It is interesting to note that this
correlation is low before the LC-phase and drops but remains finite after the
LC. Moreover, the MF predicts a uniform phase again after the LC, where there
is no population in the side modes, while the Gaussian approximation
delineates that a finite population after the LC is the correlation effect,
only.

Figure 8(d)-(f) presents the temporal evolution of the pumped mode for the
corresponding detuning in (a)-(c) at a fixed pumping rate of $\Omega_{2}=3$
starting from vacuum, i.e. $\braket{\hat{a}_{m}^{\dagger}\hat{a}_{m}}=0$. For
the uniform cases in (a),(b) the order parameter settles to a fixed value at
the long-time limit after some transient behavior, a value which is very close
to the corresponding MF value.

Figure 9: (a),(b) The spectrum $\mathcal{A}(\omega)$ of the 3-mode cavity at
$P_{\textrm{LP}}$ and $P_{\textrm{HP}}$ in Fig. 4. Using
$2^{\textrm{nd}}$-cumulant approximation (solid lines) reveals a frequency
pulling towards lower frequencies as compared to MF calculation (dashed
lines). The effect is more pronounced in the high density phase where the
corrections to the MF are bigger.

In contrary, the order parameter of the LC phase shows sustainable
oscillations in time, after a short transient period, corroborating with the
LC-phase features. It is interesting to note that, unlike the MF-predicted LC
however, where the time-periodic behavior was only observable in the side-
modes, here the coupled correlations lead to oscillatory behaviors in all
correlations.

In Fig. 9, we compare the spectral functions of the modes obtained from the
Keldysh approach based on the $2^{\textrm{nd}}$-cumulant results, with the
ones determined from the MF-values. There is no qualitative change in the
response and the peak swap between the LP and HP is still present. The
$2^{\textrm{nd}}$-cumulant approximation reveals corrections to the eigen-
frequencies of the system highlighted as a frequency pulling towards lower
values. The effect is more pronounced in the HP, in agreement with the
presence of stronger corrections to the MF results where the Gaussian
approximation predicts a non-vanishing population for the un-pumped modes in
contrary to zero MF-values.

##  V Keldysh Formalism

The exact diagonalization of the Liouvillian super-operator suffers from the
finite size effects, i.e. the truncation of the Fock space, while the MF
approach ignores quantum fluctuations. To include the correlation effects
while approaching the thermodynamic limit (the typical validity range of a MF
treatment), we employ Keldysh formalism Sieberer et al. (2016).

We readily write the rotated Keldysh action, associated with the Hamiltonian
(II) and subjected to the dissipation (2), in the quantum and classical fields
$a_{m_{c,q}}$, for $m=1,2,3$ as

$S_{k}(\vec{\alpha}_{c,q})=\sum_{m=1}^{3}S_{0}^{m}+S_{\textrm{int}}+S_{\textrm{drive}}+S_{\gamma}\,,$
(11)

where $S_{0}^{m},S_{\textrm{int}},S_{\textrm{drive}},S_{\gamma}$ are defined
in Appendix B. Here, we only highlight that the symmetry $a_{1}\leftrightarrow
a_{3}$ is still present in the action (11).

We perform the saddle point approximation via $\partial
S_{k}/\partial\alpha_{c,q}=0$, set the quantum fields to zero
$\alpha_{m_{q}}=0$, and obtain the EoM for the classical fields
$\alpha_{c_{m}}$ [cf. Appendix B]. These equations coincide with the mean-
field ones upon rotating back to physical fields, i.e.,
$\alpha_{m_{c}}=\sqrt{2}\alpha_{m}$.

To go beyond mean-field, we study the fluctuations around the MF stationary-
states as $\alpha_{m}=\alpha_{m_{c}}+\delta\alpha_{m}$. We expand the Keldysh
action in Eq. (11) and retain terms up to second order in fluctuations, i.e.
only the Gaussian parts. Therefore, the Gaussian action can be written in the
normal form as

$S_{G}=\int_{\omega}\mathbf{\delta\Phi^{\dagger}(\omega)}M(\omega)\mathbf{\delta\Phi(\omega)},$
(12)

where we went to Fourier space with the 3-mode, 12 component, nambu-spinor
$\mathbf{\delta\Phi(\omega)}$, and the matrix $M(\omega)$ is given by

$M_{12\times 12}=\frac{1}{2}\begin{pmatrix}[0]_{6\times
6}&[G^{A}(\omega)]_{6\times 6}^{-1}\\\ [G^{R}(\omega)]_{6\times
6}^{-1}&[P^{K}]_{6\times 6}\end{pmatrix}\,,$ (13)

where $[G^{A}(\omega)]^{-1},[G^{R}(\omega)]^{-1},P^{K}(\omega)$ are the
inverse of the advanced, retarded Green’s functions, and the Keldysh
component, respectively [cf. Appendix B.2].

Using Gaussian integration we determine the single-mode action and the
associated Green’s functions. We here report solely the uniform phase case,
$\alpha_{2}\neq 0,\,\alpha_{1,3}=0$ hence, the single-mode Green’s functions
has simple analytical expressions

$\displaystyle G_{R}^{2}$
$\displaystyle=-\frac{4(\omega+\Delta_{2}+i\gamma_{2}+U_{0}n_{2})}{4(\omega+\Delta_{2}+i\gamma_{2})(-\omega+\Delta_{2}-i\gamma_{2})+8U_{0}\Delta_{2}n_{2}+3U_{0}^{2}n_{2}^{2}}$
(14) $\displaystyle G_{R}^{m}$
$\displaystyle=\frac{4(\omega+\Delta_{\bar{m}}+i\gamma_{\bar{m}}+U_{0}n_{2})}{4(\omega+\Delta_{\bar{m}}+U_{0}n_{2}+i\gamma_{\bar{m}})(\omega-\Delta_{m}-U_{0}n_{2}+i\gamma_{m})+U_{0}^{2}n_{2}^{2}}\,,$
(15)

where $n_{2}=|\alpha_{2}|^{2}$ and $(m,{\bar{m}})=(1,3),(3,1)$.

From the poles of the Green’s function we obtain the eigenvalues

$\displaystyle\omega^{2}$
$\displaystyle=\pm\frac{1}{2}\sqrt{4\Delta_{2}^{2}+8\Delta_{2}U_{0}n_{2}+3U_{0}^{2}n_{2}^{2}}+i\gamma_{2}$
(16) $\displaystyle\omega^{m}$
$\displaystyle=\frac{\Delta_{m}-\Delta_{\bar{m}}+i\gamma_{m}+i\gamma_{\bar{m}}}{2}\pm\frac{1}{2}\sqrt{(\Delta_{m}+\Delta_{\bar{m}}+i\gamma_{m}-i\gamma_{\bar{m}})^{2}+4U_{0}n_{2}(\Delta_{m}+\Delta_{\bar{m}}+i\gamma_{m}-i\gamma_{\bar{m}})+3U_{0}^{2}n_{2}^{2}}$

where $\omega^{i}$ is associated with the $i^{\textrm{th}}$ cavity mode and
$i=1,2,3$.

First, we notice that the pumped mode behaves as a single driven Kerr-
oscillator whereas the symmetric modes $a_{1,3}$, even though empty and not
pumped, present a non-trivial Green’s function, and their Green’s functions
are “coupled”. Second, the eigenvalues we obtained from the Green’s functions
coincide, upon the substitutions $\alpha_{m}=\sqrt{2}\alpha_{m}^{MF}$,
$\gamma_{0}=\gamma_{1}=\gamma_{2}=\gamma_{3}$,
$\Delta_{1}+\Delta_{3}=2(\Delta_{2}-\delta_{D})$, with those obtained from the
stability matrix [cf. Appendix A].

##  VI Experimental Realizations

We consider a multi-mode spherical Fabry-Perot cavity, as sketched in Fig. 1.
Depending on the mirrors radii of curvature and the cavity length, the cavity
leads to an effective harmonic trap for the photons and sets the bare cavity
spectrum, i.e. $\omega_{1,2,\cdots}$ Siegman (1986). Due to the finite
reflectively from the mirrors, there will be a cavity decay rate $\gamma_{n}$
corresponding to the $n^{\textrm{th}}$-mode, which can be approximated with a
fixed rate $\gamma_{0}$, over a finite frequency range of our interest. The
self- and cross-mode interactions are induced via a non-linear medium inside
the cavity, ranging from ultracold atomic gas trapped along the cavity axis
Kollár et al. (2015); Vaidya et al. (2018); Guo et al. (2019), to excitons
Carusotto and Ciuti (2013) or highly-excited Rydberg excitons in solid state.
By re-writing the dynamics in terms of the polaritons, i.e. the hybrid
particles of cavity modes and the material degrees of freedom, we can get the
final dynamics of Eq. II. In brief, $\hat{a}_{n},\hat{a}_{n}^{\dagger}$
represent the annihilation and creation operator of a polariton in the
$n^{\textrm{th}}$-cavity mode, and $U_{0}$ is the effective interaction
between polaritons after approximating the two-body interaction, or the third-
order optical non-linearity, with the short-range contact interaction. By
tuning the relative frequency between the _atoms_ the cavity modes, the
interaction can can be attractive or repulsive. A laser with frequency
$\omega_{L}$ coherently drives the $2^{\textrm{nd}}$-cavity mode. Note that
due to the orthogonality of the Gauss-Lauguerre modes of spherical cavities,
excitation of a specific mode is possible. Although here we only focused on a
canonical cavity-QED setup, recent developments of the multi-mode
superconducting cavities with parametrically-induced self- and cross-Kerr non-
linearity allows ones to engineer similar Hamiltonian using the circuit-QED
toolbox, as well Hung et al. (2021).

##  VII Conclusion and Outlook

For a long time it has been speculated if the limit cycle, a recurring
phenomenon in many non-linear classical and quantum many-body systems, can be
classified as a time crystal. Although the semi-classical MF treatment of
several many-body systems suggests the appearance of periodic long-time
behavior, however, the robustness of these periodic dynamics, against quantum
fluctuations has been an open question. In this work, we investigated this
question in a canonical yet important system, i.e. a multi-mode cavity with
self and cross-Kerr non-linearity.

We used the Keldysh formalism and extended the MF results via a higher-order
cumulant expansion. Both of these methods confirmed the robustness of the LC
phase phase and the overall topology of the phase diagram. Additionally, we
revealed corrections to the MF results, e.g. the change of the correlation
spectra, the non-vanishing population of the side modes, and the oscillating
cross-correlations.

The results and the methods presented in this work can be employed in the
study of a versatile group of the many-body systems to understand the effects
of higher order correlations and investigate the robustness of dissipation-
stabilized MF phases, including the limit cycle. It is straightforward to
extend this work to study the dynamical behavior of such systems subject to a
parametric drive with a modulated amplitude, and investigate the system
transition from a stationary dynamics, to a limit cycle, to a chaotic phase as
a function of the modulation frequency and depth.

Another interesting direction is to extend the system size either by including
more modes in a single cavity or considering an array of coupled cavities with
self- and cross-Kerr non-linearity. When subject to a two-photon parametric
drive, instead of a single-photon as studied here, the $Z_{2}$ symmetry of the
modes are preserved and new phase transitions could appear due to the
spontaneous breaking of this symmetry Bartolo et al. (2016); Puri et al.
(2017). It would be interesting to examine if for some parameter range this
symmetry breaking can be accompanied with a limit-cycle phase. Since such a
system, at certain limits, can approximate a spin Ising Hamiltonian, it is
interesting to investigate the notion of the LC phase on the magnetization.
Finally, one can explore the ultimate behavior of cascaded quantum systems
where the limit cycle output of the first sub-system behaves as a parametric
drive for the second one. This could be thought as an alternative approach for
creating time crystals using discrete time-symmetric drive.

## Acknowledgment

MS would like to thank O. Zilberberg for fruitful discussions and acknowledges
financial support from the Swiss National Science Foundation through the Grant
No. PP00P2_163818. HA acknowledges the financial supports from Baden-
Württemberg Stiftung Eliteprogram award and the Purdue University Startup
fund. SFY would like to thank the NSF through the CUA PFC grant PHY-1734011
and the AFOSR via FA9550-19-1-0233, and the NSF via PHY-1912607.

## References

  * Tindall et al. (2019) J. Tindall, B. B. . J. R. Coulthard, and D. Jaksch, Physical Review Letters 123, 030603 (2019).
  * Mur-Petit et al. (2020) J. Mur-Petit, A. Relano, R. A. Molina, and D. Jaksch, Physical Review Research 3, 024 (2020).
  * Soriente et al. (2021) M. Soriente, T. L. Heugel, K. Arimitsu, R. Chitra, and O. Zilberberg, Physical Review Research 3, 023100 (2021).
  * Rodriguez et al. (2016) S. R. K. Rodriguez, A. Amo, I. Sagnes, L. L. Gratiet, E. Galopin, A. Lemaitre, and J. Bloch, Nature Communications 7, 11887 (2016).
  * Cao et al. (2016) B. Cao, K. W. Mahmud, and M. Hafezi, Physical Review A 94, 063805 (2016).
  * Biondi et al. (2017) M. Biondi, G. Blatter, H. E. Türeci, and S. Schmidt, Physical Review A 96, 043809 (2017).
  * Collodo et al. (2019) M. C. Collodo, A. Potocnik, S. Gasparinetti, J.-C. Besse, M. Pechal, M. Sameti, M. J. Hartmann, A. Wallraff, and C. Eichler, Physical Review Letters 122, 183601 (2019).
  * Landa et al. (2020) H. Landa, M. Schiro, and G. Misguich, Physical Review Letters 124, 043601 (2020).
  * Leib and Hartmann (2014) M. Leib and M. J. Hartmann, Physical Review Letters 112, 223603 (2014).
  * Pikovsky et al. (2001) A. Pikovsky, M. Rosenblum, and J. Kurths, _Synchronization: A Universal Concept in Nonlinear Sciences_ (Cambridge University Press, London, England, 2001).
  * Strogatz (2001) S. H. Strogatz, Nature 410, 268 (2001).
  * Eneriz et al. (2019) H. Eneriz, D. Z. Rossatto, F. A. Cardenas-Lopez, E. Solano, and M. Sanz, Scientific Reports 9, 19933 (2019).
  * Ludwig et al. (2008) M. Ludwig, B. Kubala, and F. Marquardt, New Journal of Physics 10, 095013 (2008).
  * Rodrigues and Armour (2010) D. A. Rodrigues and A. D. Armour, Physical Review Letters 104, 053601 (2010).
  * Qian et al. (2012) J. Qian, A. A. Clerk, K. Hammerer, and F. Marquardt, Physical Review Letters 109, 253601 (2012).
  * Nation (2013) P. D. Nation, Physical Review A 88, 053828 (2013).
  * Lörch et al. (2014) N. Lörch, J. Qian, A. Clerk, F. Marquardt, and K. Hammerer, Physical Review X 4, 011015 (2014).
  * Lee et al. (2011) T. E. Lee, H. Häffner, and M. C. Cross, Physical Review A 84, 031402 (2011).
  * Jin et al. (2013) J. Jin, D. Rossini, R. Fazio, M. Leib, and M. J. Hartmann, Physical Review Letters 110, 163605 (2013).
  * Jin et al. (2014) J. Jin, D. Rossini, M. Leib, M. J. Hartmann, and R. Fazio, Physical Review A 90, 023827 (2014).
  * Heinrich et al. (2010) G. Heinrich, M. Ludwig, J. Qian, B. Kubala, and F. Marquardt, Physical Review Letters 107, 043603 (2010).
  * Buca et al. (2019) B. Buca, J. Tindall, and D. Jaksch, Nature Communications 10, 1730 (2019).
  * Alaeian et al. (2021) H. Alaeian, G. Giedke, I. Carusotto, R. Löw, and T. Pfau, Physial Review A 103, 0137712 (2021).
  * Owen et al. (2017) E. T. Owen, J. Jin, D. Rossini, R. Fazio, and M. J. Hartmann, New Journal of Physics 20, 045004 (2017).
  * Chan et al. (2015) C.-K. Chan, T. E. Lee, and S. Gopalakrishnan, Physical Review A 91, 051601(R) (2015).
  * Bhaseen et al. (2012) M. J. Bhaseen, J. Mayoh, B. D. Simons, and J. Keeling, Physical Review A 85, 013817 (2012).
  * Gambetta et al. (2019) F. M. Gambetta, I. Lesanovsky, and W. Li, Physical Review A 100, 022513 (2019).
  * Takemura et al. (2020) N. Takemura, M. Takiguchi, and M. Notomi, Optics Express 28, 27657 (2020).
  * Tindall et al. (2020) J. Tindall, C. S. Munoz, B. Buca, and D. Jaksch, New Journal of Physics 22, 013026 (2020).
  * Scarlatella et al. (2019) O. Scarlatella, A. A. Clerk, and M. Schiro, New Journal of Physics 21, 043040 (2019).
  * Soriente et al. (2020) M. Soriente, R. Chitra, and O. Zilberberg, Phys. Rev. A 101, 023823 (2020).
  * Heyl et al. (2013) M. Heyl, A. Polkovnikov, and S. Kehrein, Phys. Rev. Lett. 110, 135704 (2013).
  * Ferri et al. (2021) F. Ferri, R. Rosa-Medina, F. Finger, N. Dogra, M. Soriente, O. Zilberberg, T. Donner, and T. Esslinger, _Emerging dissipative phases in a superradiant quantum gas with tunable decay_ (2021), eprint 2104.12782.
  * Reiter et al. (2020) F. Reiter, T. L. Nguyen, J. P. Home, and S. F. Yelin, Physical Review Letters 125, 233602 (2020).
  * Schleich (2001) W. Schleich, _Quantum Optics in Phase Space_ (Wiley-VCH, Weinheim, Germany, 2001).
  * Bulla et al. (2003) R. Bulla, N.-H. Tong, and M. Vojta, Physical Review Letters 91, 170601 (2003).
  * Schollwöck (2005) U. Schollwöck, Review of Modern Physics 77, 259 (2005).
  * Sieberer et al. (2013) L. M. Sieberer, S. D. Huber, E. Altman, and S. Diehl, Physical Review Letters 110 (2013).
  * Aoki et al. (2014) H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Review of Modern Physics 86, 779 (2014).
  * Täuber and Diehl (2014) U. C. Täuber and S. Diehl, Physical Review X 4, 021010 (2014).
  * Mascarenhas et al. (2015) E. Mascarenhas, H. Flayac, and V. Savona, Physical Review A 92, 022116 (2015).
  * Jin et al. (2016) J. Jin, A. Biella, O. Viyuela, L. Mazza, J. Keeling, R. Fazio, and D. Rossini, Physical Review X 6, 031011 (2016).
  * Mathey and Diehl (2020) S. Mathey and S. Diehl, Physical Review B 102, 134307 (2020).
  * Arndt and Hassler (2021) L. Arndt and F. Hassler, Physical Review A 103, 023506 (2021).
  * Sieberer et al. (2016) L. M. Sieberer, M. Buchhold, and S. Diehl, Reports on Progress in Physics 79, 096001 (2016).
  * Siegman (1986) A. E. Siegman, _Lasers_ (University Science Books, 1986).
  * Kollár et al. (2015) A. J. Kollár, A. T. Papageorge, K. Baumann, M. A. Armen, and B. L. Lev, New Journal of Physics 17, 043012 (2015).
  * Vaidya et al. (2018) V. D. Vaidya, Y. Guo, R. M. Kroeze, K. E. Ballantine, A. J. Kollár, J. Keeling, and B. L. Lev, Phys. Rev. X 8, 011002 (2018).
  * Guo et al. (2019) Y. Guo, R. M. Kroeze, V. D. Vaidya, J. Keeling, and B. L. Lev, Phys. Rev. Lett. 122, 193601 (2019).
  * Carusotto and Ciuti (2013) I. Carusotto and C. Ciuti, Rev. Mod. Phys. 85, 299 (2013).
  * Hung et al. (2021) J. S. C. Hung, J. H. Busnaina, C. W. S. Chang, A. M. Vadiraj, I. Nsanzineza, E. Solano, H. Alaeian, E. Rico, and C. M. Wilson, _Quantum simulation of the bosonic creutz ladder with a parametric cavity_ (2021), eprint 2101.03926.
  * Bartolo et al. (2016) N. Bartolo, F. Minganti, W. Casteels, and C. Ciuti, Physical Review A 94, 033841 (2016).
  * Puri et al. (2017) S. Puri, C. K. Andersen, A. L. Grimsmo, and A. Blais, Nature Communications 8, 15785 (2017).
  * Kubo (1962) R. Kubo, Journal of the Physical Society of Japan 17, 1100 (1962).

##  Appendix A The mean-field EoM, dynamical stability, and covariance matrix

Directly from the master equation of Eq. 3 we obtain the Heisenberg-Langevin
equations of motion for the field operators as

$\displaystyle~{}\frac{d}{dt}\hat{a}_{1}$
$\displaystyle=-i\left(\Delta_{1}-i\gamma_{1}\right)\hat{a}_{1}-iU_{0}\left[\left(\hat{a}_{1}^{\dagger}\hat{a}_{1}+2\hat{a}_{2}^{\dagger}\hat{a}_{2}+2\hat{a}_{3}^{\dagger}\hat{a}_{3}\right)\hat{a}_{1}+\hat{a}_{2}^{2}\hat{a}_{3}^{\dagger}\right]+\sqrt{2\gamma_{1}}~{}\hat{\xi}_{1}(t)\,,$
(17) $\displaystyle\frac{d}{dt}\hat{a}_{2}$
$\displaystyle=-i\left(\Delta_{2}-i\gamma_{2}\right)\hat{a}_{2}-iU_{0}\left[\left(\hat{a}_{2}^{\dagger}\hat{a}_{2}+2\hat{a}_{1}^{\dagger}\hat{a}_{1}+2\hat{a}_{3}^{\dagger}\hat{a}_{3}\right)\hat{a}_{2}+2\hat{a}_{2}^{\dagger}\hat{a}_{1}\hat{a}_{3}\right]-i\Omega_{2}+\sqrt{2\gamma_{2}}~{}\hat{\xi}_{2}(t)\,,$
$\displaystyle\frac{d}{dt}\hat{a}_{3}$
$\displaystyle=-i\left(\Delta_{3}-i\gamma_{3}\right)\hat{a}_{3}-iU_{0}\left[\left(\hat{a}_{3}^{\dagger}\hat{a}_{3}+2\hat{a}_{2}^{\dagger}\hat{a}_{2}+2\hat{a}_{1}^{\dagger}\hat{a}_{1}\right)\hat{a}_{3}+\hat{a}_{2}^{2}\hat{a}_{1}^{\dagger}\right]+\sqrt{2\gamma_{3}}~{}\hat{\xi}_{3}(t)\,,$

where $\hat{\xi}_{m}(t)$ is a zero-mean noise operator with the correlations
of
$\braket{\hat{\xi}_{m}(t_{1})\hat{\xi}_{n}^{\dagger}(t_{2})}=\delta_{mn}~{}\delta(t_{1}-t_{2})$
and $\braket{\hat{\xi}_{m}^{\dagger}(t_{1})\hat{\xi}_{n}(t_{2})}=0$.

Ignoring the quantum correlations in the aforementioned EoM and employing the
factorization assumption, one can get a set of coupled non-linear equations
for describing the MF of $\braket{\hat{a}_{m}}=\alpha_{m}e^{i\phi_{m}}$ as

$\displaystyle~{}\frac{d}{dt}\braket{\hat{a}_{1}}$
$\displaystyle=-i\left(\Delta_{1}-i\gamma_{1}\right)\braket{\hat{a}_{1}}-iU_{0}\left[\left(|\braket{\hat{a}_{1}}|^{2}+2|\braket{\hat{a}_{2}}|^{2}+2|\braket{\hat{a}_{3}}|^{2}\right)\braket{\hat{a}_{1}}+\braket{\hat{a}_{2}}^{2}\braket{\hat{a}_{3}}^{*}\right]\,,$
(18) $\displaystyle\frac{d}{dt}\braket{\hat{a}_{2}}$
$\displaystyle=-i\left(\Delta_{2}-i\gamma_{2}\right)\braket{\hat{a}_{2}}-iU_{0}\left[\left(|\braket{\hat{a}_{2}}|^{2}+2|\braket{\hat{a}_{1}}|^{2}+2|\braket{\hat{a}_{3}}|^{2}\right)\braket{\hat{a}_{2}}+2\braket{\hat{a}_{2}}^{*}\braket{\hat{a}_{1}}\braket{\hat{a}_{3}}\right]-i\Omega_{2}\,,$
$\displaystyle\frac{d}{dt}\braket{\hat{a}_{3}}$
$\displaystyle=-i\left(\Delta_{3}-i\gamma_{3}\right)\braket{\hat{a}_{3}}-iU_{0}\left[\left(|\braket{\hat{a}_{3}}|^{2}+2|\braket{\hat{a}_{2}}|^{2}+2|\braket{\hat{a}_{1}}|^{2}\right)\braket{\hat{a}_{3}}+\braket{\hat{a}_{2}}^{2}\braket{\hat{a}_{1}}^{*}\right]\,.$

As can be seen, while the coherent drive restricts
($2\phi_{2}-\phi_{1}-\phi_{3}$), hence breaking the U(1)-symmetry of this
field, there is no additional constraint on the phase of the un-pumped modes.
This phase freedom leads to the emergence of the LC-phase. The stagnation
points of the aforementioned EoM is determined as LHS = 0. When dynamics are
contractive at a particular stagnation point, field operators can be
linearized around that MF with a fluctuation vector $[\delta\hat{\Phi}]$ as

$~{}[\delta\hat{\Phi}]=\begin{bmatrix}\delta\hat{a}_{1}&\delta\hat{a}_{2}&\delta\hat{a}_{3}&\delta\hat{a}_{1}^{\dagger}&\delta\hat{a}_{2}^{\dagger}&\delta\hat{a}_{3}^{\dagger}\end{bmatrix}^{T}\,,$
(19)

where the superscript $T$ means the matrix transpose.

To evaluate the stability of these solutions we employ the dynamical stability
analysis. $\mathcal{M}$, i.e. the Bogoliubov matrix of the small excitation,
has the following structure

$~{}\mathcal{M}=\begin{bmatrix}R&S\\\ S^{*}&R^{*}\end{bmatrix}\,,$ (20)

where

$\displaystyle
R=-i\begin{bmatrix}(\Delta_{1}-i\gamma_{1})+2U_{0}\mathcal{N}&2U_{0}\left(\braket{a_{1}}\braket{a_{2}}^{*}+\braket{a_{2}}\braket{a_{3}}^{*}\right)&2U_{0}\braket{a_{1}}\braket{a_{3}}^{*}\\\
2U_{0}\left(\braket{a_{1}}^{*}\braket{a_{2}}+\braket{a_{2}}^{*}\braket{a_{3}}\right)&(\Delta_{2}-i\gamma_{2})+2U_{0}\mathcal{N}&2U_{0}\left(\braket{a_{1}}\braket{a_{2}}^{*}+\braket{a_{2}}\braket{a_{3}}^{*}\right)\\\
2U_{0}\braket{a_{1}}^{*}\braket{a_{3}}&2U_{0}\left(\braket{a_{1}}^{*}\braket{a_{2}}+\braket{a_{2}}^{*}\braket{a_{3}}\right)&\left(\Delta_{3}-i\gamma_{3}\right)+2U_{0}\mathcal{N}\end{bmatrix},$

and

$\displaystyle
S=-iU_{0}\begin{bmatrix}\braket{a_{1}}^{2}&2\braket{a_{1}}\braket{a_{2}}&\left(\braket{a_{2}}^{2}+2\braket{a_{1}}\braket{a_{3}}\right)\\\
2\braket{a_{1}}\braket{a_{2}}&\left(\braket{a_{2}}^{2}+2\braket{a_{1}}\braket{a_{3}}\right)&2\braket{a_{2}}\braket{a_{3}}\\\
\left(\braket{a_{2}}^{2}+2\braket{a_{1}}\braket{a_{3}}\right)&2\braket{a_{2}}\braket{a_{3}}&\braket{a_{3}}^{2}\end{bmatrix}\,,$
(21)

for $\mathcal{N}=\alpha_{1}^{2}+\alpha_{2}^{2}+\alpha_{3}^{2}$ being the mean-
value of the total number of the photons in the cavity.

When MF is dynamically-stable the Bogoliubov matrix $\mathcal{M}$ is negative-
definite. For the stationary-state of the stable solutions we can define the
covariance matrix as
$\Gamma_{a}(\omega)=\braket{[\delta\hat{\Phi}(\omega)][\delta\hat{\Phi}(\omega)]^{\dagger}}$,
with the following entries

$~{}\Gamma_{a_{mn}}(\omega)=\mathcal{F}\\{\lim_{t\longrightarrow\infty}\braket{\delta\hat{\Phi}_{m}(t+\tau)\delta\hat{\Phi}_{n}^{\dagger}(t)}\\}_{\tau}\,.$
(22)

From this equation it is clear that the spectrum of the correlation matrix is
closely related to the noise operator spectrum. The diagonal entries of the
covariance matrix are related to the cavity transmission spectrum (auto-
correlations), and the off-diagonal entries signify the cross-correlations
between the modes hence, are related to the entanglement between the modes.

###  A.1 uniform phase

Within this phase, only the $2^{\textrm{nd}}$-mode has a non-zero MF and
$\alpha_{1,3}=0$, making $R,S$ diagonal and anti-diagonal matrices,
respectively. The eigenvalues of the dynamical stability matrix $\mathcal{M}$
for identical loss rates of the cavity modes $\gamma_{0}$, are

$\displaystyle\lambda_{1,2}$ $\displaystyle=-\gamma_{0}\pm
i\sqrt{3U_{0}^{2}n_{2}^{2}+4U_{0}\Delta_{2}n_{2}+\Delta_{2}^{2}}\,,$ (23)
$\displaystyle\lambda_{3,4,5,6}$ $\displaystyle=-\gamma_{0}\pm
i\left(\frac{\Delta_{1}-\Delta_{3}}{2}\pm\sqrt{3U_{0}^{2}n_{2}^{2}+4U_{0}(\Delta_{2}-\delta_{D})n_{2}+(\Delta_{2}-\delta_{D})^{2}}\right)\,,$

where the dispersion parameter $\delta_{D}$ is defined as
$2\delta_{D}=2\Delta_{2}-(\Delta_{1}+\Delta_{3})=2\omega_{2}-\omega_{1}-\omega_{3}$.

As can be seen, for $U_{0}\Delta_{2}\geq 0$ and
$U_{0}(\Delta_{2}-\delta_{D})\geq 0$, matrix $\mathcal{M}$ is negative
definite in the uniform phase hence, stagnation points are dynamically
_stable_.

The many-body system is within the uniform phase when $\Omega_{2}$ is either
small or large. If the former, the cross-interaction compared to the self-
interaction is so small that the pair generation process cannot start. In the
latter, the number of particles in the pumped mode is so large that the self-
interaction shifts the pumped mode out of resonance by several $\gamma_{0}$
such that the inter-modal scattering ceases. In other words, in the extreme of
a strong pump, the large self-interaction dominates all many-body interactions
hence, pushing the system to the single-body dynamics again.

On the single-mode branch only, one gets

$~{}\frac{dn_{2}}{d\Omega_{2}}=\frac{2\Omega_{2}}{3U_{0}^{2}n_{2}^{2}+4U_{0}\Delta_{2}n_{2}+\left(\Delta_{2}^{2}+\gamma_{2}^{2}\right)}\,.$
(24)

When $U_{0}\Delta_{2}\leq 0$, there might be points that the slope first
diverges and later changes the sign. These turning points in the MF dynamics
signify the existence of multiple MF attractors. More specifically, the
boundaries of the uniform phase, can be determined as

$~{}n_{2}=\frac{-2U_{0}\Delta_{2}\mp\sqrt{U_{0}^{2}\left(\Delta_{2}^{2}-3\gamma_{2}^{2}\right)}}{3U_{0}^{2}}\,.$
(25)

The above equation determines that for $\delta_{D}=0$ the multi-stability
exists when $\Delta_{2}\geq\gamma_{2}\sqrt{3}$.

A similar argument for a dispersive cavity, i.e. $\delta_{D}\neq 0$, shows
that the multi-stability can exist for $U_{0}(\Delta_{2}-\delta_{D})\leq 0$
and $|\Delta_{2}-\delta_{D}|\geq\gamma_{2}\sqrt{3}$.

The covariance matrix within this uniform phase has the following general form

$\displaystyle\Gamma_{a}(\omega)=\begin{bmatrix}\braket{\hat{a}_{1}(\omega)\hat{a}_{1}^{\dagger}(-\omega)}&0&0&0&0&\braket{\hat{a}_{1}(\omega)\hat{a}_{3}(-\omega)}\\\
0&\braket{\hat{a}_{2}(\omega)\hat{a}_{2}^{\dagger}(-\omega)}&0&0&\braket{\hat{a}_{2}(\omega)\hat{a}_{2}(-\omega)}&0\\\
0&0&\braket{\hat{a}_{3}(\omega)\hat{a}_{3}^{\dagger}(-\omega)}&\braket{\hat{a}_{3}(\omega)\hat{a}_{1}(-\omega)}&0&0\\\
0&0&\braket{\hat{a}_{1}^{\dagger}(\omega)\hat{a}_{3}^{\dagger}(-\omega)}&\braket{\hat{a}_{1}^{\dagger}(\omega)\hat{a}_{1}(-\omega)}&0&0\\\
0&\braket{\hat{a}_{2}^{\dagger}(\omega)\hat{a}_{2}^{\dagger}(-\omega)}&0&0&\braket{\hat{a}_{2}^{\dagger}(\omega)\hat{a}_{2}(-\omega)}&0\\\
\braket{\hat{a}_{3}^{\dagger}(\omega)\hat{a}_{1}^{\dagger}(-\omega)}&0&0&0&0&\braket{\hat{a}_{3}^{\dagger}(\omega)\hat{a}_{3}(-\omega)}\end{bmatrix}\,.$

The structure of this matrix implies that the $2^{\textrm{nd}}$-mode has no
correlation with the side modes. That physically is consistent with the
picture of these un-pumped modes not being populated through the parametric
process via the $2^{\textrm{nd}}$-mode.

Consequently, the covariance matrix for the quadratures gets the following
block diagonal form

$\Gamma_{\textrm{sym}}(\omega)=\begin{bmatrix}\Gamma_{2}(\omega)&\textbf{0}_{4\times
4}\\\ \textbf{0}_{2\times 2}&\Gamma_{\pm}(\omega)\end{bmatrix}\,.$ (26)

The above form emphasizes again that the pumped mode has no correlation with
the other two side modes.

###  A.2 Limit-Cycle phase

As discussed before, while the U(1)-symmetry of the un-driven system is broken
by a coherent drive, the un-pumped modes still have some phase freedom
remained since $\Phi_{0}=2\phi_{2}-\phi_{1}-\phi_{3}$, is the only constraint
imposed by the coherent pump.

Within this phase there is no stationary state and the long-time limit of the
$1^{\textrm{st}},3^{\textrm{rd}}$ MF has an oscillatory behavior as $e^{\pm
it\omega_{\textrm{LC}}}$. In other words, due to the aforementioned phase
constraint, if $\hat{a}_{1}$ oscillates as $e^{it\omega_{\textrm{LC}}}$,
$\hat{a}_{3}$ should vary as $e^{-it\omega_{\textrm{LC}}}$. Going back to the
Eq. 17, one can see that this oscillation can be interpreted in terms of a re-
normalized detuning of the parametrically-populated modes as
$\tilde{\Delta}_{1,3}=\Delta_{1,3}\pm\omega_{\textrm{LC}}$.

##  Appendix B Rotated Keldysh action

In this appendix, we give the explicit expressions for the rotated Keldysh
actions used in the Eq. (11) of the main text.

The Keldysh action of the $m^{\textrm{th}}$-mode with the self-interaction
reads as

$\displaystyle S_{0}^{m}=\alpha_{m_{c}}^{*}\left(i\frac{\partial}{\partial
t}-\Delta_{m}\right)\alpha_{m_{q}}+\alpha_{m_{q}}^{*}\left(i\frac{\partial}{\partial
t}-\Delta_{m}\right)\alpha_{m_{c}}-\frac{U_{0}}{2}\left(|\alpha_{m_{c}}|^{2}+|\alpha_{m_{q}}|^{2}\right)\left(\alpha_{m_{c}}\alpha_{m_{q}}^{*}+\alpha_{m_{c}}^{*}\alpha_{m_{q}}\right)\,.$
(27)

The action due to the cross-term interactions between the modes reads as

$\displaystyle S_{\textrm{int}}=$ $\displaystyle-
U_{0}\left[\left(|\alpha_{1_{c}}|^{2}+|\alpha_{1_{q}}|^{2}\right)\left(\alpha_{2_{c}}\alpha_{2_{q}}^{*}+\alpha_{2_{c}}^{*}\alpha_{2_{q}}+\alpha_{3_{c}}\alpha_{3_{q}}^{*}+\alpha_{3_{c}}^{*}\alpha_{3_{q}}\right)\right]$
(28) $\displaystyle-
U_{0}\left[\left(|\alpha_{3_{c}}|^{2}+|\alpha_{3_{q}}|^{2}\right)\left(\alpha_{2_{c}}\alpha_{2_{q}}^{*}+\alpha_{2_{c}}^{*}\alpha_{2_{q}}+\alpha_{1_{c}}\alpha_{1_{q}}^{*}+\alpha_{1_{c}}^{*}\alpha_{1_{q}}\right)\right]$
$\displaystyle-
U_{0}\left[\left(|\alpha_{2_{c}}|^{2}+|\alpha_{2_{q}}|^{2}\right)\left(\alpha_{1_{c}}\alpha_{1_{q}}^{*}+\alpha_{1_{c}}^{*}\alpha_{1_{q}}+\alpha_{3_{c}}\alpha_{3_{q}}^{*}+\alpha_{3_{c}}^{*}\alpha_{3_{q}}\right)\right]$
$\displaystyle-\frac{U_{0}}{2}\left[\left(\alpha_{2_{c}}^{*^{2}}+\alpha_{2_{q}}^{*^{2}}\right)\left(\alpha_{1_{c}}\alpha_{3_{q}}+\alpha_{1_{q}}\alpha_{3_{c}}\right)+2\alpha_{2_{c}}^{*}\alpha_{2_{q}}^{*}\left(\alpha_{1_{c}}\alpha_{3_{c}}+\alpha_{1_{q}}\alpha_{3_{q}}\right)\right]+c.c.\,.$
(29)

The coherent drive action is

$S_{drive}=-\Omega_{2}\sqrt{2}\left(\alpha_{2_{q}}+\alpha_{2_{q}}^{*}\right)\,.$
(30)

And finally the action due to the Lindblad dissipator reads as

$S_{\gamma}=i\sum_{m=1}^{3}\gamma_{m}\left(2|\alpha_{m_{q}}|^{2}+\alpha_{m_{c}}\alpha_{m_{q}}^{*}-\alpha_{m_{c}}^{*}\alpha_{m_{q}}\right).$
(31)

The action is stationary at the saddle point, determined via $\partial
S_{k}/\partial\alpha_{c,q}=0$, which leads to $\alpha_{m_{q}}=0$ and the
following equations of motion for $\alpha_{m_{c}}$

$\displaystyle\frac{d}{dt}\alpha_{1_{c}}$
$\displaystyle=-i\left(\Delta_{1}-i\gamma_{1}\right)\alpha_{1_{c}}-i\frac{U_{0}}{2}\left[\left(|\alpha_{1_{c}}|^{2}+2|\alpha_{2_{c}}|^{2}+2|\alpha_{3_{c}}|^{2}\right)\alpha_{1_{c}}+\alpha_{2_{c}}^{2}\alpha_{3_{c}}^{*}\right]\,,$
(32) $\displaystyle\frac{d}{dt}\alpha_{2_{c}}$
$\displaystyle=-i\left(\Delta_{2}-i\gamma_{2}\right)\alpha_{2_{c}}-i\frac{U_{0}}{2}\left[\left(|\alpha_{2_{c}}|^{2}+2|\alpha_{1_{c}}|^{2}+2|\alpha_{3_{c}}|^{2}\right)\alpha_{2}+2\alpha_{2_{c}}^{*}\alpha_{1_{c}}\alpha_{3_{c}}\right]-i\Omega_{2}\sqrt{2}\,,$
$\displaystyle\frac{d}{dt}\alpha_{3_{c}}$
$\displaystyle=-i\left(\Delta_{3}-i\gamma_{3}\right)\alpha_{3_{c}}-i\frac{U_{0}}{2}\left[\left(|\alpha_{3_{c}}|^{2}+2|\alpha_{2_{c}}|^{2}+2|\alpha_{1_{c}}|^{2}\right)\alpha_{3_{c}}+\alpha_{2_{c}}^{2}\alpha_{1_{c}}^{*}\right]\,.$

Comparing the EoM for $\alpha_{m_{c}}$ of the action saddle points with Eq. 17
for the MFs, one can see that $\alpha_{m_{c}}=\sqrt{2}\braket{\hat{a}_{m}}$,
which further clarifies the meaning of the classical fields in terms of the
mean-fields.

###  B.1 Approximated Gaussian Action

As mentioned in the text, the effects of quantum fluctuations
$\delta\alpha_{m_{c,q}}$, can be included by expending the Keldysh action
around the saddle points. Since the Hamiltonian of Eq. II is quartic, the
expansion has terms linear - quartic in fluctuations, in general. Here, we
report the explicit form of the Keldysh action up to second order in
fluctuations, i.e. only the Gaussian parts

$\displaystyle S_{k}^{(2)}$
$\displaystyle=\sum_{m=1}^{3}\mathbf{\delta\alpha_{m_{c}}^{*}}\left[\left(i\frac{\partial}{\partial
t}-\Delta_{m}\right)-U_{0}\left(|\alpha_{1_{c}}|^{2}+|\alpha_{2_{c}}|^{2}+|\alpha_{3_{c}}|^{2}\right)-i\gamma_{m}\right]\mathbf{\delta\alpha_{m_{q}}}+i2\gamma_{m}\mathbf{\delta\alpha_{m_{q}}^{*}}\mathbf{\delta\alpha_{m_{q}}}$
(33)
$\displaystyle+\sum_{m=1}^{3}\mathbf{\delta\alpha_{m_{q}}^{*}}\left[\left(i\frac{\partial}{\partial
t}-\Delta_{m}\right)-U_{0}\left(|\alpha_{1_{c}}|^{2}+|\alpha_{2_{c}}|^{2}+|\alpha_{3_{c}}|^{2}\right)+i\gamma_{m}\right]\mathbf{\delta\alpha_{m_{c}}}$
$\displaystyle-
U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)\left(\mathbf{\delta\alpha_{1_{c}}^{*}}\mathbf{\delta\alpha_{2_{q}}}+\mathbf{\delta\alpha_{2_{c}}}\mathbf{\delta\alpha_{1_{q}}^{*}}+\mathbf{\delta\alpha_{2_{c}}^{*}}\mathbf{\delta\alpha_{3_{q}}}+\mathbf{\delta\alpha_{3_{c}}}\mathbf{\delta\alpha_{2_{q}}^{*}}\right)-U_{0}~{}\alpha_{1_{c}}\alpha_{2_{c}}\left(\mathbf{\delta\alpha_{1_{c}}^{*}}\mathbf{\delta\alpha_{2_{q}}^{*}}+\mathbf{\delta\alpha_{2_{c}}^{*}}\mathbf{\delta\alpha_{1_{q}}^{*}}\right)+c.c.$
$\displaystyle-
U_{0}~{}\alpha_{2_{c}}\alpha_{3_{c}}\left(\mathbf{\delta\alpha_{2_{c}}^{*}}\mathbf{\delta\alpha_{3_{q}}^{*}}+\mathbf{\delta\alpha_{3_{c}}^{*}}\mathbf{\delta\alpha_{2_{q}}^{*}}\right)-U_{0}~{}\alpha_{1_{c}}\alpha_{3_{c}}^{*}\left(\mathbf{\delta\alpha_{1_{c}}^{*}}\mathbf{\delta\alpha_{3_{q}}}+\mathbf{\delta\alpha_{3_{c}}}\mathbf{\delta\alpha_{1_{q}}^{*}}\right)+c.c.$
$\displaystyle-\frac{U_{0}}{2}\left(\alpha_{1_{c}}^{2}\mathbf{\delta\alpha_{1_{c}}^{*}}\mathbf{\delta\alpha_{1_{q}}^{*}}+\alpha_{1_{c}}^{*^{2}}\mathbf{\delta\alpha_{1_{c}}}\mathbf{\delta\alpha_{1_{q}}}\right)-\frac{U_{0}}{2}\left(\alpha_{3_{c}}^{2}\mathbf{\delta\alpha_{3_{c}}^{*}}\mathbf{\delta\alpha_{3_{q}}^{*}}+\alpha_{3_{c}}^{*^{2}}\mathbf{\delta\alpha_{3_{c}}}\mathbf{\delta\alpha_{3_{q}}}\right)$
$\displaystyle-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)\left(\mathbf{\delta\alpha_{1_{c}}}\mathbf{\delta\alpha_{3_{q}}}+\mathbf{\delta\alpha_{3_{c}}}\mathbf{\delta\alpha_{1_{q}}}+\mathbf{\delta\alpha_{2_{c}}}\mathbf{\delta\alpha_{2_{q}}}\right)+c.c.\,$

###  B.2 $M$ matrix coefficients

In this appendix, we give the explicit form of the entries of the matrix $M$
in Eq. (13)

$[G^{A}(\omega)]^{-1}$ reads as

$\begin{smallmatrix}\left(\omega-\Delta_{1}-U_{0}N-i\gamma_{1}\right)&-\frac{U_{0}}{2}\alpha_{1_{c}}^{2}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}\alpha_{2_{c}}&-U_{0}\alpha_{1_{c}}\alpha_{3_{c}}^{*}&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)\\\
-\frac{U_{0}}{2}\alpha_{1_{c}}^{*^{2}}&\left(-\omega-\Delta_{1}-U_{0}N+i\gamma_{1}\right)&-U_{0}\alpha_{1_{c}}^{*}\alpha_{2_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}^{*}\alpha_{3_{c}}\\\
-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-U_{0}\alpha_{1_{c}}\alpha_{2_{c}}&\left(\omega-\Delta_{2}-U_{0}N-i\gamma_{2}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{2_{c}}\alpha_{3_{c}}\\\
-U_{0}\alpha_{1_{c}}^{*}\alpha_{2_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&\left(-\omega-\Delta_{2}-U_{0}N+i\gamma_{2}\right)&-U_{0}\alpha_{2_{c}}^{*}\alpha_{3_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)\\\
-U_{0}\alpha_{1_{c}}^{*}\alpha_{3_{c}}&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-U_{0}\alpha_{2_{c}}\alpha_{3_{c}}&\left(\omega-\Delta_{3}-U_{0}N-i\gamma_{3}\right)&-\frac{U_{0}}{2}\alpha_{3_{c}}^{2}\\\
-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}\alpha_{3_{c}}^{*}&-U_{0}\alpha_{2_{c}}^{*}\alpha_{3_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-\frac{U_{0}}{2}\alpha_{3_{c}}^{*^{2}}&\left(-\omega-\Delta_{3}-U_{0}N+i\gamma_{3}\right)\end{smallmatrix}$

$[G^{R}(\omega)]^{-1}$ reads as

$\begin{smallmatrix}\left(\omega-\Delta_{1}-U_{0}N+i\gamma_{1}\right)&-\frac{U_{0}}{2}\alpha_{1_{c}}^{2}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}\alpha_{2_{c}}&-U_{0}\alpha_{1_{c}}\alpha_{3_{c}}^{*}&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)\\\
-\frac{U_{0}}{2}\alpha_{1_{c}}^{*^{2}}&\left(-\omega-\Delta_{1}-U_{0}N-i\gamma_{1}\right)&-U_{0}\alpha_{1_{c}}^{*}\alpha_{2_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}^{*}\alpha_{3_{c}}\\\
-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-U_{0}\alpha_{1_{c}}\alpha_{2_{c}}&\left(\omega-\Delta_{2}U_{0}N+i\gamma_{2}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{2_{c}}\alpha_{3_{c}}\\\
-U_{0}\alpha_{1_{c}}^{*}\alpha_{2_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&\left(-\omega-\Delta_{2}-U_{0}N-i\gamma_{2}\right)&-U_{0}\alpha_{2_{c}}^{*}\alpha_{3_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)\\\
-U_{0}\alpha_{1_{c}}^{*}\alpha_{3_{c}}&-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{2}+2\alpha_{1_{c}}\alpha_{3_{c}}\right)&-U_{0}\left(\alpha_{1_{c}}^{*}\alpha_{2_{c}}+\alpha_{2_{c}}^{*}\alpha_{3_{c}}\right)&-U_{0}\alpha_{2_{c}}\alpha_{3_{c}}&\left(\omega-\Delta_{3}-U_{0}N+i\gamma_{3}\right)&-\frac{U_{0}}{2}\alpha_{3_{c}}^{2}\\\
-\frac{U_{0}}{2}\left(\alpha_{2_{c}}^{*^{2}}+2\alpha_{1_{c}}^{*}\alpha_{3_{c}}^{*}\right)&-U_{0}\alpha_{1_{c}}\alpha_{3_{c}}^{*}&-U_{0}\alpha_{2_{c}}^{*}\alpha_{3_{c}}^{*}&-U_{0}\left(\alpha_{1_{c}}\alpha_{2_{c}}^{*}+\alpha_{2_{c}}\alpha_{3_{c}}^{*}\right)&-\frac{U_{0}}{2}\alpha_{3_{c}}^{*^{2}}&\left(-\omega-\Delta_{3}-U_{0}N-i\gamma_{3}\right)\end{smallmatrix},$

where in the above equations
$N=|\alpha_{1_{c}}|^{2}+|\alpha_{2_{c}}|^{2}+|\alpha_{3_{c}}|^{2}$ is the MF
total number of particles in the cavity.

And finally we have $P^{K}$ as

$\begin{pmatrix}i2\gamma_{1}[I]_{2\times 2}&0&0\\\ 0&i2\gamma_{2}[I]_{2\times
2}&0\\\ 0&0&i2\gamma_{3}[I]_{2\times 2}\end{pmatrix},$ (34)

for $[I]_{2\times 2}$ being an identity matrix so rank 2.

###  B.3 Single mode action

We can determine the single-mode actions of each mode by Gaussian integration.
Considering the steady state with empty mode $a_{1,3}$, Figs. 5(a)-(b),
simplifies the analytical expressions for the single mode Green’s function to

$\displaystyle G_{A}^{2}$
$\displaystyle=-\frac{4(\omega+\Delta_{2}-i\gamma_{2}+U_{0}\left|\alpha_{2}\right|^{2})}{4(\omega+\Delta_{2}-i\gamma_{2})(-\omega+\Delta_{2}+i\gamma_{2})+8U_{0}\Delta_{2}\left|\alpha_{2}\right|^{2}+3U_{0}^{2}\left|\alpha_{2}\right|^{4}}$
(35) $\displaystyle G_{A}^{m}$
$\displaystyle=\frac{4(\omega+\Delta_{\bar{m}}-i\gamma_{\bar{m}}+U_{0}\left|\alpha_{2}\right|^{2})}{4(\omega+\Delta_{\bar{m}}+U_{0}\left|\alpha_{2}\right|^{2}-i\gamma_{\bar{m}})(\omega-\Delta_{m}-U_{0}\left|\alpha_{2}\right|^{2}-i\gamma_{m})+U_{0}^{2}\left|\alpha_{2}\right|^{4}}\,,$
(36)

where $(m,{\bar{m}})=(1,3),(3,1)$. From the poles of the Green’s function we
obtain the eigenvalues

$\displaystyle\epsilon^{2}$
$\displaystyle=\pm\frac{1}{2}\sqrt{4\Delta_{2}^{2}+8\Delta_{2}U_{0}\left|\alpha_{2}\right|^{2}+3U_{0}^{2}\left|\alpha_{2}\right|^{4}}+i\gamma_{2}$
(37) $\displaystyle\epsilon^{m}$
$\displaystyle=\frac{1}{2}(\pm\sqrt{(\Delta_{m}+\Delta_{\bar{m}}+i\gamma_{m}-i\gamma_{\bar{m}})^{2}+4U_{0}\left|\alpha_{2}\right|^{2}(\Delta_{m}+\Delta_{\bar{m}}+i\gamma_{m}-i\gamma_{\bar{m}})+3U_{0}^{2}\left|\alpha_{2}\right|^{4}}+\Delta_{m}-\Delta_{\bar{m}}+i\gamma_{m}+i\gamma_{\bar{m}})$
(38)

Interestingly, the pumped mode behaves as a single driven Kerr-oscillator
whereas the symmetric modes $a_{1,3}$, even though empty, present a non-
trivial Green’s function.

###  B.4 Spectral function

Figure 10: Spectral function, $\mathcal{A}(\omega)$, of the three-mode
harmonic cavity uniform phases, for fixed detuning $\Delta_{2}=+5$ and
$U_{0}=-1$. Green, red and blue correspond to the $1^{\textrm{st}}$-,
$2^{\textrm{nd}}$\- and $3^{\textrm{rd}}$-mode, respectively. LP phase (top)
and HP (bottom) at points $S_{1,2,3}$ in Fig. 4. The LP shows no particular
features whereas the HP shows the typical peak inversion characterizing a
dissipation stabilized excited state.

Using a Keldysh action approach we readily get access to the Green’s function
of the system and associated dynamical observable. In Fig. 10, we compare the
spectral functions, $\mathcal{A}(\omega)=-2\text{Im}[G^{R}(\omega)]$, of the
three-modes in the LP (top) and HP (bottom) phases for different pump
strengths at fixed detuning $\Delta_{2}=+5$, points $S_{1,2,3}$ in Fig. 4.
From the top panel we see that the LP phase has the response of a “ground
state”, i.e., a positive (negative) peak at positive (negative) frequencies,
see Sec. III. We note that peaks at negative frequencies are unresolved due to
the scale resolution. On the other hand, the HP phase features a peak swap
with a positive (negative) peak at negative (positive) frequencies, features
of a stabilized excited state [cf. Sec. III]. We note that the two uniform
phases have different occupation of the $2^{\textrm{nd}}$-mode, see Fig. 5 LP
and HP, nevertheless their spectral functions highlight a more profound
difference between them than just a different mode occupation. Interestingly,
the peak swap is also present in the response of the empty side modes and we
can think of it as being in presence of a normal phase to excited normal phase
transition.

##  Appendix C Beyond MF and $2^{\textrm{nd}}$-order Cumulants

Following the approach of employing cumulant expansion for exploring the many-
body dynamics and phase transitions as detailed in Kubo (1962); Reiter et al.
(2020) we perform an extension to the second order, i.e. writing the three-
body correlations appearing in EoM of Eq. 17 as multiplications of the two-
body and single-body moments. In what follows we detail the approach for the
single-mode Kerr cavity and compare the results with ones obtained from the MF
and the full density matrix solution. The extension to the 3-mode case
directly follows the approach here to derive the EoM for 15 independent first-
and second-order correlations.

###  C.1 single-mode Kerr cavity

For the single-mode driven cavity we have

$\displaystyle\frac{d}{dt}\braket{\hat{a}}=$
$\displaystyle-i\left(\Delta-i\gamma\right)\braket{\hat{a}}-iU_{0}\braket{\hat{a}^{\dagger}\hat{a}^{2}}-i\Omega_{2}=-i\left(\Delta-i\gamma+2U_{0}\braket{\hat{a}^{\dagger}\hat{a}}-2U_{0}\braket{\hat{a}^{\dagger}}\braket{\hat{a}}\right)\braket{\hat{a}}-iU_{0}\braket{\hat{a}^{2}}\braket{\hat{a}^{\dagger}}-i\Omega_{2}\,,$
(39)

where we assumed a vanishing $3^{\textrm{rd}}$-order cumulant hence replacing
the three-body correlation of $\braket{\hat{A}\hat{B}\hat{C}}$ in terms of the
two-body correlation and the single-body expectation values as follow

$\braket{\hat{A}\hat{B}\hat{C}}=\braket{\hat{A}\hat{B}}\braket{\hat{C}}+\braket{\hat{A}\hat{C}}\braket{\hat{B}}+\braket{\hat{B}\hat{C}}\braket{\hat{A}}-2\braket{\hat{A}}\braket{\hat{B}}\braket{\hat{C}}$

Note that the above relation is basically the same as Wick’s theorem results
for the Gaussian states.

Therefore, we have

$\braket{\hat{a}^{\dagger}\hat{a}^{2}}=2\left(\braket{\hat{a}^{\dagger}\hat{a}}-\braket{\hat{a}^{\dagger}}\braket{\hat{a}}\right)\braket{\hat{a}}+\braket{\hat{a}^{2}}\braket{\hat{a}^{\dagger}}\,.$
(40)

As can be seen from Eq. 40, the two-body correlations are needed to describe
the dynamics of $\braket{\hat{a}}$.

$\displaystyle\frac{d}{dt}\braket{\hat{a}^{\dagger}\hat{a}}$
$\displaystyle=-2\gamma\braket{\hat{a}^{\dagger}\hat{a}}+i\Omega_{2}\left(\braket{\hat{a}}-\braket{\hat{a}^{\dagger}}\right)\,,$
(41) $\displaystyle\frac{d}{dt}\braket{\hat{a}^{2}}$
$\displaystyle=-i2\Omega_{2}\braket{\hat{a}}-i2\left(\Delta-i\gamma-2.5U_{0}\right)\braket{\hat{a}^{2}}-i2U_{0}\braket{\hat{a}^{3}\hat{a}^{\dagger}}\,.$

The last term of the second equation involves the four-body correlation of
$\braket{\hat{a}^{3}\hat{a}^{\dagger}}$, rendering the sets of EoM to an non-
closed set.

$\displaystyle\braket{\hat{A}\hat{B}\hat{C}\hat{D}}$
$\displaystyle=\braket{\hat{A}\hat{B}}\braket{\hat{C}\hat{D}}+\braket{\hat{A}\hat{C}}\braket{\hat{B}\hat{D}}+\braket{\hat{A}\hat{D}}\braket{\hat{B}\hat{C}}-2\braket{\hat{A}}\braket{\hat{B}}\braket{\hat{C}}\braket{\hat{D}}\,.$
(42)

Finally, the desired 4-body correlation of the single-mode Kerr cavity reads
as

$\braket{\hat{a}^{3}\hat{a}^{\dagger}}=3\braket{\hat{a}\hat{a}^{\dagger}}\braket{\hat{a}^{2}}-2\braket{\hat{a}}^{3}\braket{\hat{a}^{\dagger}}\,.$
(43)

As can be seen, the last equation is in terms of the single-body expectation
values, i.e. $\braket{\hat{a}},\braket{\hat{a}^{\dagger}}$, and two-body
correlations as $\braket{\hat{a}^{2}},\braket{\hat{a}^{\dagger}\hat{a}}$,
only.

Finally, the EoM for the Gaussian approximation of the solution reads as

$\displaystyle\frac{d}{dt}\braket{\hat{a}}$
$\displaystyle=-i\left(\Delta-i\gamma+2U_{0}\braket{\hat{a}^{\dagger}\hat{a}}-2U_{0}\braket{\hat{a}^{\dagger}}\braket{\hat{a}}\right)\braket{\hat{a}}-iU_{0}\braket{\hat{a}^{2}}\braket{\hat{a}^{\dagger}}-i\Omega_{2}\,,$
(44) $\displaystyle\frac{d}{dt}\braket{\hat{a}^{\dagger}\hat{a}}$
$\displaystyle=-2\gamma\braket{\hat{a}^{\dagger}\hat{a}}+i\Omega_{2}\left(\braket{\hat{a}}-\braket{\hat{a}^{\dagger}}\right)\,,$
$\displaystyle\frac{d}{dt}\braket{\hat{a}^{2}}$
$\displaystyle=-i2\Omega_{2}\braket{\hat{a}}-i2\left(\Delta-i\gamma+0.5U_{0}+3U_{0}\braket{\hat{a}^{\dagger}\hat{a}}\right)\braket{\hat{a}^{2}}+i4U_{0}\braket{\hat{a}}^{3}\braket{\hat{a}^{\dagger}}\,.$

Figure 11 compares the results of the MF-theory (blue lines) with the the full
density matrix solution (red line) and the one obtained from the Gaussian
approximation (green line). As can be seen the Gaussian approximation not only
matches with the full density matrix solutions far away from the phase
transition point but also captures the transition behavior from the LP to the
HP. Specifically, the MF bistability region and the associated hysteresis
behavior disappear and the two branches are connected via a rapid but
continuous change. The remained discrepancy between the
$2^{\textrm{nd}}$-cumulant results and the density matrix within the
transition region is related to the non-Gaussian nature of the cavity state
due to the quartic interaction. We note that the multistable behavior does not
appear in the analytic solution and the quantum solution is unique, while the
semiclassical approach gives multiple dynamically stable solutions.

Figure 11:  Comparison between MF (blue line), $2^{\textrm{nd}}$-cumulant
(green line) i.e. the Gaussian approximation, and the full density matrix (red
line) solutions in a single-mode cavity at $\Delta=+5$. For the MF solutions,
the solid and dashed lines correspond to the stable and unstable solutions,
respectively. $U_{0}$ = -1 in all calculations.