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

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


I. INTRODUCTION
Understanding the quantum phase transitions (QPT) of many-body systems and their time evolution are among the main goals of modern physics.Historically, phase transitions have been introduced and investigated in the thermodynamics limit of closed systems where they thermalize in the steady-state.In recent years, however, phase transitions and critical phenomena in drivendissipative many-body quantum systems have emerged as a major field of research.Extensive studies have been exploring the new classes of quantum phase transitions in the non-equilibrium steady-state (NESS) of open quantum systems, on account of experimental realization of dissipative quantum simulators [1][2][3].
Experimental platforms, such as cavity arrays, superconducting circuits, and exciton-polaritons, have put forward versatile testbeds to examine the interplay between (in)coherent drive, dissipation, and interaction on NESS phases.This includes the multi-stability and crystallization in driven-dissipative nonlinear resonator arrays [4][5][6][7] and spins [8], and synchronized switching in the arrays of coupled Josephson junctions [9], to name a few.Limitcycle (LC), aka synchronization, is an intriguing phase of open quantum systems, where the dynamics traverses a closed trajectory in the phase space [10][11][12].
In this article, we ask the question of what are the elements that appear typically in the phase diagram of a dissipative non-equilibrium system.The focus is on directly experimentally realizable systems that allow in complexity at least for limit cycles.The theoretical methodology will be built upon typical tools, such as mean-field theory, and expanded as far as necessary to show and confirm the complex results.The story we tell in this article, thus, introduces first the physical system, which consists of a dissipative multi-mode cavity shown in Fig. 1 which is driven on one mode and connected to two side modes via a Kerr nonlinearity.Depending on the parameters, a complex phase diagram with various steady states, bistability, or limit cycles ensues.As a benchmark, this system is compared to the same system without the dissipation and a much simpler phase diagram that consists only of two types of steady-state.This is a complex interacting system where meanfield theory is not automatically expected to provide the correct results.We therefore support and correct the results using two different methods, (i) employing a thermodynamic-limit Keldysh formalism and (ii) including higher-order correlations.Both these methods confirm the phase diagram and provide the quantitative corrections.The results and the methods presented in this work enable one to study the quantum properties of NESS beyond the semi-classical approximation and shed light on the emergence of dissipative phase transitions in open quantum systems.
The paper is structured as follows.In section II we introduce our model of a driven three-mode cavity with self and mutual Hubbard-type interaction subject to singlephoton (coherent) drive.When the cavity-environment can be described by a Markovian process, the joint density matrix of the cavity evolves according to a Lindblad master equation.In section III we describe the physics of the closed system which we later use as a benchmark to highlight the profound effects of dissipation.Section IV presents the numerical results of an attractive interaction and details the MF phase diagram, explains the uniform and multi-stable regions for the pumped mode, and the emergence of the LC-phase for un-pumped side modes.Moreover, we elaborate on the spectral features of various phases obtained from Keldysh Green's functions, and showcase dissipation-stabilized cross-overs in the open system.The beyond-MF results employing the 2 nd -cumulant approach is presented in sub-section IV C where we compare the Gaussian approximation of the cavity states with their corresponding MF results.The implications of quantum correlation on LC is detailed in this sub-section.In section V we briefly summarize our Keldysh formalism and the relevant Green's function presented in this work.Finally, section VII concludes our results and presents an outlook for follow-up works.

II. THE MODEL
We consider the dynamics of a three-mode lossy cavity with self-and cross-Kerr non-linearity, i.e. a Hubbard interaction, subject to a single-photon coherent drive and a single-photon loss [23].In the laser rotated frame, the conserved dynamics of such system is given via the fol-lowing Hamiltonian ( = 1) where the 2 nd -mode is subjected to a single-photon coherent drive at the rate of Ω 2 and the frequency of ω L .âm , â † m is the annihilation and creation operators of the m th -cavity mode, respectively, ∆ m = ω m − ω L is the detuning of the m th cavity mode from the coherent drive, and U 0 is the interaction rate.
For U 0 ≥ 0 the system energy increases with increasing the particle number hence, a repulsive interaction.Similarly, when U 0 ≤ 0 the system energy decreases with increasing the particle number and the interaction is attractive.
In the presence of an incoherent single-photon loss with the rate of 2γ m from the cavity, the Markov-Born approximation leads to the following Lindblad dissipator for the cavity density operator The time evolution of the multi-mode cavity density operator ρ(t) is determined via the following master equation In this work, we are interested in the long-time solution of this density matrix, where all the transient dynamics are over.This is obtained by numerically integrating the coupled equations of motions (EoM) [cf.Appendix A] implemented by the 4 th -order Runge-Kutta method.
When the driving laser is swept from the red to the blue detuning, the nonlinear dynamics sets in as a parametric amplification process and photons are created in the side modes as correlated pairs.As shown in Appendix A, for an attractive interaction, a multi-stability region exists for the red-detuned driving fields satisfying ∆ 2 ≥ √ 3γ 2 .Similarly, for repulsive interactions a multi-stability phase exists for the blue-detuned coherent drives when ∆ 2 ≤ − √ 3γ 2 .For some pumping rates and detuning, the amplification moves into a regime of self-sustained oscillations, aka limit cycle, due to the nonlinear self-and cross-Kerr coupling.In the next section, we detail the phase diagram of the aforementioned system assuming an attractive interaction.

