# Dissipative phase transitions in $n$-photon driven quantum nonlinear
resonators

Fabrizio Minganti<EMAIL_ADDRESS>Institute of Physics, Ecole
Polytechnique Fédérale de Lausanne (EPFL), CH-1015 Lausanne, Switzerland
Center for Quantum Science and Engineering, Ecole Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Vincenzo Savona
<EMAIL_ADDRESS>Institute of Physics, Ecole Polytechnique Fédérale de
Lausanne (EPFL), CH-1015 Lausanne, Switzerland  Center for Quantum Science and
Engineering, Ecole Polytechnique Fédérale de Lausanne (EPFL), CH-1015
Lausanne, Switzerland  Alberto Biella<EMAIL_ADDRESS>Pitaevskii
BEC Center, CNR-INO and Dipartimento di Fisica, Università di Trento, I-38123
Trento, Italy

(August 28, 2024)

###### Abstract

We investigate and characterize the emergence of finite-component dissipative
phase transitions (DPTs) in nonlinear photon resonators subject to $n$-photon
driving and dissipation. Exploiting a semiclassical approach, we derive
general results on the occurrence of second-order DPTs in this class of
systems. We show that for all odd $n$, no second-order DPT can occur while,
for even $n$, the competition between higher-order nonlinearities determines
the nature of the criticality allowing for second-order DPTs to emerge only
for $n=2$ and $n=4$. As pivotal examples, we study the full quantum dynamics
of three- and four-photon driven-dissipative Kerr resonators, confirming the
prediction of the semiclassical analysis on the nature of the transitions. The
stability of the vacuum and the typical timescales needed to access the
different phases are also discussed. We also show a first-order DPT where
multiple solution emerge around zero, low, and high-photon number. Our results
highlight the crucial role played by strong and weak symmetries in triggering
critical behaviors, providing a Liouvillian framework to study the effect of
high-order nonlinear processes in driven-dissipative systems, that can be
applied to problems in quantum sensing and information processing.

##  1 Introduction

Nonlinear bosonic systems, such as optical cavities, polaritonic systems,
optomechanical resonators, and superconducting circuits, represent an
extremely rich and versatile tool to explore and simulate nonequilibrium
quantum physics [1, 2, 3]. These systems are intrinsically open, meaning that
particle, energy, and correlations can be gained or lost through the coupling
with the environment [4]. Drives are then applied to these systems, bringing
them out of their thermal equilibrium, and compensating for the losses induced
by the environment. As a result, the complex interplay between driving,
dissipation, and Hamiltonian terms gives rise to a nontrivial open-system
dynamics, leading to nonequilibrium stationary states whose properties differ
from those of closed quantum systems at equilibrium [5, 6, 7].

The symmetries of the drive and dissipators play a fundamental role in
determining both the nature of the steady state and the dynamical properties
of a quantum system [8, 9, 10]. For instance, in a nonlinear photonic cavity
the possibility to exploit nonlinear and engineered pumping schemes, in the
presence of moderate single-particle dissipation, opened venues to the
generation and stabilization of nonclassical states [11, 12, 13]. A pivotal
example in this field is the use of two-photon drives to generate, stabilize,
and control photonic Schrödinger cat states [14, 15], that have been proposed
as a fundamental building block of quantum computing devices [16]. Beyond
their interest for quantum information, parametric processes have been at the
center of intense research, leading to the exploration of their properties
both in classical [17, 18] and quantum configurations [11, 19, 20, 21].

The study and characterization of dissipative phase transitions (DPTs), and
their peculiarities, has been the focus of a huge theoretical and experimental
research, especially concerning the connection of DPTs to multimodality and
metastability [22]. In this scenario, two main distinctions have been drawn in
the characterization of DPTs [10]. First-order DPTs are abrupt and
discontinuous changes in the properties of the system’s steady state [23].
These have been associated with hysteresis and critical slowing down [24, 11,
25], which allow to observe the emergence of metastable dynamics and to study
its competition with the other typical timescales of the system. Key to
understand second-order DPTs – where the steady state transitions
continuously, but it is characterized by a divergent response – are symmetries
and their breaking [11, 26, 27]. In particular, DPTs can be associated with
weak and strong symmetry breaking [28, 29]. Second-order DPTs have been shown
to be key for several technological tasks. The cross-fertilization between
quantum information processing and open system criticality led to innovative
ideas to protect quantum information [28, 30], enhance quantum sensing [31,
32, 33, 34], and review laser theory [35, 36, 37].

In the panorama of DPTs, the parametrically-driven (or two-photon) Kerr
resonator has attracted considerable interest [11, 10, 26, 27, 38]. Indeed, it
provided an ideal test model that displays both first- and second-order
criticalities in different regions of the parameter space [11], and represents
one of the few cases for which a steady-state can be analytically found [39,
40]. Thus, DPTs of this model were extensively studied [11, 19, 26, 38], in
connection with the spectral properties of the Liouvillian [10] and more
exotic phenomena, such as exceptional points and parity-time symmetry breaking
[41]. These findings represented the natural extension of the well known
results about first-order transition in the coherently-driven (one-photon)
Kerr resonator [24, 11, 25, 42, 43, 44, 45], pioneered by Drummond and Walls
[46], and showed how the presence of multi-photon driving and dissipation can
drastically modify the physics of nonlinear bosonic resonators [11, 26, 28].
Remarkably, all these results obtained at the single-resonator level provided
a guideline to investigate emergent phenomena in more complex lattice
architectures [27, 47, 48, 49, 50], allowing drawing important parallelisms.

In this work, we advance these ideas, and explore the DPTs of nonlinear
photonic resonators in the presence of parametric $n$-photon drive and losses,
going beyond the aforementioned $n=1,2$ cases. Beside the fundamental
theoretical interest, our research is further motivated by the nontrivial
achievement of higher-order photon pumping scheme exploiting strong
nonlinearities in superconducting circuits [51, 52, 53], making the study of
these phase transition timely. We provide general results showing that second-
order DPTs can only emerge in even-driven Kerr-like resonators, and we point
out the strong technical limitations to their realizations in $n>4$ driving
schemes. We detail how the nature of DPTs is deeply connected to the presence
of strong and weak symmetries and can be analyzed exploiting a general
theoretical framework based on the structure and on the spectral properties of
the Liouvillian superoperator. We also demonstrate that it is possible to
witness nontrivial phenomena such as multistability, emerging from the
competition between drive, dissipation, and nonlinearity in these models.
Thanks to these results, we discuss possible technological limitation that can
affect protocols based on $n>2$ driving schemes.

We provide a detailed numerical analysis of the full quantum model for the
$n=3$ and $n=4$ cases. In the former, we confirm that the system can only
undergo a first-order phase transitions accompanied by the symmetry breaking
of the discrete weak ${Z}_{3}$ symmetry, as the system parameters are properly
varied. For the 4-photon driven resonator, we show that both a first- and
second-order DPT can occur, accompanied by a breaking of the ${Z}_{4}$,
possibly incurring in multistability. We discuss analogies and differences
between the strong and weak symmetric cases. We verify these results also
within the Liouvillian theory, and we derive a precise phase diagram.

The paper is structured as follows. In Sec. 2 we introduce the model, the
master equation governing the driven-dissipative dynamics, the symmetry
properties of the problem, and their consequences on the Liouvillian spectrum.
In Sec. 3 we discuss the emergence of DPTs in the semiclassical limit, while
in Secs. 4 and 5 we study the full quantum dynamics for $n=3,4$ resonators,
respectively. Finally, in Sec. 6 we draw our conclusions and discuss some
future perspectives.

##  2 The model

We consider a bosonic $n$-driven nonlinear resonator, whose Hamiltonian reads

$\hat{H}_{n}=\sum_{m=1}^{m_{\rm
max}}\frac{U_{m}}{m}\left(\hat{a}^{\dagger}\right)^{m}\hat{a}^{m}+G_{n}\left[\hat{a}^{n}+\left(\hat{a}^{\dagger}\right)^{n}\right],$
(1)

where $\hat{a}$ ($\hat{a}^{\dagger}$) is the bosonic annihilation (creation)
operator. The interaction strengths $U_{m}$ sets the scale of $m$-photon
processes. For instance, $U_{1}$ characterizes the energy of one photon in the
resonator (in the frame rotating at the pump frequency) and rescales the term
$\hat{a}^{\dagger}\hat{a}$, $U_{2}$ is a standard Kerr interaction, and so on.
As detailed also in Appendix A, for a $n$-photon process we should consider at
least processes up to the order $m_{\rm max}=\left\lfloor n/2+1\right\rfloor$,
where $\left\lfloor A\right\rfloor$ indicates the integer part of the number
$A$. As we will see in the following, the high-order $U_{m}$s play a
fundamental role in determining the nature of the DPTs, and for this reason
need to be included in a minimal model. $G_{n}$, instead, represents the
$n$-photon drive amplitude.

Given the dissipative nature of the system, and within the Born and Markov
approximations, the system’s dynamics is ruled by a Lindblad master equation
reading (hereafter we set $\hbar=1$)

$\partial_{t}\hat{\rho}(t)=\mathcal{L}[\hat{\rho}(t)]=-i[\hat{H}_{n},\hat{\rho}(t)]+\gamma\mathcal{D}[\hat{a}]+\eta_{n}\mathcal{D}[\hat{a}^{n}],$
(2)

with
$\mathcal{D}[\hat{O}]=\hat{O}\hat{\rho}(t)\hat{O}^{\dagger}-\\{\hat{O}^{\dagger}\hat{O},\hat{\rho}(t)\\}/2$.
The first term in Eq. (2) rules the coherent (unitary) part of the dynamics,
and follows from Eq. (1), upon an appropriate rescaling of $U_{m}$ due to the
dressing of the cavity eigenmodes by the the environment (Lamb-shift-like
terms [4]). The second and the third terms in Eq. (2) account for the
incoherent one- and $n$-photon losses, respectively. While one-photon
dissipation is an unavoidable feature in any photonic resonator, emerging from
the coupling of the cavity modes with the electromagnetic vacuum, $n$-photon
losses naturally emerge as a byproduct of the engineered processes leading to
$n$-photon drive. Notice that other $m$-photon (with $m\neq 1,\,n$)
dissipative processes can be safely neglected, because their emergence is
linked to the presence of engineered $m$-photon exchanges, and here we are
considering only a single drive acting at each time. Although the Hamiltonian
in Eq. (1) is quite general and platform-independent, we provide a brief
discussion on how such terms can emerge in a superconducting circuit
implementation in Appendix A.

###  2.1 Liouvillian spectrum, symmetries, and their breaking

Our analysis will mainly focus on the steady states $\hat{\rho}_{\rm
ss}^{(k)}$, i.e., that density matrices which do not evolve any more under the
action of the Lindblad master equation (2), defined by

