Dissipative phase transitions in n -photon driven quantum nonlinear resonators

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 and allows 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 solutions emerge around zero, low, and high-photon numbers. Our results highlight the crucial role played by strong and weak symmetries in triggering critical behaviors, providing a Liouvillian framework to study the effects of high-order nonlinear processes in driven-dissipative systems, that can be applied to problems in quantum sensing and information processing.

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 dictates the system's dynamics and determines its steady states, whose properties can 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 in 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 have been the focus of a vast 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 characterizing DPTs [10].First-order DPTs are discontinuous changes in the properties of the system's steady state as a function of a control parameter [23].These have been associated with hysteresis and critical slowing down [11,24,25], which allow one to observe the emergence of metastable dynamics and to study its competition with the other typical timescales of the system.Key to understanding second-order DPTswhere the steady state transitions continuously, but it is characterized by a divergent response functionare symmetries [11,26,27].In particular, DPTs can be associated with weak and strong symmetry breaking [28,29].Second-order DPTs are useful 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 parametricallydriven (or two-photon) Kerr resonator has attracted considerable interest [10,11,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 DPTs in the coherently-driven (one-photon) Kerr resonator [11,24,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 singleresonator level provided a guideline to investigate emergent phenomena in more complex lattice architectures [27,[47][48][49][50].The fundamental theoretical interest in studying the properties of single or few resonators [10,24,26,[51][52][53][54], but with n-photon driving schemes, is also strongly motivated by the recent achievement of higher-order photon pumping exploiting strong nonlinearities in superconducting circuits [55][56][57].

Summary of the main results
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.We provide a general criterion for arbitrary n to argue the nature of criticalities for this ubiquitous class of models in quantum optics.To this aim, we prove a no-go theorem in the semiclassical limit, predicting the emergence and the nature of DPTs in this class of systems.We derive a general recipe to connect the order of the DPTs and the presence of certain Hamiltonian terms, highlighting the importance that strong and weak symmetry has in constraining the system dynamics, and we point out the technical limitations in witnessing second-order DPTs for n > 4 driving schemes.
We then test the semiclassical predictions by performing a detailed numerical analysis of the full quantum model for the n = 3 and n = 4 cases.We confirm the role of symmetries predicted by the semiclassical theory, and we analyze the DPTs within the theoretical framework of the spectral properties of the Liouvillian superoperator.For n = 3, we confirm the semi-classical prediction and show that the system can only undergo a first-order phase transition accompanied by the symmetry breaking of the discrete weak Z 3 symmetry, as the system parameters are scaled towards the thermodynamic limit.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 , showing that 9 states can be stabilized across a firstorder DPT, making these systems possible candidates for associative memories [21].
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.

The model
We consider a bosonic n-driven nonlinear resonator whose Hamiltonian reads where â (â † ) 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 â † â, U 2 is a standard Kerr interaction, and so on.As detailed also in Appendix A, for a n-photon drive we should consider at least processes up to the order m max = ⌊n/2 + 1⌋, where ⌊A⌋ 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 (Gorini-Kossakowski-Sudarshan) Lindblad master equation [58,59] reading (hereafter we set ℏ = 1) 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 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 ̸ = 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.

Liouvillian spectrum, symmetries, and their breaking
Our analysis will mainly focus on the steady states ρ(k) ss , i.e., that density matrices that do not evolve anymore under the action of the Lindblad master equation (2), defined by k is an index that 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 ρss .DPTs occur when the steady state of an open quantum system can display a nonanalytical behavior as a function of a generic parameter ζ [10,23].As DPTs cannot occur in a finite-size system, one needs to investigate the so-called thermodynamic limit, formally defined as L → ∞.While in a lattice system, one can think at L as the number of sites, in the case under consideration scaling L towards the thermodynamic limit implies a rescaling of the system parameters, as detailed in Sec.3.1.The non-analytical change of the steady-state is then witnessed by the expectation value of some operator ô as ζ crosses the critical point ζ c .As such, we say that there is a phase transition of order M if [10] The n-photon driven Kerr resonator explicitly displays a Z n symmetry1 .That is, the transformation â → â e i2πk/n , k = 0, 1, . . ., 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 Ẑn = e i2πâ † â/n acts.One speaks of strong symmetries if Ẑn commutes with both the Hamiltonian and the jump operators, i.e.: In this case, Z n implies the existence of a corresponding conserved quantity ⟨ Ẑn ⟩ t ≡ Tr[ρ(t) Ẑn ] = const.The system will display n independent steady states, each one characterized by a different value of ⟨ Ẑn ⟩ ss ≡ lim t→∞ ⟨ Ẑn ⟩ t .In our case, such a condition is fulfilled if and only if γ = 0 (i.e., the photons are never lost individually).The presence of a strong symmetry implies that there exist two superoperators 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 ⟨ Ẑn ⟩ t changes in time [8,9].However, the superoperator Z n = Ẑn • Ẑn commutes with the Liouvillian, i.e. [L, As a consequence of the conditions in Eqs. ( 7) and ( 8), strong and weak symmetries constrain the structure of the Liouvillian 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 L, we can introduce its eigenvalues λ i and right eigenoperators ρi , defined via the relation where Re [λ i ] ≤ 0, ∀i represents the decay rates induced by the dissipative dynamics [4,61].

DPTs and weak symmetries
The presence of a weak symmetry allows for refining the discussion on the spectral properties of the system.The eigenvalues z Since each eigenstate of L must also be eigenstate of Z n , we can introduce the "quantum number" k, such that, for a weak Z n symmetry, We sort the eigenvalues in such a way that 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 For this reason, the eigenvalues λ j (k) and eigenoperators ρ(k) j describe the whole physics within each of the Liouvillian symmetry sectors.
Weak symmetries fix the structure of the eigenoperators, which, on the number (Fock) basis, read (12) where mod(p − q, n) indicates the modulo operation.In other words, ρ(k) j must be an operator containing only elements such that (m − n) is either k, or k ± n, or k ± 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 characterized 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, ρ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 ρss ∝ ρ(0) 0 .Furthermore, the discontinuous behavior of the steady state in Eq. ( 4) is signaled 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 λ (k) m becomes exactly zero, both in its real and imaginary parts, as a function of the parameter ζ.In finite-size systems, phase transitions cannot be observed, and λ 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 ρss ∝ ρ(0) 0 and the eigenoperator ρ(0) 1 display level touching (detail can be found in Ref. [10]).More specifically, λ (0) 1 = 0 at the critical point, and the minimum in λ (0) 1 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 does not evolve anymore under the action of the Liouvillian.In this case, the phase transition is associated with λ becoming and remaining zero in a whole region where the symmetry is broken.For instance, in the case of a Z 2 breaking, λ (1) 0 = 0 after the transition, while for Z 3 one has λ The corresponding states ρ(k) 0 , belonging to different symmetry sectors with respect to ρ(0) 0 , allow constructing the symmetry-breaking steady states.Indeed, by choosing the correct superposition of the form ρj

DPTs and strong symmetries
In the case of a strong symmetry, any eigenoperator is characterized by two quantum numbers 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 L = k R , and the coherences, for which k L ̸ = k R .Consequently, the symmetry sectors are For each of the population sectors, there must exist a well-defined steady state ρ(k) ss ∝ ρ(k,k) 0 (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 the spontaneous symmetry breaking accounts for 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 ρ(k,k) becomes zero in the thermodynamic limit.A spontaneous symmetry breaking, instead, implies that the eigenoperators = 0. Spontaneous symmetry breaking thus implies that quantum superpositions between the states composing , 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].

Semiclassical analysis of the nphoton driven resonator
The equation of motion for the expectation value of the observable â evolving under Eq.( 2) is Due to the presence of non-quadratic terms, these equations of motion cannot be closed, leading to a hierarchy of coupled equations.

The thermodynamic limit and finitecomponent phase transitions
We now introduce the dimensionless parameter L such that ) and we will consider the thermodynamic limit L → ∞.In such a limit (G n ) α U m and (G n ) β η n are constants [for α and β such that (−n/2 + 1)α − m + 1 = 0 and (−n/2 + 1)β − 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. [52] for the n = 1 case.The semiclassical (coherent state) approximation amounts to assuming that the state of the resonator is coherent, i.e., where â |α(t)⟩ = α(t) |α(t)⟩.Accordingly, the equation of motion for the rescaled coherent field α(t) = ⟨â⟩ / √ L leads to a generalized driven-dissipative Gross-Pitaevskii-like equation Equation Eq. ( 18) is independent of L and the photon number scales as N = |α| 2 ∝ L confirming that L → ∞ corresponds to a well defined thermodynamic limit with an infinite number of photons.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 [11,24,36,51,52,54,62,63]. In general, we expect the semiclassical approximation (17) to be valid and predictive in the L → ∞ limit, and far from the critical points where nonlinear processes inducing quantum fluctuations cannot be neglected.This assumption, corroborated by the previously-cited extensive literature corpus, bears resemblance to the mean-field approximation in all-to-all connected twolevel systems, where a similar approximation becomes valid in the limit of an infinitely large number of systems [64,65].

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 , η n , G n }.
Despite the simplification introduced by the semiclassical approximation, Eq. ( 18) cannot be yet analytically solved.At the steady state, i.e., ∂ t α = 0, Eq. ( 18) reads (19) In general, Eq. ( 19) gives rise to multiple solutions for the photonic field α.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 α (obtained by considering m 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 highdensity state) accompanied by a hysteresis region with multiple stable solutions in the thermodynamic limit [11,24,52].
If n ≥ 2 Eq. ( 19) always admits the solution α = α vac = 0.However, the other solutions α 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 (20) where we selected the positive branch of the square root since, up to a phase, one can always choose G n ∈ R + 3 .

Second-order phase transitions and behavior around α = 0
For the class of systems under consideration, a second-order phase transition occurs when the state changes from N = ⟨â † â⟩ = 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 We note that the limit N → 0 + must be taken since Eq. ( 20) is defined only for N ̸ = 0.At this specific value of G n the system is thus allowed to pass from the semiclassical solution α vac = 0 to another stable solution with α ̸ = 0. Notice that Eq. (20) admits at most three possible behaviors around N = 0 as sketched in Fig. 2: • The curve G n intersects 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.
3 This amounts to a change in the initial condition by sending â → âe iφ 0 • The curve G n intersects the zero with a negative derivative (blue line in Fig. 2).This resembles the S-like shape of bistability in one-photondriven 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 Notice that, in the case of a vertical-tangent point, higher-order derivatives need to be computed.

Universal features and a semiclassical no-go theorem
From the remarks in the previous section and using Eqs.(22a) and (22b), 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.
Semiclassical no-go theorem.Consider a nphoton driven-dissipative resonator, with nonvanishing Kerr nonlinearity, governed by the Lindbladian (2).Then, according to the semiclassical equations of motion: (a) a second-order DPT never occurs for odd n; 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 ̸ = 0 or γ ̸ = 0 is given by The case where U 1 , γ = 0, the expansion leads to To prove (a), we consider odd-n, and from Eq. ( 23) we get and therefore the condition (22a) for the occurrence of a second-order DPT is never satisfied.If U 1 , γ = 0, instead, Eq. (24) gives We have therefore proven the statement (a).
Let us now consider the case of even n.From Eq. ( 23) we find that for U 1 ̸ = 0 or γ ̸ = 0 a second order DPT is possible only for n = 2, with a critical point given by Higher n results in 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 = γ = 0. From Eq. ( 24) follows that, for n = 4, a second-order DPT can can place also for Therefore, condition (i) in Eq. ( 22a) is satisfied in the case U 1 = γ = 0.As for condition (ii) in Eq. (22b), we have which can be satisfied choosing U 2 and U 3 with the same sign.Finally, one can easily show that for n > 4, Eq. ( 24) gives G n (N ) = ∞, demonstrating the impossibility of a DPT, thus proving (e).■ Before dealing with the analysis of the full quantum results, let us remark that γ = 0 is impossible to achieve in actual realizations.For many practical purposes, however, one can consider system "sizes" L where, to a reasonable approximation, the role of γ can be neglected, and thus the approximation γ = 0 faithfully recovers the results of finite-time experiments.Furthermore, the detuning terms can be easily manipulated, therefore making it possible to approximately fulfill the condition γ = U 1 = 0 necessary to witness the second-order DPT for the n = 4 case.
We also notice that the mechanism enabling secondorder 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. (28).Finally, we stress that, although second-order DPTs could, in principle, also emerge for even n > 4, these would require setting  U 2 = 0, which, contrarily to detuning U 1 , cannot be easily manipulated.
To demonstrate this fact, we compute ∆n, defined as the minimal possible jump between the vacuum and a non-zero stable semiclassical solution, as a function of both U 2 and U 3 .To find that state, one spans all values of G n and finds the stable state with the minimal number of photons.The results, reported in Fig. 3 for n = 6, show that for all nonzero U 2 , the system always displays a finite jump associated with a first-order DPT, confirming that the system has no second-order transition.

Stability of the vacuum across a first-order DPT and convergence radius of the semiclassical solution
In this section, we show that, within the semiclassical picture, the solution α vac = 0 is always asymptotically stable for n > 2 if γ ̸ = 0 (i.e., in the presence of a first-order DPT).
Consider α = α vac + δα, where δα ∈ C is a small perturbation (|δα| ≪ 1) around the vacuum solution.Plugging the above parametrization into Eq.( 18), and expanding it at the first order in δα, we get where δ⃗ α = (Re[δα], Im[δα]) ⊺ and is the so-called stability matrix.The solutions of Eq. ( 31) are given by δ⃗ α(t) = exp(−λ ± t)δ⃗ α(0), where 1 are the eigenvalues of M. Thus, it is straightforward to conclude that for n > 2 the vacuum solution α vac = 0 is always stable at a semiclassical level for finite single-photon losses since For n = 2 the vacuum gets unstable when which implies Re[λ + ] > 0. Equation ( 34) has significant consequences since it implies that, contrary to the n = 1, 2 case, the semiclassical dynamics never triggers a transition from the vacuum to high-density solutions if γ ̸ = 0.However, as we will see 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 γ = 0, the vacuum is marginally stable, and higherorder perturbation theory is needed to assess the stability of the vacuum.
This stability analysis evaluates the stability of the vacuum solution to a weak perturbation (within the linear approximation) but provides no information about the effect of non-infinitesimal perturbations of the vacuum.To this extent, one can perform an explicit numerical study of Eq. ( 18) considering a set of initial states α(t = 0) = r e iθ (with r ∈ R + and θ ∈ [0, 2π]), where we scan both r and θ.One then studies the evolution of α(t), which can either converge back to the vacuum or reach one of the other stable solutions.One then defines the convergence radius r max as the maximal value of the radius r for which the initial solution converges to the vacuum for all the possible values of the phase θ (see Sec.4.2 for an explicit example).

Multistability of solutions with different number of photons
The solution around N = 0 predicts either a firstor 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 solutions.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 numbers.
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 numbers [c.f.Fig. 4 (b)].This translates in the presence of multiple local minima (or maxima) of the function G n (N ), as shown in Fig. 4 (a).Therefore, one can study the function 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 favor the solution with more or fewer photons.Since the relevance of each interaction term can change in different occupation regimes, several solutions can emerge.The stability of the semiclassical solutions in the presence of quantum fluctuations needs to be numerically assessed.

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.Throughout 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 [66], in conjunction with the algorithm detailed in Appendix B.
Here, we consider the three-photon-driven Kerr resonator governed by the master equation with We focus on the γ ̸ = 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.

Semiclassical vs quantum solution
First, we analyze the photon number as a function of the driving strength G 3 .In Fig. 5, 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 ρss ≃ |0⟩ ⟨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 for L = 15 and the parameters in Fig. 5.The solid black line is the result of full quantum simulations.The semiclassical predictions 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).The results of the eigendecomposition are plotted with solid blue and dashed red lines.coherent states, i.e., where |α 1,2,3 ⟩ are coherent states with the same number of photons and a relative phase difference of ±2π/3, i.e., 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.

Analysis of the first-order transition
To confirm the presence of a first-order phase transition, we plot in Fig. 5(b) the Liouvillian eigenvalue λ (0) 1 associated with the slowest relaxation rate in the steady-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 ρ(0) 1 associated with such a state.According to Ref. [10], in the critical region one can use the eigendecomposition of ρ(0) where ρ± 1 represent the density matrices of the metastable states.As such, we expect that, in the thermodynamic limit, ρ± 1 recover the two stable solutions of the semiclassical theory.We show the eigendecomposition in Fig. 6.We indeed find that the semiclassical approximation qualitatively recovers the results of the eigendecomposition.As discussed in Sec.3.2.3, the semiclassical analysis predicts the presence of a stable vacuum in the whole symmetry-broken region.This is also shown in Fig. 7, where we plot the convergence radius r max of the semiclassical approximation.r max remains finite for all values of the drive, despite it decreasing at larger drives.Nonetheless, this decreases as a function of the pump amplitude G 3 .Within this picture, it is the presence of rare and collective quantum fluctuations that trigger the jump between the otherwise semiclassically stable solutions (faithfully representing the quantum state at both sides of the first-order transition).This phenomenology is typical of first-order DPTs (see, e.g., Refs.[25,47]).
This analysis is confirmed both by the eigendecomposition in Fig. 6 (the vacuum remains long-lived even far from the transition point) and from the spectral analysis in Fig. 5(b).Considering larger values of L results in slower timescales.We confirm the scaling towards the thermodynamic limit of the Liouvillian gap λ (0) 1 in Fig. 8.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. 5(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.

Spontaneous symmetry breaking
The spontaneous symmetry breaking implies that, for strong enough pumping, each of the state |α i ⟩ ⟨α i | in Eq. ( 38) becomes a steady state of the system, since they are not eigenstates of Z 3 [10].We confirm this  picture in Fig. 5(c), where we show that also λ (1) 0 becomes zero.As expected, we obtain an identical result for λ (2) 0 (not shown).This implies that the system's coherent state has become metastable.

Four-photon Kerr resonator
Here, we consider the four-photon drivendissipative Kerr resonator, reading with (42)

Strong symmetry and second-order phase transition
We start by considering the strong symmetric case γ = 0 with U 1 = 0.For this set of parameters, 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. 9.We recall that since the system has a strong Z 4 symmetry, the number of Liouvillian sectors is 4 × 4, being characterized by the two quantum numbers k L and k R .The 4 sectors with k L = k R describe the evolution of the populations, while the remaining 12 with k L ̸ = k 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 ∈ [0, 3] and increase the thermodynamic parameter L. Calling ρ(j) ss ∝ ρ(j,j,) 0 the steady state in each symmetry sector, ⟨â † â⟩ j = Tr â † âρ j ss are plotted in Figs.9(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 where j labels the symmetry sector and |j⟩ is the Fock state with j photons.For large drive, instead, the system transition towards where the Schrödinger cats |K i ⟩ are where |α n ⟩ are coherent states such that α n = e iπn/2 α and 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 corroborates 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 λ (j,j) 1 in Figs.9(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 λ (j,j+1) 0 for the sectors (j, j + 1) (where j + 1 = 0 if j = 3) in Figs.9(i-l).These represent the decay rate of coherences between the sector j and j +1, and their closure indicates the possibility of retaining everlasting coherences.In this case, we observe that, after the critical point, these eigenvalues 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 ̸ = k (not shown).This is associated with a spontaneous breaking of the strong Z 4 symmetry, because it results in

Weak symmetry and multistability
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. 10(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 solutions, the full quantum simulation is characterized by a single first-order DPT, from the vacuum to the highestpopulated manifold.Indeed, if we analyze the Liouvillian gap λ Liouvillian gap associated with a first-order DPT.If, however, we also consider the second eigenvalue λ (0) 2 as in Fig. 10(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 are characterized by two slow timescales.
We corroborate this phenomenon by analyzing the symmetry sectors responsible for spontaneous symmetry breaking.In Fig. 10(c), we plot λ (1) 0 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 occurrence of the first-order transition.Furthermore, we also observe a second slow timescale for this symmetry sector, i.e., λ (1) 1 in Fig. 10(d).These slow timescales represent the fact that there exist multiple symmetrybroken states, and there is a slow rate at 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 states of the dynamics, either the vacuum or the one at large photon number, there exists 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. 11(a) we use the eigendecomposition to express the eigenoperators associated with the slowest eigenvalues as 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 Liouvil-

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, determine 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 transi-tions 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 [67,68] and Gaussian expansions [49].This is not the case for n > 2 because the drive and n-photon dissipation can, in principle, introduce non-Gaussian correlation above the coherent-state solution [53].The emergence of these 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 the 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 has been recently proposed [28,30], detuning and critical phenomena may play a key role in storing quantum information.The metastability of the vacuum may also prove an obstacle to a rapid and reliable initialization of bosonic qubits.This work paves the way for 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.

A Interaction, nonlinearities, and nphoton drives in superconducting circuits
Let us consider a standard LC resonator characterized by the Hamiltonian: where ω 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 (49) where ϕ is the flux coordinate of the circuit at the junction and ϕ 0 is the magnetic flux quantum.In most implementations, the expansion in Eq. ( 49) can be stopped at the ϕ 4 order.If the Josephson junction belongs to a single resonator, by substituting ϕ/ϕ 0 ∝ â + â † , and discarding counter-rotating terms, which are out of resonance, one obtains the Kerr resonator Hamiltonian, reading 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 â + â † ) can emerge by, e.g., capacitive coupling an incoming wave-guide with the cavity, higher order drive requires 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 ϕ = k ϕ k , where ϕ k represents the flux coordinate of each one of the modes.For instance, two-photon drives can be achieved by standard four-wave mixing, rewriting ϕ = ϕ a + ϕ b + ϕ c , where ϕ a ∝ â + â † is the field within the resonator, and ϕ b ∝ b + b † and ϕ c ∝ ĉ + ĉ † 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 b (ĉ) with a c-number oscillating at the driving frequency ω b (ω c ) via an adiabatic elimination, reading b → be iω b t (ĉ → ce iωct ).All in all, discarding again out-of-resonance terms, the Hamiltonian for the resonator resulting from the fourth-order expansion of the potential cos(ϕ) would result in a nonlinear Hamiltonian ĤNL , reading ĤNL = ĤKerr + G 2 â2 e 2iωpt + â † 2 e −2iωpt .(51)

Figure 2 :
Figure 2: (a) Possible different behaviors of Gn as a function of N according to Eq. (20).The marker indicates the minima of the function Gn, 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 Gn.

Figure 3 :
Figure 3: In a system with n = 6 photon drive, analysis of the minimum possible jump in the photon number ∆n as a function of U2 and U3.This is defined in the main text as the stable semiclassical state with the smallest nonzero amplitude for all values of G6.Parameters: γ = 0, U1 = 0, and η4/U4 = 0.1.

Figure 5 :
Figure 5: 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 ⟨â † â⟩, renormalized by scaling parameter L. Panel (b): Real part of λ (0) 1 , i.e., the Liouvillian gap in the same symmetry sector as the steady state, inducing the first-order transition.The three vertical lines indicate the scaling values studied in Fig. 8. Panel (c): Real part of λ (1) 0 , i.e., the Liovillian eigenvalue signaling the spontaneous symmetry breaking.Parameters: U1/γ = −20, U2/γ = 10, η3/γ = 1.The cutoff Nc, also used in Figs. 6 and 8, is chosen to ensure that, between the simulation with Nc and Nc + 10, all quantities differ by less than 1%.The cutoffs are: Nc = 30 for L = 1; Nc = 95 for L = 5; Nc = 140 for for L = 10; Nc = 190 for L = 15.

Figure 7 :
Figure 7: a function of the drive amplitude, the convergence radius rmax, numerically computed according to the procedure in Sec.3.2.3.The dashed vertical line indicates where the high-photon solution becomes stable.Below this value, the convergence radius becomes infinite.Parameters as in Fig. 5.

Figure 8 :
Figure 8: Scaling towards the thermodynamic limit for the three vertical lines in Fig. 5(b), demonstrating the presence of a first-order DPT and the emergent stability of the vacuum.

Figure 9 :
Figure 9: Study of the strongly-symmetric four-photon Kerr resonator, and 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: γ = U1 = 0, U2 = 10η4, U3 = η4.The cutoff Nc is chosen to ensure that, between the simulation with Nc and Nc + 10, all quantities differ by less than 1%.The cutoffs are: Nc = 45 for L = 1; Nc = 75 for L = 2; Nc = 150 for for L = 5; Nc = 290 for L = 10.

Figure 10 :
Figure 10: 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 smallest Liouvillian eigenvalues in the k = 1 sector, demonstrating the presence of SSB and of a slow timescale.Parameters: U1 = 10γ, U2/γ = −25γ, U3 = 3γ, η4 = 0.1γ.The cutoff Nc, also used in Fig. 11, is chosen to ensure that, between the simulation with Nc and Nc + 10, all quantities differ by less than 1%.The cutoffs are: Nc = 40 for L = 1; Nc = 60 for L = 2; Nc = 100 for L = 5; Nc = 190 for L = 10; Nc = 280 for L = 15; Nc = 320 for L = 20.

Figure 11 :
Figure 11: 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. 10 for L = 20.