III. CLOSED SYSTEM
To draw a preliminary understanding of our model we first describe the closed system.We begin by noting that the Hamiltonian of Eq. ( 1) is unchanged under the exchange of the un-pumped modes, a 1 ↔ a 3 .Additionally, there are neither terms that mimic a single-nor a twophoton pump for these modes taken separately.Therefore, a symmetry-breaking due to population of modes a 1 and/or a 3 is only possible if both of them are populated at the same time.On the contrary, due to the single-photon pump of the 2 nd -mode, the system features a phase transition as a function of the pump strength Ω 2 .We study the mean-field energy potential landscape, H3 , for a zero occupation of the 1 st and 3 rd -mode [55].
where α 2 is the order parameter.Without loss of generality, we assume Ω 2 to be real and find the potential extrema to be given by This leaves us with a cubic equation in Re(α 2 ).We study the discriminant and find the boundaries between the regions with only one and three possible solutions where we explicitly assumed U 0 < 0 in our case.
In the parameter regime where three solutions are allowed, we find one low (LP) and one high population (HP) phase accompanied by an un-physical one, i.e. imaginary eigen-frequencies of excitations.
To understand the nature of the extrema of the potential, we find the excitation spectrum associated with the uniform case Hamiltonian (4) [3].We obtain the eigenfrequencies with eigen-vectors and their associated symplectic norms where The symplectic norm describes the nature of a state of the system.Whenever all positive(negative) frequency eigenmodes have a positive (negative) symplectic norm the state of the is the ground state.On the other hand if a physical state of the system, i.e., a state with real excitation eigenmodes, has at least one excitation mode with negative (positive) norm for positive (negative) excitation frequency, then the state is an excited state of the closed system.Additionally, a positive symplectic norm is associated with particle-like processes where an external excitation is absorbed by the system whereas a negative norm underpins a hole-like process where an excitation in the system is destroyed [3].Finally, the symplectic norm helps identifying the so called negative-mass instabilities [30].
In Fig. 2, we plot the three-mode mean-field energy potential of Eq. ( 4).Due to the single-photon pump of the 2 nd -mode, the system features a phase transition as a function of the pump strength Ω 2 between an LP and an HP phase.The phase diagram comprises of two qualitatively-distinct regions in parameter space with: I, the energy functional has only one clear extremum, the HP phase; II, three extrema exists including a "proper" LP phase ground-state, a saddle point, and a HP phase that reveals itself as a maximum, i.e. an excited-state.Two distinct regions are indicated by their respective meanfield energy potential ( H3) landscape as a function of the real part of the cavity field, Re(α2).The dark-blue region indicates the parameter regime where the LP phase is the ground state of the system and the HP represents a physically-allowed excited state.The LP phase disappears in the light-blue region where the HP is the only physical state.
In Fig. 3, we plot the order parameter |α 2 | 2 and the excitation spectra on top of the possible solutions along the vertical orange cut (a,c) and the horizontal brown cut (b,d) of Fig. 2. In both cases the order parameter |α 2 | 2 has a finite value within the entire range but in region II, two distinct phases with different populations are possible.Besides, there is a third phase whose spectrum is fully imaginary and therefore not a physical state (gray lines in Fig. 3(a),(b)).The HP exists throughout the entire phase diagram but always presents a negative (positive) symplectic norm with a positive (negative) eigenfrequency.Therefore, we confirm that it is indeed an excited-state of the closed system.Due to the presence of the negative interactions a "true ground-state" of the unbounded Hamiltonian is only possible in a limited region of the parameter space (see Fig. 2) and is identified with the LP phase, accompanied by a positive (negative) symplectic norm with a positive (negative) eigenfrequency [3,31].