$\partial_{t}\hat{\rho}_{\rm ss}^{(k)}=\mathcal{L}\hat{\rho}_{\rm
ss}^{(k)}=0.$ (3)

$k$ is an index which labels these steady states, and in the present analysis
will be solely tied to the presence of a strong symmetry (see below).
Otherwise, the steady state is unique and will be simply called
$\hat{\rho}_{\rm ss}$. In a thermodynamic limit, formally defined as
$L\to\infty$, the steady-state of an open quantum system can display a
nonanalytical behavior as a function of a generic parameter $\zeta$ [23, 10].
A transition is then the nonanalytical change of some operator $\hat{o}$ as
$\zeta$ crosses the critical point $\zeta_{c}$. Formally, we say that there is
a phase transition of order $M$ if [10]

$\lim_{\zeta\to\zeta_{c}}\left|\frac{\partial^{M}}{\partial\zeta^{M}}\lim_{L\to\infty}\mathrm{Tr}\\!\left[\hat{\rho}_{\rm
ss}^{(k)}(\zeta,L)\hat{o}\right]\right|=+\infty.$ (4)

The $n$-photon driven Kerr resonator explicitly displays a ${Z}_{n}$ symmetry
111We will use the notation ${Z}_{n}$ for the symmetry group, $\hat{Z}_{n}$
for the operator associated with such a symmetry, and $\mathcal{Z}_{n}$ for
the corresponding superoperator.. That is, the transformation

$\hat{a}\to\hat{a}\ e^{i2\pi k/n},\quad k=0,\,1,\,\dots,\,n,$ (5)

leaves the master equation (2) unchanged. However, one can define two types of
symmetries in open quantum systems [8, 9]. For the model under consideration,
these are defined according to the way the operator
$\hat{Z}_{n}=e^{i2\pi\hat{a}^{\dagger}\hat{a}/n}$ acts. One speaks of strong
symmetries if $\hat{Z}_{n}$ commutes with both the Hamiltonian and the jump
operators, i.e.:

$[\hat{Z}_{n},\hat{H}]=[\hat{Z}_{n},\hat{a}]=[\hat{Z}_{n},\hat{a}^{n}]=0.$ (6)

In this case, ${Z}_{n}$ imply the existence of a corresponding conserved
quantity $\langle\hat{Z}_{n}\rangle_{t}\equiv{\rm
Tr}[\hat{\rho}(t)\hat{Z}_{n}]={\rm const}$. The system will display $n$
independent steady states, each one characterized by a different value of
$\langle\hat{Z}_{n}\rangle_{\rm
ss}\equiv\lim_{t\to\infty}\langle\hat{Z}_{n}\rangle_{t}$. In our case, such a
condition is fulfilled if and only if $\gamma=0$ (i.e., the photons are never
lost one-by-one). The presence of a strong symmetry implies that there exist
two superoperators $\mathcal{Z}_{n}^{\rm
L}=\hat{Z}_{n}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\hat{\mathds{1}}$
and $\mathcal{Z}_{n}^{\rm
R}=\hat{\mathds{1}}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\hat{Z}_{n}$
222 The
$\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}$
notation for superoperators indicates that, if
$\mathcal{S}=\hat{A}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\hat{C}$,
then $\mathcal{S}\hat{B}=\hat{A}\hat{B}\hat{C}$. Details can be found Details
in Ref. [54]. , such that

$[\mathcal{L},\mathcal{Z}_{n}^{\rm L,\,R}]=0.$ (7)

A weak symmetry, instead, does not respect the conditions in Eq. (6), and as
such the symmetry of the model does not entail a conserved quantity, meaning
that $\langle\hat{Z}_{n}\rangle_{t}$ changes in time [9, 8]. However, the
superoperator
$\mathcal{Z}_{n}=\hat{Z}_{n}\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}}\hat{Z}_{n}$
commutes with the Liouvillian, i.e.

$[\mathcal{L},\mathcal{Z}_{n}]=0.$ (8)

As a consequence of the conditions in Eqs. (7) and (8), strong and weak
symmetries constrain the structure of the Liouvillian $\mathcal{L}$ and of its
spectrum. A compact and convenient way to discuss symmetries and phase
transitions is via the spectral properties of the Liouvillian [10]. Given any
Liouvillian $\mathcal{L}$, we can introduce its eigenvalues $\lambda_{i}$ and
right eigenoperators $\hat{\rho}_{i}$, defined via the relation

$\mathcal{L}\hat{\rho}_{i}=\lambda_{i}\hat{\rho}_{i},$ (9)

where $\mathbb{R}\mathrm{e}\left[\lambda_{i}\right]\leq 0,\forall i$
represents the decay rates induced by the dissipative dynamics [4, 55].

Figure 1: Sketch of the structure of (one of the) steady-state density matrix
(matrices) $\hat{\rho}_{0}^{(0)}$ ($\hat{\rho}_{0}^{(0,0)}$) for a weak (a, b)
and a strong (c, d) $Z_{3}$ (a,c) and $Z_{4}$ (b,d) symmetries. White
indicates that the matrix element is zero. Comparing the weak and strong
symmetric cases, one notices the effect of one-photon dissipation: the
breaking of the strong symmetry by $\gamma$ results in a mixing of the
populations, and thus in more nonzero elements. Nonetheless, being an
incoherent process, no coherences between different symmetry sectors can be
retained, thus resulting in the coarser steady state structure. Hamiltonian
parameters: (a, c) $U_{1}=5\eta_{3}$, $U_{2}=3\eta_{3}$, $G_{3}=9\eta_{3}$;
(b, d) $U_{1}=5\eta_{4}$, $U_{2}=3\eta_{4}$, $U_{2}=\eta_{4}/5$,
$G_{4}=20\eta_{4}$. Dissipation: (a) $\gamma=\eta_{3}$, (b) $\gamma=\eta_{4}$

####  2.1.1 DPTs and weak symmetries

The presence of a weak symmetry allows refining the discussion on the spectral
properties of the system. The eigenvalues $z_{n}^{(k)}$ of $\mathcal{Z}_{n}$
are the $n$ roots of the unity [indeed, $(\mathcal{Z}_{n})^{n}=1$], that is
$z_{n}^{(k)}=e^{2i\pi k/N}$ for $k=0,\,1,\dots n-1$. Since all the eigenstate
of $\mathcal{L}$ must also be eigenstate of $\mathcal{Z}_{n}$, we can
introduce the “quantum number” $k$, such that, for a weak $Z_{n}$ symmetry,

$\mathcal{Z}_{n}\hat{\rho}_{i}^{(k)}=z_{n}^{(k)}\hat{\rho}_{i}^{(k)},\quad\mathcal{L}\hat{\rho}_{i}^{(k)}=\lambda_{i}^{(k)}\hat{\rho}_{i}^{(k)}.$
(10)

We sort the eigenvalues in such a way that
$\lvert\mathbb{R}\mathrm{e}\left[\lambda_{0}\right]^{(k)}\rvert<\lvert\mathbb{R}\mathrm{e}\left[\lambda_{1}\right]^{(k)}\rvert<\ldots<\lvert\mathbb{R}\mathrm{e}\left[\lambda_{n}\right]^{(k)}\rvert$.
The presence of a symmetry thus implies that the Liouvillian does not mix
eigenoperators with different values of $k$, and therefore the Liouvillian can
be partitioned (block-diagonalized) into different symmetry sectors
$\mathcal{L}_{k}$, i.e.,

$\mathcal{L}=\bigoplus_{k}\mathcal{L}_{k}.$ (11)

For this reason, the eigenvalues ${\lambda_{j}}^{(k)}$ and eigenoperators
$\hat{\rho}_{j}^{(k)}$ describe the whole physics within each of the
Liouvillian symmetry sectors.

Weak symmetries fix the structure of the eigenoperators, which, in the number
(Fock) basis, read

$\hat{\rho}_{j}^{(k)}=\sum_{p,q}c_{p,q}\ket{p}\bra{q}\,,\quad\operatorname{mod}(p-q,n)=k,$
(12)

where $\operatorname{mod}(p-q,n)$ indicates the modulo operation. In other
words, $\hat{\rho}_{j}^{(k)}$ must be an operator containing only elements
such that $(m-n)$ is either $k$, or $k\pm n$, or $k\pm 2n$, etc. For example,
for a $Z_{2}$ symmetry, this implies $(m-n)$ either even or odd, and therefore
the eigenoperators of the Liouvillian must be characterised by a checkerboard-
like structure. We show a typical steady-state structure for a weak $Z_{3}$
and $Z_{4}$ symmetries in Figs. 1(a) and (b), respectively. As also
demonstrated in Ref. [9], in the case of a weak symmetry, $\hat{\rho}_{\rm
ss}$ is generally unique and thus must belong to the $k=0$ symmetry sector of
the Liouvillian. For this reason, for any finite number of photons in the
system, the $n$-photon driven Kerr resonator with weak $Z_{n}$ symmetry will
admit a unique steady state $\hat{\rho}_{\rm ss}\propto\hat{\rho}_{0}^{(0)}$.

Furthermore, the discontinuous behavior of the steady state in Eq. (4) is
signalled by the Liouvillian spectral properties. In the thermodynamic limit,
a second eigenoperator, which is stationary under the action of the
Liouvillian, emerges. Accordingly, an eigenvalue $\lambda_{m}^{(k)}$ becomes
exactly zero, both in its real and imaginary parts, as a function of the
parameter $\zeta$. In finite-size systems, phase transitions cannot be
observed, and $\lambda_{m}^{(k)}\neq 0$ if $m\neq 0$ and $k\neq 0$.
Nevertheless, the study of the Liouvillian spectral properties provides much
useful information about the scaling and nature of the transition [47].

Within this formalism and notation, a first order phase transition can be seen
as a change in the $k=0$ symmetry sector, where the steady state
$\hat{\rho}_{\rm ss}\propto\hat{\rho}_{0}^{(0)}$ and the eigenoperator
$\hat{\rho}_{1}^{(0)}$ display level touching (detail can be found in Ref.
[10]). More specifically, $\lambda_{1}^{(0)}=0$ at the critical point, and the
the minimum in $\lambda_{1}^{(0)}$ reaches zero as the system scales towards
the thermodynamic limit.

A spontaneous symmetry breaking of ${Z}_{n}$, instead, means the emergence of
$n-1$ states, each one belonging to a different $k$-symmetry sector, that do
not evolve any more under the action of the Liouvillian. In this case, the
phase transition is associated to $\lambda_{0}^{(1)},\dots\lambda_{0}^{(n-1)}$
becoming and remaining zero in a whole region where the symmetry is broken.
For instance, in the case of a $\mathcal{Z}_{2}$ breaking,
$\lambda_{0}^{(1)}=0$ after the transition, while for $\mathcal{Z}_{3}$ one
has $\lambda_{0}^{(1)}=\lambda_{0}^{(2)}=0$. The corresponding states
$\hat{\rho}_{0}^{(k)}$, belonging to different symmetry sectors with respect
to $\hat{\rho}_{0}^{(0)}$, allow constructing the symmetry-breaking steady
states. Indeed, by choosing the correct superposition of the form
$\hat{\tilde{\rho}}_{j}=\sum_{k}^{n-1}c_{i,k}\hat{\rho}_{0}^{(k)}$, one can
obtain well-defined density matrices such that
$\mathcal{Z}_{n}\hat{\tilde{\rho}}_{j}\neq z_{n}^{(k)}\hat{\tilde{\rho}}_{j}$
but $\mathcal{L}\hat{\tilde{\rho}}_{j}=0$.

####  2.1.2 DPTs and strong symmetries

In the case of a strong symmetry any eigenoperator is characterized by two
quantum numbers $(k_{\rm L},k_{\rm R})$, such that $\mathcal{Z}_{n}^{\rm
L,\,R}\hat{\rho}_{i}^{(k_{\rm L},k_{\rm R})}=e^{\pm 2i\pi k_{\rm
L,\,R}/n}\hat{\rho}_{i}^{(k_{\rm L},k_{\rm R})}$, where, again, $k_{\rm
L,\,R}=0,1,\dots n$. We deduce that

$\begin{split}\hat{\rho}_{i}^{(k_{\rm L},k_{\rm
R})}&=\sum_{p,q}c_{p,q}\ket{p}\bra{q}\,,\\\ \operatorname{mod}(p,n)=k_{\rm
L},&\quad\operatorname{mod}(q,n)=k_{\rm R}.\end{split}$ (13)

Such a structure is shown in Figs. 1(c) and (d) for two steady states of the
$(0,0)$ symmetry sector for $Z_{3}$ and $Z_{4}$ symmetries.

Notice now that we can define two different types of eigenoperators: those
which describe the evolution of populations, for which $k_{\rm L}=k_{\rm R}$,
and the coherences, for which $k_{\rm L}\neq k_{\rm R}$. Consequently, the
symmetry sectors are $\mathcal{L}_{k_{\rm L},k_{\rm R}}$, i.e.,

$\mathcal{L}=\bigoplus_{k_{\rm L},k_{\rm R}}\mathcal{L}_{k_{\rm L},k_{\rm
R}}.$ (14)

For each of the population sectors, there must exist a well-defined steady
state $\hat{\rho}_{\rm ss}^{(k)}\propto\hat{\rho}_{0}^{(k,k)}$(trace one,
Hermitian, and positive semidefinite matrix which does not evolve under the
action of the Liouvillian), while coherences are always traceless matrices.
Accordingly, the definition of the phase transition and of the spontaneous
symmetry breaking takes into account the presence of multiple disconnected
eigenspaces.

A first-order DPT occurs in the population sectors, and it is associated with
the presence of an eigenoperator $\hat{\rho}_{1}^{(k,k)}$ whose eigenvalue
$\lambda_{1}^{(k,k)}$ becomes zero in the thermodynamic limit.

A spontaneous symmetry breaking, instead, implies that the eigenoperators
$\hat{\rho}_{0}^{(k_{\rm L},k_{\rm R})}$ acquire an eigenvalue
$\lambda_{1}^{(k_{\rm L},k_{\rm R})}=0$. Spontaneous symmetry breaking thus
implies that quantum superpositions between the states composing
$\hat{\rho}_{\rm ss}^{(k_{\rm L})}\propto\hat{\rho}_{0}^{(k_{\rm L},k_{\rm
L})}$ and $\hat{\rho}_{\rm ss}^{(k_{\rm R})}\propto\hat{\rho}_{0}^{(k_{\rm
R},k_{\rm R})}$, i.e., two steady states of different symmetry sectors, can be
maintained indefinitely. Indeed, not only the populations do not evolve, but
also the coherences remain stationary. In this regard, DPTs accompanied by
spontaneous breaking of strong symmetries bear a closer resemblance with
Hamiltonian transitions, and this is the reason for their use in quantum
information [28, 30].

##  3 Semiclassical analysis of the $n$-photon driven resonator

The equation of motion for the expectation value of the observable $\hat{a}$
evolving under Eq. (2) is

$\begin{split}\partial_{t}\langle\hat{a}\rangle_{t}&=-i\sum_{m=1}^{m_{\rm
max}}U_{m}\langle\left(\hat{a}^{\dagger}\right)^{m-1}\hat{a}^{m}\rangle_{t}-inG_{n}\langle\left(\hat{a}^{\dagger}\right)^{n-1}\rangle_{t}\\\
&\quad-\frac{\gamma}{2}\langle\hat{a}\rangle-\frac{n\eta_{n}}{2}\langle\left(\hat{a}^{\dagger}\right)^{n-1}\hat{a}^{n}\rangle_{t}.\end{split}$
(15)

Due to the presence of non-quadratic terms, these equations of motion cannot
be closed, leading to a hierarchy of coupled equations.

###  3.1 The thermodynamic limit and finite-component phase transitions

We now introduce the dimensionless parameter $L$ such that

$G_{n}=\tilde{G}_{n}/\sqrt{L^{n-2}},\,U_{m}=\tilde{U}_{m}/L^{m-1},\,\eta_{n}=\tilde{\eta}_{n}/L^{n-1},$
(16)

and we will consider the thermodynamic limit $L\to\infty$. In such a limit
$(G_{n})^{\alpha}U_{m}$ and $(G_{n})^{\beta}\eta_{n}$ are constants [for
$\alpha$ and $\beta$ such that $(-n/2+1)\alpha-m+1=0$ and
$(-n/2+1)\beta-n+1=0$], but the number of excitations diverges. Such a
rescaling of the system parameters can be seen as the generalization of the
scaling proposed in Ref. [56] for the $n=1$ case. The semiclassical (coherent
state) approximation amounts to assuming that the state of the resonator is
coherent, i.e.,

$\hat{\rho}(t)=\ket{\alpha(t)}\bra{\alpha(t)},$ (17)

where $\hat{a}\ket{\alpha(t)}=\alpha(t)\ket{\alpha(t)}$. Accordingly, the
equation of motion for the rescaled coherent field
$\tilde{\alpha}(t)=\braket{\hat{a}}/\sqrt{L}$ leads to a generalized driven-
dissipative Gross-Pitaevskii-like equation

$\begin{split}\partial_{t}\tilde{\alpha}&=\left[-i\sum_{m}\tilde{U}_{m}|\tilde{\alpha}|^{2(m-1)}-\frac{n}{2}\tilde{\eta}_{n}|\tilde{\alpha}|^{2(n-1)}\right.\\\
&\quad\left.-\frac{\gamma}{2}\right]\tilde{\alpha}-in\tilde{G}_{n}\left(\tilde{\alpha}^{*}\right)^{n-1}.\end{split}$
(18)

Equation Eq. (18) is independent of $L$ and the photon number scales as
$N=|\alpha|^{2}\propto L$ confirming that $L\to\infty$ corresponds to a well
defined thermodynamic limit with an infinite number of photons. In general, we
expect the semiclassical approximation (17) to be valid and predictive in the
$L\to\infty$ limit, and far from the critical points where nonlinear processes
inducing quantum fluctuations cannot be neglected.

The parameter $L$ allows introducing the idea of finite-component phase
transitions — where the thermodynamic limit is replaced by a scaling of the
system parameters [24, 57, 56, 11, 58, 59, 36, 60].

###  3.2 Analysis of the transition properties

Given the invariance of Eq. (18) to the transformations in Eq. (16), in the
following analysis we will work with the bare quantities
$\\{U_{m},\eta_{n},G_{n}\\}$.

Despite the simplification introduced by the semiclassical approximation, Eq.
(18) cannot be yet analytically solved. At the steady state, i.e.
$\partial_{t}\alpha=0$, Eq. (18) reads

$\left(\sum_{m=1}^{m_{\rm
max}}U_{m}N^{m-1}-i\frac{\gamma+n\eta_{n}N^{n-1}}{2}\right)\alpha=nG_{n}\left(\alpha^{*}\right)^{n-1}.$
(19)

In general, Eq. (19) gives rise to multiple solutions for the photonic field
$\alpha$. The onset of new stable solutions of Eq. (19) can be associated with
the emergence of phase transitions. For example when $n=1$ the cubic equation
for $\alpha$ (obtained by considering $m_{\rm max}=2$) gives rise to the well-
known S-shaped curve for the photon number [46] signaling the presence of a
first-order phase transition (between a low- and high-density state)
accompanied by an hysteresis region with multiple stable solutions in the
thermodynamic limit [56, 11, 24].

If $n\geq 2$ Eq. (19) always admits the solution $\alpha=\alpha_{\rm vac}=0$.
However, the other solutions $\alpha$ of Eq. (19) cannot be analytically
found. Our strategy is thus to solve an inverse problem. Being interested in
studying the emergence of criticalities as the driving strength is varied, by
multiplying both sides by their complex conjugate, one finally obtains the
equation

$G_{n}(N)=\sqrt{\frac{4\left(\sum_{m}U_{m}N^{m-1}\right)^{2}+{\left(\gamma+n\eta_{n}N^{n-1}\right)^{2}}}{4n^{2}N^{n-2}}},$
(20)

where we selected the positive branch of the square root since, up to a phase,
one can always choose $G_{n}\in\mathbb{R}^{+}$ 333This amounts to a change the
initial condition by sending $\hat{a}\to\hat{a}e^{i\varphi_{0}}$.

####  3.2.1 Stability of the vacuum

In this section we show that, within the semiclassical picture, the solution
$\alpha_{\rm vac}=0$ is always asymptotically stable for $n>2$ if $\gamma\neq
0$.

Consider $\alpha=\alpha_{\rm vac}+\delta\alpha$, where
$\delta\alpha\in\mathbb{C}$ is a small perturbation ($|\delta\alpha|\ll 1$)
around the vacuum solution. Plugging the above parametrization into Eq. (18),
and expanding it at the first order in $\delta\alpha$, we get

$\partial_{t}({\delta\vec{\alpha}})=\mathsf{M}\cdot{\delta\vec{\alpha}},$ (21)

where $\delta\vec{\alpha}=({\rm Re}[\delta\alpha],{\rm
Im}[\delta\alpha])^{\intercal}$ and

$\mathsf{M}=\begin{pmatrix}-\gamma/2&-2\delta_{n,2}G_{2}+U_{1}\\\
-2\delta_{n,2}G_{2}-U_{1}&-\gamma/2\end{pmatrix}$ (22)

is the so-called stability matrix. The solutions of Eq. (21) are given by
$\delta\vec{\alpha}(t)=\exp(-\lambda_{\pm}t)\delta\vec{\alpha}(0)$, where
$\lambda_{\pm}=-\gamma/2\pm\sqrt{4\delta_{n,2}(G_{2})^{2}-U_{1}^{2}}$ are the
eigenvalues of $\mathsf{M}$. Thus it is straightforward to conclude that for
$n>2$ the vacuum solution $\alpha_{\rm vac}=0$ is always stable at a
semiclassical level for finite single-photon losses since