IV. OPEN SYSTEM
For the open system, we determine the phase diagram by initializing the system in many randomized initial conditions and examining the dynamical stability of the end points using the Bogoliubov matrix spectrum [cf.Appendix A].In this respect, the phase boundaries can be thought of as dynamical phase transitions separating distinct long-time asymptotic behaviors [32].
In addition to multi-stable phases, where the final state depends on the initial conditions, it is also possible to find regions of parameter space where no stable fixed point exists.In such cases the system may be attracted to some time-dependent solutions such as limit cycles [56], as found in other coupled nonlinear systems.In particular, we search for the complete set of stable attractors of the long-time dynamics, including the fixed points, the multistable coexistence phases, and time-dependent trajectories.The possible steady-states for the evolution of the system under Eq.( 3) are phases where either the only populated mode is the pumped mode, (α n=2 = 0, α n =2 = 0), called the uniform phase, or phases where the side modes, a 1 , a 3 , are also populated and host a limit-cycle.
In Fig. 4 we plot the three-mode cavity open phase di-agram.It comprises four different regions: (I) the white region that presents only one uniform stable steady-state which we identify as the LP phase, (II) the dark blue region of tri-stability between two uniform phases namely, the LP and the HP one, and a non-uniform solution where the un-pumped modes showcase a limit-cycle behavior, (III) the blue region of bi-stability between the LP and the HP uniform phases, and finally, (IV) the light blue region that presents a HP uniform phase.At the boundary between region (I) and (IV), denoted as the striped region, an exceptional-point appears and the LP and HP are smoothly connected.
When compared to the phase diagram of the closed system [cf.Fig. 2], it becomes clear that the dissipation has a marked impact.It lifts the boundary of the closed phase diagram (red dashed line), stabilizes an excited state of the closed system, the HP, in region (IV), renders the LP as the only attractor of the dynamics in region (I), and leads to the emergence of limit-cycles.

A. Steady-state
In this section, we focus on different points of the open phase diagram in order to clearly expose the nature of the underlying attractors including their stability.As discussed in Appendix A, the possible fixed points of the EoM are determined as d dt âm = 0 in Eq.A2.Among the possible fixed points in the long-time limit, the semiclassical dynamics eventually evolves towards some stable fixed points for all initial conditions, known as the steady-states of the system.
We use the mean-field approximation, i.e. we factorize higher order moments as â † m ân ≈ â † m ân , and obtain the equations of motion for the MF order parameters α m , with m = 1, 2, 3.When U 0 ∆ 2 ≥ 0 the EoM has one stable solution only as (α n=2 = 0, α n =2 = 0), whereas when U 0 ∆ 2 ≤ 0, there might exist several stable steady-state solutions [cf.Appendix A].We self-consistently determine the steady-state of the system which, in general, may allow limit-cycle solutions.For these solutions, the long-time limit of the steady-state has the general form of e itωLC .Physically, the LC solutions and their associated frequencies (ω LC ) can be thought as the frequency of the parametrically generated pair, i.e., the un-pumped modes, via the parametric Kerr process.
To draw a more comprehensive picture of the steadystate behavior of the system, in Figure 5(a)-(c) we show the mean-field occupation of each cavity mode, n m = |α m | 2 , versus the pumping rate Ω 2 at various detunings, ∆ 2 = -5, 0, +5 (vertical dashed lines in Fig. 4), for an attractive interaction of U 0 = −1, and bare cavity spacing as ∆ 1,3 = ∆ 2 ∓ 1.The red lines show the behavior of the 2 nd -mode while the blue lines correspond to the 1 st -and the 3 rd -mode populations.The solid and dashed lines signify the stable and unstable solutions, respectively.For U 0 ∆ 2 ≥ 0, panels (a)-(b), there is only one stable steady-state, i.e. the uniform phase HP, but when U 0 ∆ 2 < √ 3 we observe the regions of multi-stability and transitions to the non-uniform phases, Fig. 5(c).
In Fig. 5(d)-(f) we show three exemplary MF time traces for randomized initial conditions corresponding to the points T 1,2,3 , i.e.Ω 2 = 3, in Fig. 4.They are obtained from the direct integration of the EoM and illustrate the implications of a limit-cycle phase and its difference with the uniform one.For each detuning, the upper row shows the temporal behavior of n 2 = |α 2 | 2 while the lower panel shows the real part of the 1 st -mode order parameter, i.e.Re(α 1 ).
For T 1 and T 2 , corresponding to ∆ 2 = −5, 0 respectively, there is only a uniform phase where α 1,3 = 0 and α 2 = 0 [cf.Appendix A].Accordingly, in Fig. 5(d),(e) we see that all time traces converge to only one non-zero value (HP) for the pumped mode and a zero value for the side modes, independent of the initial conditions.
On the other hand, for T 3 corresponding to ∆ 2 = +5, where a multi-stability is predicted [cf.Appendix A], one can see that the time traces in the upper panel of Fig. 5 and (e) the system always evolves towards a uniform steadystate with only the pumped mode populated (HP).In (f) there are three possible steady-states, two uniform phases LP, HP (red and brown) and a limit-cycle (LC, orange).In the limit-cycle case the side-modes show an oscillatory behavior at the frequency ωLC, as depicted with light blue traces in the bottom panel (f).All modes have the same decay rate γ0 to which other rate parameters are normalized.The time is in units of 1/γ0.U0 = −1, Ω2 = 3 in all temporal calculations.The side modes are equally-spaced around the pumped mode as ω1,3 = ω2 ∓ 1.
From those three phases only the LC corresponds to a non-zero values for the 1 st , 3 rd modes (orange line in upper panel of the Fig. 5(f)), where their order parameter Re(α 1 ) shows a periodic long-time behavior (ligh blue in the lower panel of Fig. 5(f)).