${\rm Re}[\lambda_{\pm}]=-\frac{\gamma}{2}<0.$ (23)

For $n=2$ the vacuum gets unstable when

${\rm Re}\left[\sqrt{4(G_{2})^{2}-U_{1}}\right]>\frac{\gamma}{2},$ (24)

which implies ${\rm Re}[\lambda_{+}]>0$. Equation (24) has important
consequences since it implies that, contrary to the $n=1,2$ case, the
semiclassical dynamics never triggers a transition from the vacuum to a high-
density solutions if $\gamma\neq 0$. However, as we will see in in Secs. 4 and
5, quantum fluctuation in finite-size systems can make the vacuum solution
unstable and allow for the onset of phase transitions.

Finally, we note that in the case of strong symmetry $\gamma=0$, the vacuum is
marginally stable, and higher-order perturbation theory is needed to assess
the stability of the vacuum.

####  3.2.2 Second-order phase transitions and behavior around $\alpha=0$

Figure 2: (a) Possible different behaviors of $G_{n}$ as a function of $N$
according to Eq. (20). The marker indicates the minima of the function
$G_{n}$, while the hatching indicates the unphysical solutions $N<0$. (b) By
simply inverting the plot, we can gain information on $N$ as a function of
$G_{n}$.

For the class of systems under consideration a second-order phase transition
occurs when the state changes from $N=\braket{\hat{a}^{\dagger}\hat{a}}=0$ to
$N>0$ continuously as a function of the driving strength $G_{n}$ [11]. In
other words, if a second-order DPT occurs, semiclassicaly the critical point
must correspond to a solution of Eq. (20) where

$G_{n}^{(c)}\equiv\lim_{N\to 0^{+}}G_{n}(N).$ (25)

We note that the limit $N\to 0^{+}$ must be taken since Eq. (20) is defined
only for $N\neq 0$. At this specific value of $G_{n}$ the system is thus
allowed to pass from the semiclassical solution $\alpha_{\rm vac}=0$ to
another stable solution with $\alpha\neq 0$.

Notice that Eq. (20) admits at most three possible behaviors around $N=0$ as
sketched in Fig. 2:

  * •

The curve $G_{n}$ intersect the zero with a positive derivative (red line in
Fig. 2). In this case, the system can undergo a second-order DPT, passing
continuously from the zero solution to a nonzero one.

  * •

The curve $G_{n}$ intersect the zero with a negative derivative (blue line in
Fig. 2). This is resembling the S-like shape of bistability in one-photon
driven systems. Since the photon number should monotonically increase by
increasing the photon drive, this solution can never be stable, and therefore
the system can only undergo a first-order DPT.

  * •

The curve $G_{n}$ never intersect the zero for a finite value of $N$ (green
line in Fig. 2). Also in this case the system can never experience a second-
order DPT.

We conclude that, a necessary (but not sufficient) condition to observe
second-order DPTs, according to the semiclassical theory, is

${\rm(i)}\quad 0<G_{n}^{(c)}<\infty.$ (26a)
${\rm(ii)}\quad\left.\frac{\partial G_{n}(N)}{\partial
N}\right|_{G_{n}=G_{n}^{(c)}}\geq 0.$ (26b)

Notice that, in the case of a vertical-tangent point, higher-order derivatives
need to be computed.

####  3.2.3 Universal features and a no-go theorem

From the remarks in the previous section and using Eqs. (26a) and (26b), we
can already draw some important conclusions about the nature of DPTs in this
class of systems. In particular, we formulate the following no-go theorem.

No-go theorem. Consider a $n$-photon driven-dissipative resonator, with
nonvanishing Kerr nonlinearity, governed by the Lindbladian (2) then: (a) a
second-order DPT never occurs for odd $n$; (b) If $\gamma\neq 0$ (weak
symmetry), $n=2$ is the only case with a second-order DPT; (c) If $\gamma=0$,
and $U_{1}\neq 0$, again $n=2$ is the only cases where a second-order DPT can
emerge; (d) A DPT for $n=4$ can be found only if $U_{1}=\gamma=0$. (e) For
$U_{2}\neq 0$, no second-order DPTs can occur if $n>4$.

###### Proof.

The semiclassical solutions of the stationary Gross-Pitaevskii equation (19)
must satisfy Eq. (20). The behaviour of this function around $N=0^{+}$ for
$U_{1}\neq 0$ or $\gamma\neq 0$ is given by

$G_{n}(N)\simeq\frac{1}{2n}N^{\frac{2-n}{2}}\sqrt{4U_{1}^{2}+\gamma^{2}}\left[1+\mathcal{O}(N)\right].$
(27)

The case where $U_{1},\gamma=0$, the expansion leads to

$G_{n}(N)\simeq\frac{|U_{2}|}{n}N^{\frac{4-n}{2}}\left[1+\mathcal{O}(N)\right].$
(28)

To prove (a), we consider odd-$n$, and from Eq. (27) we get

$G_{n}^{(c)}=\begin{cases}0\quad&\text{if $n=1$}\\\ \infty\quad&\text{if
$n=3,5,\dots$}\end{cases}$ (29)

and therefore the condition (26a) for the occurrence of a second-order DPT is
never satisfied. If $U_{1},\gamma=0$, instead, Eq. (28) gives

$G_{n}^{(c)}=\begin{cases}0\quad&\text{if $n=1,3$}\\\ \infty\quad&\text{if
$n=5,7,\dots$}\end{cases}.$ (30)

We have therefore proven the statement (a).

Let us now consider the case of even $n$. From Eq. (27) we find that for
$U_{1}\neq 0$ or $\gamma\neq 0$ a second order DPT is possible only for $n=2$,
with a critical point given by

$G_{2}^{(c)}=\frac{\sqrt{4U_{1}^{2}+\gamma^{2}}}{4}.$ (31)

Higher $n$ results in $G_{n}^{(c)}=\infty$. Condition (26b) reads

$\left.\frac{\partial G_{2}(N)}{\partial
N}\right|_{G_{2}=G_{2}^{(c)}}=\frac{U_{1}U_{2}+4\gamma\eta_{2}}{2\sqrt{4U_{1}^{2}+\gamma^{2}}}$
(32)

and thus it can be satisfied for an appropriate choice of the parameters.
These equations proves (b) and (c).

Let us now assume $U_{1}=\gamma=0$. From Eq. (28) follows that, for $n=4$, a
second-order DPT can can place also for

$G_{4}^{(c)}=\frac{|U_{2}|}{4}.$ (33)

Therefore, condition (i) in Eq. (26a) is satisfied in the case
$U_{1}=\gamma=0$. As for condition (ii) in Eq. (26b), we have

$\left.\frac{\partial G_{4}(N)}{\partial
N}\right|_{G_{4}=G_{4}^{(c)}}=\frac{U_{3}}{4}\operatorname{Sign}\left(U_{2}\right),$
(34)

which can be satisfied choosing $U_{2}$ and $U_{3}$ with the same sign.
Finally, one can easily show that for $n>4$, Eq. (28) gives $G_{n}(N)=\infty$,
demonstrating the impossibility of a DPT, thus proving (e). ∎

Before dealing with the analysis of the full quantum results, let us remark
that, while $\gamma=0$ is impossible to achieve in actual realizations, for
many practical purposes one can consider system “sizes” $L$ where, to a
reasonable approximation, the role of $\gamma$ can be neglected, and thus the
approximation $\gamma=0$ faithfully recovers the results of finite-time
experiments. Furthermore, the detuning terms can be easily manipulated,
therefore making it possible to approximately fulfil the condition
$\gamma=U_{1}=0$ necessary to witness the second-order DPT for the $n=4$ case.

We also notice that the mechanism enabling second-order DPTs for $n=4$ (i.e.
the fact that $U_{2}$ and $U_{3}$ have the same sign) is the same behavior
displayed by the two-photon Kerr resonator, where this role is played by the
detuning $U_{1}$ and the two-photon interaction potential $U_{2}$ in Eq. (32).
Finally we stress that, although second-order DPTs could in principle emerge
also for even $n>4$, these would require setting $U_{2}=0$, that, contrarily
to detuning $U_{1}$, cannot be easily manipulated.

Figure 3: Multistability according to the semiclassical analysis in a four-
photon driven resonator ($n=4$), where we fixed $U_{1}=10\gamma$,
$U_{2}=-25\gamma$, $U_{3}=3\gamma$, $\eta_{4}=0.1\gamma$.

####  3.2.4 Multistability of solutions with different number of photons

The solution around $N=0$ predicts either a first- or a second-order phase
transition describing the passage of the system from the vacuum to a nonzero
population phase. This analysis does not predict the behavior far from $N=0$,
and nothing prevents several Hamiltonian terms from competing with each other,
thus resulting in multiple stable solution. In particular, the presence of
this multistability would imply an overlap of “S-like” curves of the
semiclassical solution, so that for the same drive intensity, there are
multiple solutions with different photon number.

To understand which mechanism can enable multistability, let us consider again
Eq. (20). In the semiclassical formalism, multistability implies the presence
of multiple solutions at the semiclassical level with different photon number
[c.f. Fig. 3 (b)]. This translates in the presence of multiple local minima
(or maxima) of the function $G_{n}(N)$, as shown in Fig. 3 (a). Therefore, one
can study the function

$\frac{\partial G_{n}(N)}{\partial N}=0.$ (35)

The number of maxima and minima signals the presence of multiple semiclassical
solutions. And following Descartes’ rule of signs – i.e., the maximal number
of positive roots of a polynomial is the number of sign changes between
consecutive coefficients – we deduce that, a necessary condition to have
multiple solutions, is the presence of alternating signs between the various
$U_{n}$. Physically speaking, the underlying mechanism is quite
straightforward: different $U_{n}$ terms can compete with each other in
determining the energy of one photon in the system, while drive and
dissipation favour the solution with more or less photons. Since the relevance
of each interaction term can change in different occupation regimes, several
solutions can emerge. The stability of the semiclassical solutions with
respect to quantum fluctuations need to be numerically assessed.

Figure 4: Onset of a first-order phase transition for increasing $L$ (see
legend) with symmetry breaking in the three-photon Kerr resonator. Panel (a):
mean number of photons in the steady state
$\braket{\hat{a}^{\dagger}\hat{a}}$, renormalised by scaling parameter $L$.
Panel (b): Real part of $\lambda_{1}^{(0)}$, i.e., the Liouvillian gap in the
same simmetry sector as the steady state, inducing the first-order transition.
The three veritical lines indicate the scaling values studied in Fig. 6. Panel
(c): Real part of $\lambda_{0}^{(1)}$, i.e., the Liovillian eigenvalue
signalling the spontaneous symmetry breaking. Parameters: $U_{1}/\gamma=-10$,
$U_{2}/\gamma=10$, $\eta_{3}/\gamma=1$.