B. Uniform phase crossover
To illustrate the dissipation effect in modifying the closed-system phases, in Fig. 6(a,c) and (b,d), we plot the population of the 2 nd mode, |α 2 | 2 , (top) and the fluctuation eigenvalues (bottom) along the vertical (∆ 2 = +5) and horizontal cuts (Ω 2 = 0.5) of Fig. 4, respectively.In both cases, |α 2 | 2 behavior highlights a smooth cross-over from the HP to the LP between region IV and I.As can be inferred from the phase diagram, this transition traverses an exceptional point region where the imaginary part of the eigenvalues coalesce to zero and their real parts split, hosting over-and under-damped fluctuations.Using a Keldysh action approach we readily get access to the Green's function of the system and its associated dynamical observable.In Fig. 7, we compare the spectral functions, A(ω) = −2Im[G R (ω)], of the three-modes in the LP (a) and HP (b) phases for different detuning and at fixed pump strengths Ω 2 = 0.5, points P LP,HP in Fig. 4. From the left panel we see that the LP phase has the response of a "ground-state", i.e., a positive (negative) peak at the positive (negative) frequencies [cf.Sec.III].We note that peaks at negative frequencies are unresolved due to the scale resolution.On the other hand, the HP phase features a peak swap with a positive (negative) peak at the negative (positive) frequencies, a hallmark of a stabilized excited-state [cf.Sec.III and Fig. 10 for more information].We note that the two uniform phases have different occupations for the 2 nd -mode, see Fig. 5 LP and HP.Nevertheless, their spectral functions highlight a more profound difference rather than just a different mode occupation.Interestingly, the peak swap is also present in the response of the empty side modes and we can think of it as being in presence of a normal phase to excited normal phase transition [31].
Even though the LP and HP phases are not described by the same order parameters, through the study of the dynamical responses we identified a peak swap between their spectral functions and we traced it back to the LP being the ground-state of the closed system and the HP an excited-state [3,33].The dissipation stabilizes an excited-state, the HP, and lifts the boundaries of the closed system phase diagram changing its topography, see Sec.III.The spectral functions shown in Fig. 7

C. Gaussian approximation and beyond the MF results
So far, we have studied the MF-phase diagram of the open system and investigated the disspation-stabilized phases in contrast to the closed system.Moreover, using Keldysh approach, we have been able to investigate the spectral signatures of each phase and delimit the onset of a PT, beyond the MF.
As can be inferred from the Langevin EoM in Eq A1 of Appendix A, the quartic interaction leads to an infinite hierarchy of moments.Therefore, any semi-analytic or numerical calculations require a truncation of this hierarchy.The mean-field treatment ignores the correlation via the factorization approximation, hence truncating the cumulant expansions to the 1 st -order.Recent studies however, show that the inclusion of quantum correlations lead to marked deviations from the mean-field results, especially close to the phase transition points [34].So far, there have been several approaches to include the effect of higher-order correlation including exact diagonalization, diagrammatic expansions, functional renormalization group, numerical or density matrix renormalization group analysis, and phase space methods.[35][36][37][38][39][40][41][42][43][44].
In this section, we extend the EoM to include the moment dynamics up to the 2 nd -order, while assuming a vanishing 3 rd -order, to find the Gaussian approximation of the NESS with an emphasis on the robustness of the LC-phase.Further details for this specific problem and the comparison between the results of this approach with the MF and the full density matrix calculations for a single-mode Kerr cavity can be found in the Appendix C.
For panels (a),(b) corresponding to the MF uniform phase, i.e. one solution for the pumped mode and no occupation of the side modes, the Gaussian approximation results are in good agreement with the MF ones depicted in Fig. 5(a),(b).For the side-modes however, the MF predicts zero population while the 2 nd -cumulant suggests a finite but low value.This can be understood in terms of the quantum fluctuations that had been ignored in the MF, and partially resumed in the Gaussian approximation.Besides, the weak correlation between the side modes, i.e. |g 13 |, reinforces that interpretation.
For ∆ 2 = +5 however, the MF (Fig. 5(c)) and the Gaussian approximation (Fig. 8(c)) show marked differences.While the MF shows a multistable behavior, a 1 st -order phase transitions, and an LC behavior for the side modes, the Gaussian approximation results are continuous for all modes.For the side modes and before the LC-phase, the trend is quite similar to the other detuning in panels (a) and (b), i.e. a finite but low population.Unlike the MF case however, the transition to the LC-phase is not a 1 st -order as suggested by the MF and instead it shows a large but continuous change of the order parameter, â † 1 â1 .This transition is accompanied with a maximum in the correlation between the side modes,i.e.â1 â3 , a quantity that signifies the correlated photonpair generation within this phase.It is interesting to note that this correlation is low before the LC-phase and drops but remains finite after the LC.Moreover, the MF predicts a uniform phase again after the LC, where there is no population in the side modes, while the Gaussian approximation delineates that a finite population after the LC is the correlation effect, only.
Figure 8(d)-(f) presents the temporal evolution of the pumped mode for the corresponding detuning in (a)-(c) at a fixed pumping rate of Ω 2 = 3 starting from vacuum, i.e. â † m âm = 0.For the uniform cases in (a),(b) the order parameter settles to a fixed value at the long-time limit after some transient behavior, a value which is very close to the corresponding MF value.
In contrary, the order parameter of the LC phase shows sustainable oscillations in time, after a short transient period, corroborating with the LC-phase features.It is interesting to note that, unlike the MF-predicted LC how-  4. Using 2 nd -cumulant approximation (solid lines) reveals a frequency pulling towards lower frequencies as compared to MF calculation (dashed lines).The effect is more pronounced in the high density phase where the corrections to the MF are bigger.
ever, where the time-periodic behavior was only observable in the side-modes, here the coupled correlations lead to oscillatory behaviors in all correlations.
In Fig. 9, we compare the spectral functions of the modes obtained from the Keldysh approach based on the 2 nd -cumulant results, with the ones determined from the MF-values.There is no qualitative change in the response and the peak swap between the LP and HP is still present.The 2 nd -cumulant approximation reveals corrections to the eigen-frequencies of the system highlighted as a frequency pulling towards lower values.The effect is more pronounced in the HP, in agreement with the presence of stronger corrections to the MF results where the Gaussian approximation predicts a non-vanishing population for the un-pumped modes in contrary to zero MFvalues.

V. KELDYSH FORMALISM
The exact diagonalization of the Liouvillian superoperator suffers from the finite size effects, i.e. the truncation of the Fock space, while the MF approach ignores quantum fluctuations.To include the correlation effects while approaching the thermodynamic limit (the typical validity range of a MF treatment), we employ Keldysh formalism [45].
We readily write the rotated Keldysh action, associated with the Hamiltonian (1) and subjected to the dissipation (2), in the quantum and classical fields a mc,q , for m = 1, 2, 3 as where S m 0 , S int , S drive , S γ are defined in Appendix B. Here, we only highlight that the symmetry a 1 ↔ a 3 is still present in the action (11).
We perform the saddle point approximation via ∂S k /∂α c,q = 0, set the quantum fields to zero α mq = 0, and obtain the EoM for the classical fields α cm [cf.Appendix B].These equations coincide with the mean-field ones upon rotating back to physical fields, i.e., α mc = √ 2α m .
To go beyond mean-field, we study the fluctuations around the MF stationary-states as α m = α mc + δα m .We expand the Keldysh action in Eq. ( 11) and retain terms up to second order in fluctuations, i.e.only the Gaussian parts.Therefore, the Gaussian action can be written in the normal form as where we went to Fourier space with the 3-mode, 12 component, nambu-spinor δΦ(ω), and the matrix M (ω) is given by where Using Gaussian integration we determine the singlemode action and the associated Green's functions.We here report solely the uniform phase case, α 2 = 0, α 1,3 = 0 hence, the single-mode Green's functions has simple analytical expressions where n 2 = |α 2 | 2 and (m, m) = (1, 3), (3,1).
From the poles of the Green's function we obtain the eigenvalues ) where ω i is associated with the i th cavity mode and i = 1, 2, 3. First, we notice that the pumped mode behaves as a single driven Kerr-oscillator whereas the symmetric modes a 1,3 , even though empty and not pumped, present a non-trivial Green's function, and their Green's functions are "coupled".Second, the eigenvalues we obtained from the Green's functions coincide, upon the substitutions , with those obtained from the stability matrix [cf.Appendix A].