##  4 Three-photon Kerr resonator

Having discussed the general properties of DPTs, we turn now to specific
examples, to demonstrate the validity of the semiclassical analysis, and show
the quantum properties around criticality. Throughtout the next two sections,
we will diagonalize the Liouvillian superoperator. We take full advantage of
the system’s symmetry, as detailed in Appendix B, to reduce the computational
complexity and enhance the precision of the results. For the most numerically
demanding simulations, we resort to the recently-developed Arnoldi-Lindblad
method [61], in conjunction with the algorithm detailed in Appendix B.

Here, we consider the three-photon driven Kerr resonator governed by the
master equation

$\partial_{t}\hat{\rho}(t)=-i\left[\hat{H}_{3},\hat{\rho}(t)\right]+\gamma\mathcal{D}[\hat{a}]+\eta_{3}\mathcal{D}[\hat{a}^{3}]$
(36)

with

$\hat{H}_{3}=U_{1}\hat{a}^{\dagger}\hat{a}+\frac{U_{2}}{2}\left(\hat{a}^{\dagger}\right)^{2}\hat{a}^{2}+G_{3}\left[\hat{a}^{3}+\left(\hat{a}^{\dagger}\right)^{3}\right].$
(37)

We focus on the $\gamma\neq 0$ case, where the system displays a $Z_{3}$ weak
symmetry. According to the semiclassical analysis, we expect a first-order
dissipative phase transition accompanied by the spontaneous breaking of the
weak $Z_{3}$ symmetry.

###  4.1 Semiclassical vs quantum solution

First, we analyze the photon number as a function of the driving strength
$G_{3}$. In Fig. 4 we show the results of the full quantum analysis (colored
lines), and compare them to the prediction of the semiclassical analysis
(dashed black line). For a weak drive $G_{3}$, the system is in the vacuum,
and $\hat{\rho}_{\rm ss}\simeq\ket{0}\bra{0}$. Increasing the drive intensity,
the system’s photon number deviates from the vacuum and approaches the high-
photon number branch predicted by the semiclassical theory. In this symmetry
broken phase, the stationary state is well-approximated by a statistical
mixture of three coherent states, i.e.,

$\displaystyle\hat{\rho}_{\rm
ss}\simeq\frac{\ket{\alpha_{1}}\bra{\alpha_{1}}+\ket{\alpha_{2}}\bra{\alpha_{2}}+\ket{\alpha_{3}}\bra{\alpha_{3}}}{3},$
(38)

where $\ket{\alpha_{1,2,3}}$ are coherent states with same the same number of
photons and a relative phase difference of $\pm 2\pi/3$, i.e.

$\alpha_{j+1}=\alpha_{j}\ e^{i\frac{2\pi}{3}}.$ (39)

The change in the steady state population becomes more and more abrupt as we
increase the parameters $L$, demonstrating that, indeed, the phase transition
is of the first order.

###  4.2 Analysis of the first-order transition

To confirm the presence of a first-order phase transition, we plot in Fig.
4(b) the Liouvillian eigenvalue $\lambda_{1}^{(0)}$ associated with the
slowest relaxation rate in the stady-state symmetry sector. This eigenvalue
signals hysteresis and critical slowing down, and the fact that it tends to
zero in the thermodynamic limit proves the presence of a first-order DPT [10].

We then investigate the properties of the eigenoperator $\hat{\rho}_{1}^{(0)}$
associated with such a state. According to Ref. [10], in the critical region
one can use the eigendecomposition of $\hat{\rho}_{1}^{(0)}$ to recast

$\hat{\rho}_{1}^{(0)}\simeq\hat{\rho}_{1}^{+}-\hat{\rho}_{1}^{-},$ (40)

where $\hat{\rho}_{1}^{\pm}$ represent the density matrices of the metastable
states. As such, we expect that, in the thermodynamic limit,
$\hat{\rho}_{1}^{\pm}$ recover the two stable solutions of the semiclassical
theory. We show the eigendecomposition in Fig. 5. We indeed find that the
semiclassical approximation qualitatively recovers the results of the
eigendecomposition.

Figure 5: Eigendecomposition of $\hat{\rho}_{1}^{(0)}$ for $L=15$ and the
parameters in Fig. 4. The solid black line describes the result of the full
quantum solution. The prediction of the semiclassical solitutions are the
dotted light blue curve (few photon numbers, stable), the dashed red curve
(high-photon number, stable), and the dashed black curve (unstable solution).
Finally, the results of the eigendecomposition are plotted with a solid blue
line and a dashed red line. Figure 6: Scaling towards the thermodynamic limit
for the three vertical lines in Fig. 4(b), demonstrating the presence of a
first-order DPT and the emergent stability of the vacuum.

As discussed in Sec. 3.2.1, the semiclassical analysis predict the presence of
a stable vacuum in the whole symmetry broken region. This analysis is
confirmed both by the eigendecomposition in Fig. 5 (the vacuum remains long-
lived even far from the transition point), and the from the spectral analysis
in Fig. 4(b). Considering larger values of $L$ results in slower timescales.
We confirm the scaling towards the thermodynamic limit of the Liouvillian gap
$\lambda_{1}^{(0)}$ in Fig. 6. We consider a point before the transition (red
line), at the minimum of the gap (blue line), and after the transition (green
line). The same three lines correspond to the vertical lines in Fig. 4(b). In
all the cases, after an initial transient, we see an exponential closure of
the gap as a function of $L$. The green curve confirms the vacuum
metastability predicted by the semiclassical theory.

###  4.3 Spontaneous symmetry breaking

The spontaneous symmetry breaking implies that, for strong enough pumping,
each of the state $\ket{\alpha_{i}}\bra{\alpha_{i}}$ in Eq. (38) becomes a
steady state of the system, since they are not eigenstate of $\mathcal{Z}_{3}$
[10]. We confirm this picture in Fig. 4(c), where we show that also
$\lambda_{0}^{(1)}$ becomes zero. As expected, we obtain an identical result
for $\lambda_{0}^{(2)}$ (not shown). This implies that the system coherent
state have become metastable.

Figure 7: Study of the strongly-symmetric four-photon Kerr resonator, and of
the onset of a second-order dissipative phase transition. For different values
of the thermodynamic rescaling parameter $L$: (a-d) photon number in the
$(j,j)$ symmetry sector; (e-h) Liouvillian gap in the $(j,j)$ symmetry sector;
(e-h) smallest Liouvillian eigenvalue in the $(j,j+1)$ symmetry sector.
Parameters: $\gamma=U_{1}=0$, $U_{2}=10\eta_{4}$, $U_{3}=\eta_{4}$.

##  5 Four-photon Kerr resonator

Here, we consider the four-photon driven-dissipative Kerr resonator, reading

$\partial_{t}\hat{\rho}(t)=-i\left[\hat{H}_{4},\hat{\rho}(t)\right]+\gamma\mathcal{D}[\hat{a}]+\eta_{4}\mathcal{D}[\hat{a}^{4}]$
(41)

with

$\begin{split}\hat{H}_{4}=U_{1}\hat{a}^{\dagger}\hat{a}&+\frac{U_{2}}{2}\left(\hat{a}^{\dagger}\right)^{2}\hat{a}^{2}+\frac{U_{3}}{3}\left(\hat{a}^{\dagger}\right)^{3}\hat{a}^{3}\\\
&+G_{4}\left[\hat{a}^{4}+\left(\hat{a}^{\dagger}\right)^{4}\right].\end{split}$
(42)

###  5.1 Strong symmetry and second-order phase transition

We start by considering the strong symmetric case $\gamma=0$ with $U_{1}=0$.
For this set of parameter, the semiclassical analysis predicts a second-order
phase transition associated with the spontaneous breaking of a $Z_{4}$ strong
symmetry. We analyze it in Fig. 7. We recall that, since the system has a
strong $Z_{4}$ symmetry, the number of Liouvillian sectors is $4\times 4$,
being characterized by the two quantum numbers $k_{\rm L}$ and $k_{\rm R}$.
The $4$ sectors with $k_{\rm L}=k_{\rm R}$ describe the evolution of the
populations, while the remaining $12$ with $k_{\rm L}\neq k_{\rm R}$ describe
the evolution of the coherences.

First, we consider the re-scaled photon number of the steady state for each of
the symmetry sectors $(j,j)$ with $j\in[0,3]$, and increasing the
thermodynamic parameter $L$. Calling $\hat{\rho}_{\rm
ss}^{(j)}\propto\hat{\rho}_{0}^{(j,j,)}$ the steady state in each symmetry
sector,
$\langle\hat{a}^{\dagger}\hat{a}\rangle_{j}=\mathrm{Tr}\\!\left[\hat{a}^{\dagger}\hat{a}\hat{\rho}_{\rm
ss}^{j}\right]$ are plotted in Figs. 7(a-d). For low drive amplitudes, the
system is in the $Z_{n}$ symmetric vacuum. Indeed, the states need to respect
the strong symmetry condition in Eq. (13), and thus

$\hat{\rho}_{\rm ss}^{(j)}=\ket{{\rm vac}_{j}}\bra{{\rm
vac}_{j}}=\ket{j}\bra{j},$ (43)

where $j$ labels the symmetry sector and $\ket{j}$ is the Fock state with $j$
photons. For large drive, instead, the system transition towards

$\hat{\rho}_{\rm ss}^{(j)}\simeq\ket{\mathcal{K}_{j}}\bra{\mathcal{K}_{j}}$
(44)

where the Schrödinger cats $\ket{\mathcal{K}_{i}}$ are

$\ket{\mathcal{K}_{j}}=\frac{1}{\mathcal{N}}\sum_{n=0}^{3}e^{i\pi
jn/2}\ket{\alpha_{n}},$ (45)

where $\ket{\alpha_{n}}$ are coherent states such that $\alpha_{n}=e^{i\pi
n/2}\alpha$ and $\mathcal{N}$ is a normalization factor. Increasing the value
of $L$ towards the thermodynamic limit, we observe that the passage between
the $Z_{n}$ vacua and the cat states becomes sharper and sharper, but remains
continuous. This analysis corroborate the semiclassical one, and by
appropriately taking into account the system’s symmetry, we observe a second-
order DPT. To further demonstrate that, indeed, the transition is of the
second and not of the first order, we plot $\lambda_{1}^{(j,j)}$ in Figs.
7(e-h), i.e., the Liouvillian gap of the $(j,j)$ symmetry sector. We observe
no closure of the Liouvillian gap, indicating that no critical slowing down or
hysteresis occurs for the Liouvillian populations.