VI. EXPERIMENTAL REALIZATIONS
We consider a multi-mode spherical Fabry-Perot cavity, as sketched in Fig. 1.Depending on the mirrors radii of curvature and the cavity length, the cavity leads to an effective harmonic trap for the photons and sets the bare cavity spectrum, i.e. ω 1,2,••• [46].Due to the finite reflectively from the mirrors, there will be a cavity decay rate γ n corresponding to the n th -mode, which can be approximated with a fixed rate γ 0 , over a finite frequency range of our interest.The self-and cross-mode interactions are induced via a non-linear medium inside the cavity, ranging from ultracold atomic gas trapped along the cavity axis [47][48][49], to excitons [50] or highlyexcited Rydberg excitons in solid state.By re-writing the dynamics in terms of the polaritons, i.e. the hybrid particles of cavity modes and the material degrees of free-dom, we can get the final dynamics of Eq. 1.In brief, ân , â † n represent the annihilation and creation operator of a polariton in the n th -cavity mode, and U 0 is the effective interaction between polaritons after approximating the two-body interaction, or the third-order optical non-linearity, with the short-range contact interaction.By tuning the relative frequency between the atoms the cavity modes, the interaction can can be attractive or repulsive.A laser with frequency ω L coherently drives the 2 nd -cavity mode.Note that due to the orthogonality of the Gauss-Lauguerre modes of spherical cavities, excitation of a specific mode is possible.Although here we only focused on a canonical cavity-QED setup, recent developments of the multi-mode superconducting cavities with parametrically-induced self-and cross-Kerr non-linearity allows ones to engineer similar Hamiltonian using the circuit-QED toolbox, as well [51].

VII. CONCLUSION AND OUTLOOK
For a long time it has been speculated if the limit cycle, a recurring phenomenon in many non-linear classical and quantum many-body systems, can be classified as a time crystal.Although the semi-classical MF treatment of several many-body systems suggests the appearance of periodic long-time behavior, however, the robustness of these periodic dynamics, against quantum fluctuations has been an open question.In this work, we investigated this question in a canonical yet important system, i.e. a multi-mode cavity with self and cross-Kerr non-linearity.
We used the Keldysh formalism and extended the MF results via a higher-order cumulant expansion.Both of these methods confirmed the robustness of the LC phase phase and the overall topology of the phase diagram.Additionally, we revealed corrections to the MF results, e.g. the change of the correlation spectra, the non-vanishing population of the side modes, and the oscillating crosscorrelations.
The results and the methods presented in this work can be employed in the study of a versatile group of the manybody systems to understand the effects of higher order correlations and investigate the robustness of dissipationstabilized MF phases, including the limit cycle.It is straightforward to extend this work to study the dynamical behavior of such systems subject to a parametric drive with a modulated amplitude, and investigate the system transition from a stationary dynamics, to a limit cycle, to a chaotic phase as a function of the modulation frequency and depth.
Another interesting direction is to extend the system size either by including more modes in a single cavity or considering an array of coupled cavities with self-and cross-Kerr non-linearity.When subject to a two-photon parametric drive, instead of a single-photon as studied here, the Z 2 symmetry of the modes are preserved and new phase transitions could appear due to the sponta-neous breaking of this symmetry [52,53].It would be interesting to examine if for some parameter range this symmetry breaking can be accompanied with a limitcycle phase.Since such a system, at certain limits, can approximate a spin Ising Hamiltonian, it is interesting to investigate the notion of the LC phase on the magnetization.Finally, one can explore the ultimate behavior of cascaded quantum systems where the limit cycle output of the first sub-system behaves as a parametric drive for the second one.This could be thought as an alternative approach for creating time crystals using discrete time-symmetric drive.
for N = α 2 1 + α 2 2 + α 2 3 being the mean-value of the total number of the photons in the cavity.
When MF is dynamically-stable the Bogoliubov matrix M is negative-definite.For the stationary-state of the stable solutions we can define the covariance matrix as From this equation it is clear that the spectrum of the correlation matrix is closely related to the noise operator spectrum.The diagonal entries of the covariance matrix are related to the cavity transmission spectrum (auto-correlations), and the off-diagonal entries signify the cross-correlations between the modes hence, are related to the entanglement between the modes.

uniform phase
Within this phase, only the 2 nd -mode has a non-zero MF and α 1,3 = 0, making R, S diagonal and antidiagonal matrices, respectively.The eigenvalues of the dynamical stability matrix M for identical loss rates of the cavity modes γ 0 , are where the dispersion parameter δ D is defined as 2δ As can be seen, for U 0 ∆ 2 ≥ 0 and U 0 (∆ 2 − δ D ) ≥ 0, matrix M is negative definite in the uniform phase hence, stagnation points are dynamically stable.
The many-body system is within the uniform phase when Ω 2 is either small or large.If the former, the cross-interaction compared to the self-interaction is so small that the pair generation process cannot start.In the latter, the number of particles in the pumped mode is so large that the self-interaction shifts the pumped mode out of resonance by several γ 0 such that the intermodal scattering ceases.In other words, in the extreme of a strong pump, the large self-interaction dominates all many-body interactions hence, pushing the system to the single-body dynamics again.
On the single-mode branch only, one gets When U 0 ∆ 2 ≤ 0, there might be points that the slope first diverges and later changes the sign.These turning points in the MF dynamics signify the existence of multiple MF attractors.More specifically, the boundaries of the uniform phase, can be determined as The above equation determines that for δ D = 0 the multistability exists when ∆ 2 ≥ γ 2 √ 3.
A similar argument for a dispersive cavity, i.e. δ D = 0, shows that the multi-stability can exist for U 0 (∆ 2 −δ D ) ≤ 0 and The covariance matrix within this uniform phase has the following general form .
The structure of this matrix implies that the 2 nd -mode has no correlation with the side modes.That physically is consistent with the picture of these un-pumped modes not being populated through the parametric process via the 2 nd -mode.Consequently, the covariance matrix for the quadratures gets the following block diagonal form The above form emphasizes again that the pumped mode has no correlation with the other two side modes.

Limit-Cycle phase
As discussed before, while the U(1)-symmetry of the un-driven system is broken by a coherent drive, the unpumped modes still have some phase freedom remained since Φ 0 = 2φ 2 − φ 1 − φ 3 , is the only constraint imposed by the coherent pump.
Within this phase there is no stationary state and the long-time limit of the 1 st , 3 rd MF has an oscillatory behavior as e ±itωLC .In other words, due to the aforementioned phase constraint, if â1 oscillates as e itωLC , â3 should vary as e −itωLC .Going back to the Eq.A1, one can see that this oscillation can be interpreted in terms of a re-normalized detuning of the parametrically-populated modes as ∆1,3 = ∆ 1,3 ± ω LC .

Appendix B: Rotated Keldysh action
In this appendix, we give the explicit expressions for the rotated Keldysh actions used in the Eq. ( 11) of the main text.
The Keldysh action of the m th -mode with the selfinteraction reads as The action due to the cross-term interactions between the modes reads as The coherent drive action is And finally the action due to the Lindblad dissipator reads as The action is stationary at the saddle point, determined via ∂S k /∂α c,q = 0, which leads to α mq = 0 and the following equations of motion for α mc Comparing the EoM for α mc of the action saddle points with Eq.A1 for the MFs, one can see that α mc = √ 2 âm , which further clarifies the meaning of the classical fields in terms of the mean-fields.

Approximated Gaussian Action
As mentioned in the text, the effects of quantum fluctuations δα mc,q , can be included by expending the Keldysh action around the saddle points.Since the Hamiltonian of Eq. 1 is quartic, the expansion has terms linear -quartic in fluctuations, in general.Here, we report the explicit form of the Keldysh action up to second order in fluctuations, i.e.only the Gaussian parts

M matrix coefficients
In this appendix, we give the explicit form of the entries of the matrix M in Eq. (13) [G A (ω)] −1 reads as , where in the above equations is the MF total number of particles in the cavity.
And finally we have for [I] 2×2 being an identity matrix so rank 2.

Single mode action
We can determine the single-mode actions of each mode by Gaussian integration.Considering the steady state with empty mode a 1,3 , Figs.
where (m, m) = (1, 3), (3,1).From the poles of the Green's function we obtain the eigenvalues Interestingly, the pumped mode behaves as a single driven Kerr-oscillator whereas the symmetric modes a 1,3 , even though empty, present a non-trivial Green's function.