Finally, we plot the smallest Liouvillian eigenvalue $\lambda_{0}^{(j,j+1)}$
for the sectors $(j,j+1)$ (where $j+1=0$ if $j=3$) in Figs. 7(i-l). These
represent the decay rate of coherences between the sector $j$ and $j+1$, and
their closure indicate the possibility of retaining everlasting coherences. In
this case, we observe that, after the critical point, these eigenvalue
progressively become smaller, indicating that the system undergoes a second
order phase transition. We obtain similar results for the other $(j,k)$
sectors with $j\neq k$ (not shown). This is associated with a spontaneous
breaking of the strong $Z_{4}$ symmetry, because it results in

$\begin{split}&\hat{\rho}\propto\left(\ket{\mathcal{K}_{j}}+\ket{\mathcal{K}_{k}}\right)\left(\bra{\mathcal{K}_{j}}+\bra{\mathcal{K}_{k}}\right),\\\
&\mathcal{L}\hat{\rho}=0\quad{\rm but}\quad\mathcal{Z}_{4}^{\rm
L,R}\hat{\rho}\neq z_{4}^{\rm L,R}\hat{\rho}.\end{split}$ (46)

###  5.2 Weak symmetry and multistability

Figure 8: Analysis of the classically multistable system. (a) Photon number as
a function of the drive for different values of the thermodynamic scaling
parameter $L$. The black dashed line represents the semiclassical solution.
(b) Liouvillian gap and (c) second Liouvillian eigenvalue in the $k=0$ sector,
demonstrating the presence of two slow timescales. (d) Smallest (e) and second
smalles Liouvillian eigenvalues in the $k=1$ sector, demonstrating the
presence of SSB and of a slow timescale. Parameters: $U_{1}=10\gamma$,
$U_{2}/\gamma=-25\gamma$, $U_{3}=3\gamma$, $\eta_{4}=0.1\gamma$.

We now consider a weakly symmetric case in the presence of detuning $U_{1}$,
and with competing terms giving rise to multistability according to the
semiclassical solution. First, in Fig. 8(a), we compare the results of the
semiclassical analysis with those of the full quantum simulation. We find
that, although the semiclassical solution has three stable solution, the full
quantum simulation is characterized by a single first-order DPT, from the
vacuum to the highest-populated manifold. Indeed, if we analyze the
Liouvillian gap $\lambda_{1}^{(0)}$ in Fig. 8(b) we clearly see the closure of
the Liouvillian gap associated with a first-order DPT. If, however, we also
consider the second eigenvalue $\lambda_{2}^{(0)}$ as in Fig. 8(c), we see
that a second slow timescale emerges. That is, despite the presence of a
single phase transition, the dynamics of the population of the system is
characterized by two slow timescales.

We corroborate this phenomenon by analyzing the symmetry sectors responsible
for spontaneous symmetry breaking. In Fig. 8(c), we plot $\lambda_{0}^{(1)}$
showing that, indeed, this phenomenon is accompanied by the breaking of the
weak $Z_{4}$ symmetry. Noticeably, the spontaneous symmetry breaking takes
place before the occurence of the first-order transition. Furthermore, we also
observe a second slow timescale for this symmetry sector, i.e.,
$\lambda_{1}^{(1)}$ in Fig. 8(d). This slow timescales represent the fact that
there exist multiple symmetry broken states, and there is a slow rate from
which the system switches between them. We observe similar results for the
other symmetry sectors (not shown).

The picture we derive is one in which, although there are only two real steady
state of the dynamics, either the vacuum or the one at large photon number,
there exist a third metastable state to which the system can be initialized.
Such a state is characterized by a broken symmetry, but it cannot be reached
by quantum fluctuation alone.

To further demonstrate this picture, in Fig. 9(a) we use the
eigendecomposition to express the eigenoperators associated with the slowest
eigenvalues as

$\displaystyle\hat{\rho}_{1}^{(0)}=\hat{\rho}_{1}^{(0),+}-\hat{\rho}_{1}^{(0),-},\quad\hat{\rho}_{2}^{(0)}=\hat{\rho}_{2}^{(0),+}-\hat{\rho}_{2}^{(0),-}.$
(47)

As one can see, these metastables density matrices recover the results of the
semiclassical analysis, and the region in which there is a closure of these
Liouvillian eigenvalues roughly corresponds to the region of multistability
according to the semiclassical analysis.

Overall, the system displays 9 metastable coherent-like states, approximated
by $\ket{\alpha_{\rm vac}}$, $\ket{\alpha_{\rm low}e^{i\phi_{j}}}$, and
$\ket{\alpha_{\rm high}e^{i\phi_{j}}}$, with $|\alpha_{\rm vac}|<|\alpha_{\rm
low}|<|\alpha_{\rm high}|$ and $\phi_{j}\in j\pi/4$.

Figure 9: Eigendecomposition and comparison with the semiclassical solution.
(a) Photon number of the full quantum solution (black solid line) compared to
the semiclassical solution (black dashed line) and the results of the
eigendecomposition [red and blue markers correspond to the red and blue curves
in panel (b)]. (b) The two smallest Liouvillian eigenvalues, whose behavior
across the transition has been reconstructed using the continuity of the
associated eigenoperators. Parameters as in Fig. 8 for $L=20$.

##  6 Conclusions and outlook

In this work we explored the critical properties of $n$-photon driven-
dissipative nonlinear quantum resonators. We found that the symmetries of the
model, fixed by driving and dissipation, determines the nature of the phase
transitions in the steady state. We characterize such criticalities providing
general results for this class of models.

We attack the problem using a semiclassical approach valid in a well defined
thermodynamic limit with an infinite number of excitations. In such a limit
the state of the system approaches a coherent state and quantum fluctuations
are suppressed leading to a generalized version of the driven-dissipative
Gross-Pitaevskii equation. Studying its stationary properties we formulate and
prove a no-go theorem stating that no second-order phase transitions are
possible when $n$ is odd, while, for even $n$, second-order transitions can
take place only for $n=2$ and $n=4$.

We then perform a full quantum analysis of the three- and four-photon driven
Kerr resonators. We find that quantum fluctuations trigger the transition
between semiclassical solutions in the thermodynamic limit validating the
results obtained in the semiclassical limit. While the semiclassical
approximation has been proved to be reliable for $n=1,2$, for higher $n$ there
are no strong arguments supporting its validity. Indeed the systematic
inclusion of small quantum fluctuation on top of the mean-field semiclassical
solution can be obtained via truncated Wigner methods [62, 63] and Gaussian
expansions [49]. This is not the case for $n>2$ because the drive and
$n$-photon dissipation can, in principle, introduces non-Gaussian correlation
above the coherent-state solution [64]. The emergence of this dissipative
phase transitions is understood and characterized within the spectral theory
of Liouvillian highlighting the role of weak and strong symmetries.

These results could also be relevant in field of quantum technologies and
quantum information encoding. Symmetry breaking in second-order DPTs has been
demonstrated to be a resource to improve the sensitivity of quantum
measurement protocols [34, 38]. Our work proves that such kind of enhancement
can only be attained for $n=2$ or $n=4$. Furthermore, our results may pose
constraints for the exploitation of nonlinear driven resonators for the
encoding of bosonic codes. As it been recently proposed [30, 28], detuning and
critical phenomena may play a key role in the storing quantum information. The
metastability of the vacuum may also prove an obstacle towards a rapid and
reliable initialization of bosonic qubits.

This work paves the way to future intriguing research directions. Among them,
we mention the study of the dynamical properties of these systems in
connection with quantum trajectories approaches, and the emergence of chaotic
behavior in highly nonlinear quantum resonators.

## Acknowledgements

We thank G. Rastelli and L. Gravina for the useful discussions. We acknowledge
the help of A. Mercurio in the optimization of the numerical codes. This work
was supported by the Swiss National Science Foundation through Project No.
200020_185015, Provincia Autonoma di Trento, and was conducted with the
financial support of the EPFL Science Seed Fund 2021 and from PNRR MUR project
PE0000023-NQSTI.

##  Appendix A Interaction, nonlinearities, and $n$-photon drives in
superconducting circuits

Let us consider a standard LC resonator characterized by the Hamiltonian:

$\hat{H}_{\rm cav}=\omega\hat{a}^{\dagger}\hat{a},$ (48)

where $\omega$ is the resonator frequency.

Non-quadratic (i.e., interaction) terms can emerge by considering the action
of nonlinear elements. For instance, nonlinearity can be obtained by
quantizing the flux in Josephson junctions potentials of the form

$E_{J}\cos\left(\frac{\phi}{\phi_{0}}\right)\simeq
E_{J}\left(1-\frac{\phi^{2}}{2\phi_{0}^{2}}+\frac{\phi^{4}}{24\phi_{0}^{4}}-\frac{\phi^{6}}{720\phi_{0}^{6}}\right)+\dots,$
(49)

where $\phi$ is the flux coordinate of the circuit at the junction and
$\phi_{0}$ is the magnetic flux quantum. In most implementations, the
expansion in Eq. (49) can be stopped at the $\phi^{4}$ order. If the Josephson
junction belongs to a single resonator, by substituting
$\phi/\phi_{0}\propto\hat{a}+\hat{a}^{\dagger}$, and discarding counter-
rotating terms, which are out of resonance, one obtains the Kerr resonator
Hamiltonian, reading

$\hat{H}_{\rm
Kerr}=\tilde{U}_{1}\hat{a}^{\dagger}\hat{a}+\frac{U_{2}}{2}\left(\hat{a}^{\dagger}\right)^{2}\left(\hat{a}\right)^{2}.$
(50)

In $n$-photon driven systems, photons are coherently exchanged between the
resonator and a set of external fields, $n$ at the time. While single-photon
drive (i.e., of the form $\hat{a}+\hat{a}^{\dagger}$) can emerge by, e.g.,
capacitive coupling an incoming wave-guide with the cavity, higher order drive
require to be mediated by nonlinear elements. For instance, such $n$-drive
terms can be derived from the expansion in Eq. (49) if the Josephson junction
is shared by several modes, so that $\phi=\sum_{k}\phi_{k}$, where $\phi_{k}$
represents the flux coordinate of each one of the modes. For instance, two-
photon drives can be achieved by standard four waves mixing, rewriting
$\phi=\phi_{a}+\phi_{b}+\phi_{c}$, where
$\phi_{a}\propto\hat{a}+\hat{a}^{\dagger}$ is the field within the resonator,
and $\phi_{b}\propto\hat{b}+\hat{b}^{\dagger}$ and
$\phi_{c}\propto\hat{c}+\hat{c}^{\dagger}$ are auxiliary modes. If the mode
$b$ ($c$) is driven, and evolves on a timescale much faster than the typical
time scales of the $a$ mode, one can substitute the operators $\hat{b}$
($\hat{c}$) with a c-number oscillating at the driving frequency $\omega_{b}$
($\omega_{c}$) via an adiabatic elimination, reading $\hat{b}\to
be^{i\omega_{b}t}$ ($\hat{c}\to ce^{i\omega_{c}t}$ ). All in all, discarding
again out-of-resonance terms, the Hamiltonian for the resonator resulting from
the fourth-order expansion of the potential $\cos(\phi)$ would result in a
nonlinear Hamiltonian $\hat{H}_{\rm NL}$, reading

$\hat{H}_{\rm NL}=\hat{H}_{\rm
Kerr}+G_{2}\left[\hat{a}^{2}e^{2i\omega_{p}t}+\left(\hat{a}^{\dagger}\right)^{2}e^{-2i\omega_{p}t}\right].$
(51)

By passing in the frame rotating at the drive frequency, and re-absorbing the
contribution to the energy frequency in the term
$U_{1}=\tilde{U}_{1}-\omega_{p}=-\Delta$, the Hamiltonian finally reads

$\hat{H}_{n=2}=U_{1}\hat{a}^{\dagger}\hat{a}+\frac{U_{2}}{2}\left(\hat{a}^{\dagger}\right)^{2}\hat{a}^{2}+G_{2}\left[\hat{a}^{2}+\left(\hat{a}^{\dagger}\right)^{2}\right].$
(52)

Through similar procedures, high-order expansion of nonlinear terms ($k$-wave
mixing with $k>n$), can result (in principle) in $n$-photon drives. By
including such terms, one needs to include also the corresponding
nonlinearities. As detailed in the main text, such nonlinearities can play a
fundamental role in determining the nature of the transition.

##  Appendix B An efficient algorithm for block-diagonalizing the Liovuillian
in the presence of $Z_{n}$ symmetries

We introduce here a simple algorithm to block diagonalize system displaying a
$Z_{n}$ symmetry. Although we describe it for a weakly symmetric case, its
extension to a strong symmetry is straightforward.

The Liouvillian admits an abstract definition of its spectrum via Eq. (9). To
numerically obtain the eigenvalues and eigenoperators one needs to explicit
the matrix form of $\mathcal{L}$. For a finite-dimensional Hilbert space, one
can construct such a matrix via

$\begin{split}\mathcal{L}&=-i\left(\hat{H}\otimes\hat{\mathds{1}}-\hat{\mathds{1}}\otimes\hat{H}^{\mathrm{T}}\right)\\\
&\quad+\sum_{j=1}^{3}\left(\hat{L}_{j}\otimes\hat{L}^{*}_{j}-\frac{\hat{L}_{j}^{\dagger}\hat{L}_{j}\otimes\hat{\mathds{1}}+\hat{\mathds{1}}\otimes\hat{L}_{j}^{\rm
T}\hat{L}_{j}^{*}}{2}\right),\end{split}$ (53)

where $\hat{L}_{j}^{\rm T}$ represents the transpose of $\hat{L}_{j}$. The
spectrum of the Liouvillian can then be directly obtained by diagonalizing the
matrix representation of $\mathcal{L}$. For infinite dimensional spaces (i.e.,
those of bosonic systems), one need to introduce a cutoff in the Hilbert space
$N_{c}$. That is, one projects the true infinite-dimensional Hamiltonian and
jump operators onto the space spanned by the Fock states $\ket{n}$ for
$n\in[0,N_{c})$, and assumes that the matrix elements of any operator for
$n\in[N_{c},\infty)$ are zero.

Since $[\mathcal{Z}_{n},\mathcal{L}]=0$, all the $\hat{\rho}_{i}$ are
eigenoperators of $\mathcal{Z}_{N}$. And since $\mathcal{Z}_{n}$ admits $n$
different eigenvalues, it is always possible to block-diagonalize the
Liouvillian into (at least) $n$ smaller blocks. Each block $\mathcal{L}_{n}$
describes completely the physics of each symmetry sector of the full
Liouvillian $\mathcal{L}$. Normally, to put the Liouvillian in its block-
diagonal form one would construct the basis of a symmetry sector determining
the eigenoperators $\hat{\zeta}_{i}$ of $\mathcal{Z}_{n}$ and project the
Liouvillian onto the correct symmetry sector obtaining the matrix elements

$\mathcal{L}_{i,j}=\operatorname{Tr}\left[\hat{\zeta}_{i}^{\dagger}\left(\mathcal{L}\hat{\zeta}_{j}\right)\right].$
(54)

Even if in principle correct, this process is extremely slow and inefficient
since the Liouvillian is a very sparse and large matrix.

Instead of applying Eq. (54), we notice that the Fock basis is already the
basis of eigenstates of $\mathcal{Z}_{n}$, as it follows from Eq. (12). That
is, when using Eq. (53), we are using the correct basis to obtain the block-
diagonal form of the Liouvillian, simply we are considering the basis in the
wrong order. Hence, the Liouvillian is a permutation of row and columns away
from being block diagonal, and the algorithm that we seek is one which
efficiently finds the correct permutation matrix $\mathcal{P}$ which
transforms $\mathcal{L}$ into its block diagonal form, whenever such a
transformation is possible.

The main idea is to model the block diagonalization problem as an equivalent
graph-theoretic problem.

  1. 1.

$\mathcal{L}$ is written as the adjacency matrix of an undirected graph;

  2. 2.

Each block in the block diagonal form is a single connected component in the
graph; thus, the problem boils down to finding each connected component in the
graph.

  3. 3.

We then use the Breadth First/Depth First search algorithm. consecutively to
obtain the permutation matrices and the indices of the blocks. The time to
perform this task (i.e., its computational complexity) is linear in the number
of nodes in the graph.

  4. 4.

We use the permutation matrix to produce each block $\mathcal{L}_{i}$ such
that
${\mathcal{L}=\mathcal{P}\operatorname{diag}\\{\mathcal{L}_{1}\dots\mathcal{L}_{n}\\}\mathcal{P}^{\rm
T}}$.

The key factor in the numerical speedup comes from the fact that obtaining the
permutation matrix $\mathcal{P}$ requires a number of operations linear in the
number of nonzero elements of the Liouvillian, which is a very sparse matrices
[c.f. Eq. (53)].