Spectral function
Using a Keldysh action approach we readily get access to the Green's function of the system and associated dynamical observable.In Fig. 10, we compare the spectral functions, A(ω) = −2Im[G R (ω)], of the threemodes in the LP (top) and HP (bottom) phases for dif- ferent pump strengths at fixed detuning ∆ 2 = +5, points S 1,2,3 in Fig. 4. From the top panel we see that the LP phase has the response of a "ground state", i.e., a positive (negative) peak at positive (negative) frequencies, see Sec.III.We note that peaks at negative frequencies are unresolved due to the scale resolution.On the other hand, the HP phase features a peak swap with a positive (negative) peak at negative (positive) frequencies, features of a stabilized excited state [cf.Sec.III].We note that the two uniform phases have different occupation of the 2 nd -mode, see Fig. 5 LP and HP, nevertheless their spectral functions highlight a more profound difference between them than just a different mode occupation.Interestingly, the peak swap is also present in the response of the empty side modes and we can think of it as being in presence of a normal phase to excited normal phase transition.Figure 11 compares the results of the MF-theory (blue lines) with the the full density matrix solution (red line) and the one obtained from the Gaussian approximation (green line).As can be seen the Gaussian approximation not only matches with the full density matrix solutions far away from the phase transition point but also captures the transition behavior from the LP to the HP.Specifically, the MF bistability region and the associated hysteresis behavior disappear and the two branches are connected via a rapid but continuous change.The remained discrepancy between the 2 nd -cumulant results and the density matrix within the transition region is related to the non-Gaussian nature of the cavity state due to the quartic interaction.We note that the multistable behavior does not appear in the analytic solution and the quantum solution is unique, while the semiclassical approach gives multiple dynamically stable solutions.

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

10 FIG. 2 :
FIG.2: Closed phase diagram of the three-mode uniform case.Two distinct regions are indicated by their respective meanfield energy potential ( H3) landscape as a function of the real part of the cavity field, Re(α2).The dark-blue region indicates the parameter regime where the LP phase is the ground state of the system and the HP represents a physically-allowed excited state.The LP phase disappears in the light-blue region where the HP is the only physical state.

2 FIG. 3 :
FIG. 3: (a),(b) The order parameter, |α2| 2 , and (c)-(d) the excitation spectrum, ω, on top of the LP (light hues) and HP (dark hues), along the orange and brown cut lines of Fig. 2, as a function of the pump strength, Ω2, and detuning ∆2, respectively.Real (solid) and imaginary (dashed) values with green [blue] hues encode the particle-[hole-]excitations, i.e. ds 2 > 0 [ds 2 < 0].In (a), (c), at the Ωc boundary between the light blue and dark-blue regions in Fig. 2, the LP ceases to exist and the only possible state is the excited-state HP.In (b), (d), at the ∆c boundary, the LP appears as a possible state of the system.The HP always feature a norm swap with negative norm for the positive eigenmode.The gray lines indicate the third unphysical solution.

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

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

2 FIG. 7 :
FIG. 7: Spectral function of the three-mode harmonic cavity uniform phase, for the points PLP = (∆2 = +5, Ω2 = 0.5), (a), and PHP = (∆2 = −5, Ω2 = 0.5), (b), [cf.Fig. 4].The spectral function of PHP features a peak swap with respect to the one of PLP.This signals the presence of a transition between particle-to hole-like physics going from region I to IV of the open phase diagram.

FIG. 9 :
FIG. 9: (a),(b) The spectrum A(ω) of the 3-mode cavity at PLP and PHP in Fig.4.Using 2 nd -cumulant approximation (solid lines) reveals a frequency pulling towards lower frequencies as compared to MF calculation (dashed lines).The effect is more pronounced in the high density phase where the corrections to the MF are bigger.

FIG. 10 :
FIG.10: Spectral function, A(ω), of the three-mode harmonic cavity uniform phases, for fixed detuning ∆2 = +5 and U0 = −1.Green, red and blue correspond to the 1 st -, 2 nd -and 3 rd -mode, respectively.LP phase (top) and HP (bottom) at points S1,2,3 in Fig.4.The LP shows no particular features whereas the HP shows the typical peak inversion characterizing a dissipation stabilized excited state.

2 FIG. 11 :
FIG.11:Comparison between MF (blue line), 2 nd -cumulant (green line) i.e. the Gaussian approximation, and the full density matrix (red line) solutions in a single-mode cavity at ∆ = +5.For the MF solutions, the solid and dashed lines correspond to the stable and unstable solutions, respectively.U0 = -1 in all calculations.
are the inverse of the advanced, retarded Green's functions, and the Keldysh component, respectively [cf.Appendix B 2].