## References

  * [1] I. Carusotto and C. Ciuti,  _Quantum fluids of light_ , Rev. Mod. Phys. 85, 299.
  * Carusotto  _et al._ [2020] I. Carusotto, A. A. Houck, A. J. Kollár, P. Roushan, D. I. Schuster and J. Simon,  _Photonic materials in circuit quantum electrodynamics_ , Nat. Phys. 16, 268 (2020).
  * Hur  _et al._ [2016] K. L. Hur, L. Henriet, A. Petrescu, K. Plekhanov, G. Roux and M. Schiró,  _Many-body quantum electrodynamics networks: Non-equilibrium condensed matter physics with light_ , C. R. Phys. 17, 808 (2016).
  * Breuer and Petruccione [2007] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf and J. I. Cirac,  _Quantum computation and quantum-state engineering driven by dissipation_ , Nat. Phys. 5, 633 (2009).
  * Diehl  _et al._ [2008] S. Diehl, A. Micheli, A. Kantian, B. Kraus, H. P. Büchler and P. Zoller,  _Quantum states and phases in driven open quantum systems with cold atoms_ , Nat. Phys. 4, 878 (2008).
  * Diehl  _et al._ [2010] S. Diehl, A. Tomadin, A. Micheli, R. Fazio and P. Zoller,  _Dynamical Phase Transitions and Instabilities in Open Atomic Many-Body Systems_ , Phys. Rev. Lett. 105, 015702 (2010).
  * Buča and Prosen [2012] B. Buča and T. Prosen,  _A note on symmetry reductions of the Lindblad equation: transport in constrained open spin chains_ , New J. Phys. 14, 073007 (2012).
  * Albert and Jiang [2014] V. V. Albert and L. Jiang,  _Symmetries and conserved quantities in Lindblad master equations_ , Phys. Rev. A 89, 022118 (2014).
  * Minganti  _et al._ [2018] F. Minganti, A. Biella, N. Bartolo and C. Ciuti,  _Spectral theory of Liouvillians for dissipative phase transitions_ , Phys. Rev. A 98, 042118 (2018).
  * Bartolo  _et al._ [2016] N. Bartolo, F. Minganti, W. Casteels and C. Ciuti,  _Exact steady state of a Kerr resonator with one- and two-photon driving and dissipation: Controllable Wigner-function multimodality and dissipative phase transitions_ , Phys. Rev. A 94, 033841 (2016).
  * Lebreuilly  _et al._ [2017] J. Lebreuilly, A. Biella, F. Storme, D. Rossini, R. Fazio, C. Ciuti and I. Carusotto,  _Stabilizing strongly correlated photon fluids with non-Markovian reservoirs_ , Phys. Rev. A 96, 033828 (2017).
  * Biella  _et al._ [2017] A. Biella, F. Storme, J. Lebreuilly, D. Rossini, R. Fazio, I. Carusotto and C. Ciuti,  _Phase diagram of incoherently driven strongly correlated photonic lattices_ , Phys. Rev. A 96, 023839 (2017).
  * Leghtas  _et al._ [2015] Z. Leghtas, S. Touzard, I. M. Pop, A. Kou, B. Vlastakis, A. Petrenko, K. M. Sliwa, A. Narla, S. Shankar, M. J. Hatridge et al.,  _Confining the state of light to a quantum manifold by engineered two-photon loss_ , Science 347, 853 (2015).
  * Grimm  _et al._ [2020] A. Grimm, N. E. Frattini, S. Puri, S. O. Mundhada, S. Touzard, M. Mirrahimi, S. M. Girvin, S. Shankar and M. H. Devoret,  _Stabilization and operation of a Kerr-cat qubit_ , Nature 584, 205 (2020).
  * Mirrahimi  _et al._ [2014] M. Mirrahimi, M. Leghtas, V. Albert, S. Touzard, R. Schoelkopf, L. Jiang and M. Devoret,  _Dynamically protected cat-qubits: a new paradigm for universal quantum computation_ , New Journal of Physics 16, 045014 (2014).
  * Chan  _et al._ [2008] H. B. Chan, M. I. Dykman and C. Stambaugh,  _Paths of Fluctuation Induced Switching_ , Phys. Rev. Lett. 100, 130602 (2008).
  * Leuch  _et al._ [2016] A. Leuch, L. Papariello, O. Zilberberg, C. L. Degen, R. Chitra and A. Eichler,  _Parametric Symmetry Breaking in a Nonlinear Resonator_ , Phys. Rev. Lett. 117, 214101 (2016).
  * Bartolo  _et al._ [2017] N. Bartolo, F. Minganti, J. Lolli and C. Ciuti,  _Homodyne versus photon-counting quantum trajectories for dissipative Kerr resonators with two-photon driving_ , Eur. Phys. J. Spec. Top. 226, 2705 (2017).
  * Goto [2016] H. Goto,  _Universal quantum computation with a nonlinear oscillator network_ , Phys. Rev. A 93, 050301 (2016).
  * Labay-Mora  _et al._ [2022] A. Labay-Mora, R. Zambrini and G. L. Giorgi,  _Quantum associative memory with a single driven-dissipative non-linear oscillator_ ,  (2022), 10.48550/ARXIV.2205.09491.
  * Landa  _et al._ [2020] H. Landa, M. Schiró and G. Misguich,  _Multistability of Driven-Dissipative Quantum Spins_ , Phys. Rev. Lett. 124, 043601 (2020).
  * Kessler  _et al._ [2012] E. M. Kessler, G. Giedke, A. Imamoǧlu, S. F. Yelin, M. D. Lukin and J. I. Cirac,  _Dissipative phase transition in a central spin system_ , Phys. Rev. A 86, 012116 (2012).
  * Casteels  _et al._ [2016] W. Casteels, F. Storme, A. Le Boité and C. Ciuti,  _Power laws in the dynamic hysteresis of quantum nonlinear photonic resonators_ , Phys. Rev. A 93, 033824 (2016).
  * Rodriguez  _et al._ [2017] S. R. K. Rodriguez, W. Casteels, F. Storme, N. Carlon Zambon, I. Sagnes, L. Le Gratiet, E. Galopin, A. Lemaître, A. Amo, C. Ciuti et al.,  _Probing a Dissipative Phase Transition via Dynamical Optical Hysteresis_ , Phys. Rev. Lett. 118, 247402 (2017).
  * Savona [2017] V. Savona,  _Spontaneous symmetry breaking in a quadratically driven nonlinear photonic lattice_ , Phys. Rev. A 96, 033826 (2017).
  * Rota  _et al._ [2019] R. Rota, F. Minganti, C. Ciuti and V. Savona,  _Quantum Critical Regime in a Quadratically Driven Nonlinear Photonic Lattice_ , Phys. Rev. Lett. 122, 110405 (2019).
  * Lieu  _et al._ [2020] S. Lieu, M. McGinley and N. R. Cooper,  _Tenfold Way for Quadratic Lindbladians_ , Phys. Rev. Lett. 124, 040401 (2020).
  * Halati  _et al._ [2022] C.-M. Halati, A. Sheikhan and C. Kollath,  _Breaking strong symmetries in dissipative quantum systems: Bosonic atoms coupled to a cavity_ , Phys. Rev. Res. 4, L012015 (2022).
  * Gravina  _et al._ [2022] L. Gravina, F. Minganti and V. Savona,  _A critical Schrödinger cat qubit_ , arXiv:2208.04928 (2022).
  * Fernández-Lorenzo and Porras [2017] S. Fernández-Lorenzo and D. Porras,  _Quantum sensing close to a dissipative phase transition: Symmetry breaking and criticality as metrological resources_ , Phys. Rev. A 96, 013817 (2017).
  * Ilias  _et al._ [2022] T. Ilias, D. Yang, S. F. Huelga and M. B. Plenio,  _Criticality-Enhanced Quantum Sensing via Continuous Measurement_ , PRX Quantum 3, 010354 (2022).
  * Raghunandan  _et al._ [2018] M. Raghunandan, J. Wrachtrup and H. Weimer,  _High-Density Quantum Sensing with Dissipative First Order Transitions_ , Phys. Rev. Lett. 120, 150501 (2018).
  * Di Candia  _et al._ [2021] R. Di Candia, F. Minganti, K. Petrovnin, G. Paraoanu and S. Felicetti,  _Critical parametric quantum sensing_ , arXiv:2107.04503 (2021).
  * Takemura  _et al._ [2021] N. Takemura, M. Takiguchi and M. Notomi,  _Low- and high- $\beta$ lasers in the class-A limit: photon statistics, linewidth, and the laser-phase transition analogy_, J. Opt. Soc. Am. B 38, 699 (2021).
  * Minganti  _et al._ [2021] F. Minganti, I. I. Arkhipov, A. Miranowicz and F. Nori,  _Liouvillian spectral collapse in the Scully-Lamb laser model_ , Phys. Rev. Res. 3, 043197 (2021).
  * Yacomotti  _et al._ [2022] A. M. Yacomotti, Z. Denis, A. Biella and C. Ciuti,  _Quantum Density Matrix Theory for a Laser Without Adiabatic Elimination of the Population Inversion: Transition to Lasing in the Class-B Limit_ , Laser & Photonics Reviews 17, 2200377 (2022).
  * Heugel  _et al._ [2019] T. L. Heugel, M. Biondi, O. Zilberberg and R. Chitra,  _Quantum Transducer Using a Parametric Driven-Dissipative Phase Transition_ , Phys. Rev. Lett. 123, 173601 (2019).
  * Minganti  _et al._ [2016] F. Minganti, N. Bartolo, J. Lolli, W. Casteels and C. Ciuti,  _Exact results for Schrödinger cats in driven-dissipative systems and their feedback control_ , Scientific Reports 6, 26987 (2016).
  * Roberts and Clerk [2020] D. Roberts and A. A. Clerk,  _Driven-Dissipative Quantum Kerr Resonators: New Exact Solutions, Photon Blockade and Quantum Bistability_ , Phys. Rev. X 10, 021022 (2020).
  * Zhang and Baranger [2021] X. H. H. Zhang and H. U. Baranger,  _Driven-dissipative phase transition in a Kerr oscillator: From semiclassical $\mathcal{PT}$ symmetry to quantum fluctuations_, Phys. Rev. A 103, 033711 (2021).
  * Fitzpatrick  _et al._ [2017] M. Fitzpatrick, N. M. Sundaresan, A. C. Y. Li, J. Koch and A. A. Houck,  _Observation of a Dissipative Phase Transition in a One-Dimensional Circuit QED Lattice_ , Phys. Rev. X 7, 011016 (2017).
  * Fink  _et al._ [2018] T. Fink, A. Schade, S. Höfling, C. Schneider and A. Imamoǧlu,  _Signatures of a dissipative phase transition in photon correlation measurements_ , Nat. Phys. 14, 365 (2018).
  * Brookes  _et al._ [2021] P. Brookes, G. Tancredi, A. D. Patterson, J. Rahamim, M. Esposito, T. K. Mavrogordatos, P. J. Leek, E. Ginossar and M. H. Szymanska,  _Critical slowing down in circuit quantum electrodynamics_ , Science Advances 7 (2021), 10.1126/sciadv.abe9492.
  * Chen  _et al._ [2022] Q.-M. Chen, M. Fischer, Y. Nojiri, M. Renger, E. Xie, M. Partanen, S. Pogorzalek, K. G. Fedorov, A. Marx, F. Deppe et al.,  _Quantum behavior of a superconducting Duffing oscillator at the dissipative phase transition_ ,  (2022), 10.48550/ARXIV.2206.06338.
  * Drummond and Walls [1980] P. D. Drummond and D. F. Walls,  _Quantum theory of optical bistability. I. Nonlinear polarisability model_ , J. Phys. A: Math. Theor. 13, 725 (1980).
  * Vicentini  _et al._ [2018] F. Vicentini, F. Minganti, R. Rota, G. Orso and C. Ciuti,  _Critical slowing down in driven-dissipative Bose-Hubbard lattices_ , Phys. Rev. A 97, 013853 (2018).
  * Foss-Feig  _et al._ [2017] M. Foss-Feig, P. Niroula, J. T. Young, M. Hafezi, A. V. Gorshkov, R. M. Wilson and M. F. Maghrebi,  _Emergent equilibrium in many-body optical bistability_ , Phys. Rev. A 95, 043826 (2017).
  * Verstraelen  _et al._ [2020] W. Verstraelen, R. Rota, V. Savona and M. Wouters,  _Gaussian trajectory approach to dissipative phase transitions: The case of quadratically driven photonic lattices_ , Phys. Rev. Res. 2, 022037 (2020).
  * Rota and Savona [2019] R. Rota and V. Savona,  _Simulating frustrated antiferromagnets with quadratically driven QED cavities_ , Phys. Rev. A 100, 013838 (2019).
  * Svensson  _et al._ [2018] I.-M. Svensson, A. Bengtsson, J. Bylander, V. Shumeiko and P. Delsing,  _Period multiplication in a parametrically driven superconducting resonator_ , Applied Physics Letters 113, 022602 (2018).
  * Chang  _et al._ [2020] C. W. S. Chang, C. Sabín, P. Forn-Díaz, F. Quijandría, A. M. Vadiraj, I. Nsanzineza, G. Johansson and C. M. Wilson,  _Observation of Three-Photon Spontaneous Parametric Down-Conversion in a Superconducting Parametric Cavity_ , Phys. Rev. X 10, 011011 (2020).
  * Lang and Armour [2021] B. Lang and A. D. Armour,  _Multi-photon resonances in Josephson junction-cavity circuits_ , New Journal of Physics 23, 033021 (2021).
  * Carmichael [2007] H. Carmichael,  _Statistical Methods in Quantum Optics 2: Non-Classical Fields_ (Springer, Berlin, 2007).
  * Rivas and Huelga [2011] Á. Rivas and S. F. Huelga,  _Open Quantum Systems: An Introduction_ (Springer, Berlin, 2011).
  * Casteels  _et al._ [2017] W. Casteels, R. Fazio and C. Ciuti,  _Critical dynamical properties of a first-order dissipative phase transition_ , Phys. Rev. A 95, 012128 (2017).
  * Casteels and Ciuti [2017] W. Casteels and C. Ciuti,  _Quantum entanglement in the spatial-symmetry-breaking phase transition of a driven-dissipative Bose-Hubbard dimer_ , Phys. Rev. A 95, 013812 (2017).
  * Peng  _et al._ [2019] J. Peng, E. Rico, J. Zhong, E. Solano and I. L. Egusquiza,  _Unified superradiant phase transitions_ , Phys. Rev. A 100, 063820 (2019).
  * Felicetti and Le Boité [2020] S. Felicetti and A. Le Boité,  _Universal Spectral Features of Ultrastrongly Coupled Systems_ , Phys. Rev. Lett. 124, 040404 (2020).
  * Hwang  _et al._ [2018] M.-J. Hwang, P. Rabl and M. B. Plenio,  _Dissipative phase transition in the open quantum Rabi model_ , Phys. Rev. A 97, 013825 (2018).
  * Minganti and Huybrechts [2022] F. Minganti and D. Huybrechts,  _Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of time-independent and Floquet open quantum systems_ , Quantum 6, 649 (2022).
  * Risken and Vollmer [1967] H. Risken and H. D. Vollmer,  _The influence of higher order contributions to the correlation function of the intensity fluctuation in a Laser near threshold_ , Z. Physik 201, 323 (1967).
  * Risken  _et al._ [1987] H. Risken, C. Savage, F. Haake and D. F. Walls,  _Quantum tunneling in dispersive optical bistability_ , Phys. Rev. A 35, 1729 (1987).
  * Minganti  _et al._ [2023] F. Minganti, L. Garbe, A. Le Boité and S. Felicetti,  _Non-Gaussian superradiant transition via three-body ultrastrong coupling_ , Phys. Rev. A 107, 013715 (2023).