Distinctive class of dissipation-induced phase transitions and their universal characteristics

Coupling a system to a nonthermal environment can profoundly affect the phase diagram of the closed system, giving rise to a special class of dissipation-induced phase transitions. Such transitions take the system out of its ground state and stabilize a higher-energy stationary state, rendering it the sole attractor of the dissipative dynamics. In this paper, we present a unifying methodology, which we use to characterize this ubiquitous phenomenology and its implications for the open system dynamics. Specifically, we analyze the closed system’s phase diagram, including symmetry-broken phases, and explore their corresponding excitations’ spectra. Opening the system, the environment can overwhelm the system’s symmetry-breaking tendencies, and changes its order parameter. As a result, isolated distinct phases of similar order become connected, and new phase-costability regions appear. Interestingly, the excitations differ in the newly connected regions through a change in their symplectic norm, which is robust to the introduction of dissipation. As a result, by tuning the system from one phase to the other across the dissipation-stabilized region, the open system fluctuations exhibit an exceptional pointlike scenario, where the fluctuations become overdamped, only to reappear with an opposite sign in the dynamical response function of the system. The overdamped region is also associated with squeezing of the fluctuations. We demonstrate the pervasive nature of such dissipation-induced phenomena in two prominent examples, namely, in parametric resonators and in light-matter systems. Our work draws a crucial distinction between quantum phase transitions and their zero-temperature open system counterparts.


I. INTRODUCTION
In equilibrium and at zero temperature, we observe socalled quantum phase transitions (PTs), which are generally described through changes in the system's Ginzburg-Landau (GL) energy functional and its symmetries as a control parameter is tuned [1], see Fig. 1(a). New minima and maxima develop in the GL energy functional when the system's symmetries are broken and the ground state of the system changes. Correspondingly, the PTs manifest as sharp changes in the system's observables (order parameters), e.g., a second-order PT corresponds to an abrupt change in the derivative of the order parameter (fluctuations). Depending on the specific potential deformation, PTs of different orders occur.
Phase transitions also transpire in out-of-equilibrium systems, where the interplay between coherent and incoherent dynamics compete in phase space to give rise to new features [2]. Specifically, the system's potential energy is extended to a phase-space energy functional and environment-induced dissipation channels act as additional forcing terms, see equilibrium case, dissipative PTs correspond to an abrupt change in order parameters corresponding to the stable attractor bifurcation topology. In this paper, we focus on the situation where the incoherent dissipative channels overwhelm a purely coherent symmetry-breaking PT.
In this paper, we introduce a general framework with which to analyze and understand the ubiquitous characteristics of dissipation-induced PTs. Specifically, we show how dissipation lifts the boundaries of equilibrium (closed) PTs and is able to stabilize states with an energy higher than that of the ground state of the closed system. We identify telltales in the closed system's potential energy and associated fluctuations' spectrum, where dissipation is prone to induce a PT. Specifically, this situation arises at regions where the closed system weakly breaks a symmetry and the dissipation channel manages to overwhelm the symmetry-breaking PT A (second-order) Z 2 symmetry-breaking phase transitions takes place as model parameters are varied and the ground state of the system changes from having an order parameter ϑ 0 to two degenerate states with order parameters ϑ 1 and ϑ 2 . (b) Out-of-equilibrium quantum phase transition in phase space, which is spanned by the coordinates θ 1,2 . Similar to the equilibrium case, phase transitions (bifurcations) occur when sharp changes occur for the attractors as the model parameters are varied. Crucial to our work, the bifurcations can occur due to a change in the potential energy, but also due to increased dissipative forcing. (c) The system is also characterizable using a (quasi-)energy functional in phase space (θ 1 , θ 2 ), where dissipative force terms are introduced via additional arrows, which can stabilize high energy states of the closed system. appearance of the dissipation-stablized phase, regions that were distinctly isolated in the closed system phase diagram become connected. These regions have seemingly identical steady states but interestingly differ in their dynamical responses. We highlight the broad applicability of dissipation-induced phenomena by studying how they arise in two prominent out-of-equilibrium systems, namely, the Kerr parametric oscillator (KPO) and the interpolating Dicke-Tavis-Cummings model (IDTC).
Our paper is structured as follows: In Sec. II, we provide a comprehensive overview on the machinery for studying in-and out-of-equilibrium systems via (a) their mean-field solution to obtain the system's fixed points, (ii) analyzing the fixed-points' stability through linearization of the equations of motion near the fixed points, and (iii) Keldysh action formalism for deriving dynamical response functions. In Sec. III, we harness the introduced machinery to highlight key signatures of dissipation-induced PTs, and demonstrate its ubiquitous phenomenology in two unrelated paradigmatic examples, namely, the KPO [19][20][21] and the IDTC model [22,23]. We conclude and discuss our results in Sec. IV.

II. CHARACTERIZING CLOSED AND OPEN SYSTEMS
Coupling a system to an environment opens it up to dissipation channels and can lead to dissipation-induced PTs, which is the focus of this paper. In such transitions, the incoherent channels can dramatically alter the system's steady state, overruling the closed system's ground-state physics, and forming new regimes defined by stable attractors of the open system.
In this section, we provide a step-by-step guide of the general framework for analyzing dissipation-induced transitions and their properties; for a more detailed overview and pedagogical approach, we refer the reader to dedicated literature [1,11,24,25]. At every step, we highlight the relevant assumptions, and in Sec. II E discuss the range of applicability of our framework.
In Sec. III, we use this formalism to study the impact of dissipation in the presence of a Z 2 symmetry breaking. Specifically, we discuss two seemingly unrelated physics problems, a parametrically driven Kerr resonator and a coupled cavity-atom system. We show that in both cases, the dissipation can (i) overwhelm the symmetry breaking, thus (ii) connecting regions in parameter space that are otherwise disconnected, where (iii) an excited state of the system characterized by a negative mass instability is stabilized.
In Sec. II A, we succinctly review the analysis of ground and excited states of a closed system, based on Landau-Ginzburg energy functionals and fluctuation spectra [2]. The latter is particularly useful to help identify physical states that can be stabilized by dissipation. In Sec. II B, we similarly review the analysis of stable and unstable steady states of open systems in appropriate rotating pictures, using mean-field quasienergy functionals (i.e., the energy in a rotating phasespace frame), and dissipative Bogoliubov excitation spectra [26,27]. The latter results from the linear stability analysis of steady states. As we shall see in Secs. III A and III B, in concomitance with a Z 2 symmetry-breaking dissipation-induced transition, an exceptional pointlike scenario emerges, i.e., the Bogoliubov spectrum exhibits a pair of complex-conjugated eigenvalues whose imaginary parts go to zero while the real parts split [15,28]. Such behavior similarly marks squeezing in the fluctuations of the system.
Finally, using a Keldysh action approach [11,12], we derive frequency-resolved observables, such as the spectral function [A(ω)], see Sec. II D. The response functions present peaks corresponding to fluctuation resonances of the system. These resonances correspond to poles that coincide with the dissipative Bogoliubov excitation spectra, described in Sec. II B. As we show in Secs. III A and III B, a Z 2 symmetry-breaking dissipation-induced transition connects regions in parameter space, with steady states characterized by the same order parameter, but with different dynamical responses. This is traced back to the fact that the states in the disconnected regions of the closed system correspond to ground or excited states, and the dissipation smoothly connects between them.
For pedagogical clarity, we use the example of a simple harmonic oscillator to illustrate our methodology. We refer to the harmonic oscillator Hamiltonian where we seth = 1 and ω 0 is the characteristic frequency of the oscillator. The operatorâ annihilates an excitation in the oscillator, and we rewrite the Hamiltonian in matrix form for convenience.
FIG. 2. Energy quantization in a parabolic potential. (a) A harmonic potential can also be written as a sum over a particle (blue) and a hole (red) potentials, cf. Eq. (1), corresponding to the two basis solutions of the problem. (b) Changing the potential, a Z 2 symmetrybreaking transition can occur, where two minima separated by a maximum appear, each with their own excitation spectrum.

A. Closed system
We first consider the equilibrium physics of closed systems, where the Hamiltonian encloses all the information describing the system. It is commonly difficult to solve the Hamiltonian exactly, and we present here a mean-field treatment of the problem and a study of the fluctuations around it.

Mean-field solutions
When studying a system comprised of many degrees of freedom, e.g., an ensemble of harmonic oscillators, we study the mean-field (Landau-Ginzburg) energy functional, defined forâ i → α i , where α i is a complex number corresponding to the semiclassical limit (coherent state) of the operatorâ i , and i iterates over all N degrees of freedom of the system. An inspection of the resulting N-dimensional complex functional provides insightful information on the states of the system. Specifically, the extremal points of the functional α j = (α 1, j , α 2, j , . . . α N, j ) correspond to the ground state of the system, as well as to possible excited states. The ground state is defined as the lowest-energy state, see Fig. 1(a). For the harmonic oscillator Eq. (1), the Landau-Ginzburg energy functional is parabolic,H = ω 0 |α| 2 , with only one global minimum at α 0 = 0, see Fig. 2(a).

Excitation spectrum
Having identified the extremaα j of the Landau-Ginzburg energy functional, each of these extrema are characterized by an excitation spectrum defining the fluctuations around the extremal point. This spectrum is obtained via a substitutionâ i = α i, j + δâ i, j into the original Hamiltonian and by retaining terms up to second order in the fluctuation operators δâ i, j =â i − α i, j . Then, using a Bogoliubov transformation, we diagonalize the resulting respective quadratic excitation Hamiltonians,Ĥ exc , see Eq.
(2) below. The mean-field solutions are physical only if their associated (Bogoliubov) excitation energies are real. Diagonalizing the excitation Hamiltonian, we ensure that the diagonalization procedure preserves the statistics of the excitations. To this end, the proper Bogoliubov transformation [24,29,30] for bosonic systems is found by diagonalizing the dynamical matrix,

Norm inversion in the excitations
The excitations' eigenvectors further provide us with the symplectic norm [30] associated to each excitation eigenfrequency, where v j with j = 1, . . . , 2N are the eigenvectors of D exc . The symplectic norm Eq. (3) determines the nature of the excitations: It can be either positive (ds 2 > 0) or negative (ds 2 < 0) and is a measure of the particle-or holelike nature of the excitation, respectively [24,26,30]. For instance, the ground state of the system as well as other local minima of the energy functional are accompanied by particlelike (holelike) eigenmodes at positive (negative) frequencies, whereas local maxima are associated with a holelike (particlelike) eigenmode at positive (negative) frequencies. In Fig. 2(a), we sketch the particle-(blue) versus holelike (red) ladders for the harmonic oscillator excitations around the ground state. The distinction between particle-and holelike modes is crucial and it assists us in identifying excited states that correspond to maxima in the mean-field energy functional, which can become a stable steady state due to a dissipation-induced PT.
As we tune an external parameter, the shape of the meanfield energy functional can change and develop new features. For example, a (second-order) spontaneous Z 2 symmetrybreaking PT can occur when the energy landscape moves from exhibiting a single ground state to two degenerate ground states separated by a maxima, see Fig. 2(b). In such a scenario, the symmetry-conserving mean-field solution now becomes an excited state. As the system traverses the critical point, the symplectic norm of the excited state changes signs. In the closed system setting, this state cannot be reached by adiabatic ground state evolution; we will show later that dissipation can stabilize such a steady state.

B. Open system
To characterize the physics of open systems, we first find the so-called steady (or stationary) states, which correspond to the long-time behavior of the open system. To do so, we solve Liouville's master equation [2], in the long-time limit. Here ρ is the density matrix of the system and L is the Liouvillian superoperator. We obtain the 023100-3 steady states by requiring that the system's density matrix does not change in time, i.e., dρ/dt ≡ 0. Note that rotating steady states, e.g., limit cycles, fulfill a similar condition in the correct rotating frame.
In this paper, we consider cases with a strong separation of time scales between system and environment, and assume a memoryless environment with weak system-environment coupling, such that we can perform the Born-Markov and secular approximations, thus simplifying Eq. (4) to obtain Lindblad's master equation [25,31,32]: The unitary evolution generated by the system's Hamiltonian is described through the commutator, while the dissipative dynamics with rates γ α appear in Lindblad form, where the Lindblad dissipator is defined in terms of so-called Lindblad operators (or quantum jump operators)L α that describe the coupling to the environment. The anticommutator term corresponds to dissipation, while the so-called recycling or quantum jump term, 2L α ρL † α , encodes fluctuations and ensures the normalization Tr{ρ} ≡ 1 of the system's density matrix.
In the context of our pedagogical example Eq. (1) and throughout this paper, dissipation due to coupling to a zerotemperature environment takes the Lindblad form where κ ω 0 is the positive dissipation rate describing loss of excitations. We can pictorially describe the harmonic oscillator evolution under Eq. (5). The eigenstates of the closed system describe closed circular orbits along equipotential lines of the Landau-Ginzburg energy functional. We dress the Landau-Ginzburg energy functional with arrows whose directions depend on the particular dissipation channel. As the dissipation Eq. (7) destroys excitations in the system, the arrows point toward the empty resonator state. In other words, the dissipation pushes the closed system orbits inward. The net result is a spiraling evolution toward the center, which is the sole steady state for the dissipative harmonic oscillator, see Fig. 3(a).

Mean-field steady states
The master equations, Eqs. (4) and (5), evolve the state of the system. To better characterize its steady state, we use instead the time evolution of mean-field expectation values in Schrödinger's picture, where ô is the expectation value of one of the observable operators describing the system. Considering different observables, we obtain a set of coupled first-order differential equations,  (7). (a), (b) The rotating mean-field energy potential landscapeH rot (cf. Sec. II C) for negative and positive detuning . The arrows indicate the resulting Hamiltonian motion. When changes its sign, the peak in A rot passes 0,H rot changes from a valley to a hill, and the rotation around the origin (center) changes its direction. Adding dissipation (arrows) leads to a flow toward the empty state at the center.
where O = ( ô 1 , ô 2 , ô 3 , ...) is the basic set of real operator expectation values describing our system. Note that for bosonic coherent states α = â is a complex number and we split its evolution into its real and imaginary parts in Eq. (9). Taking the steady state, dO/dt ≡ 0, we solve the resulting set of algebraic equations with multiple solutions that satisfy MO = 0. We find the set of steady-state solutions and henceforth denote by O j i the expectation values corresponding to the jth physical solution of the ith operator. Note that physical solutions have real observables, i.e., Im{O j i } = 0, whereas Eqs. set (9) can yield unphysical (complex) solutions.

Stability analysis
Among the physical steady states, some are stable and some are unstable against linear fluctuations, cf. the above discussion on excitations in the closed system. Depending on the initial boundary conditions, the system will eventually evolve toward one of the former. To find which steady states are stable or not, we perform a linear expansion of Eq. (9) around each solution according toô We then define the stability matrix M j , whose entries depend on the jth solution, and which defines the fluctuation dynamics around this solution through where δÔ j = ( δô j 1 , δô j 2 , δô j 3 , . . .) is the set of the fluctuation operators' expectation values describing our system. Equation (10) defines a set of coupled linear differential equations, where the eigenvalues i of the stability matrix M j determine the stability of the jth solution. This results from the fact that the eigenmodes of M j evolve in time according to e t . Then, the real part of plays the role of a lifetime of the excitation mode, whereas the imaginary part is its frequency. In other words, if at least one eigenvalue has a positive real part, the fluctuations diverge and the solution is unstable.
. Correspondingly, the eigenvalues show a mode softening but no unstable behavior, i.e., the imaginary parts go to zero for ω 0 = κ/2, while the real parts always remain negative but split from the value −κ. The variance associated with the real part is flat along the entire region, while the imaginary part variance increases but remains finite. (b) Same as (a) for the harmonic underdamped oscillator in the rotating frame. In the weak dissipation limit, as a function of detuning , the imaginary part of the fluctuation eigenvalues inverts at = 0, while the real parts remain negative and unchanged. The variance remains unchanged, marking the fact that no phase transition is occurring. (c) Underdamped harmonic oscillator response (spectral function) [cf. Eq. (15)] for positive (light green) and negative (dark green) detunings.

Variance
In similitude to equilibrium statistical physics, PTs in open systems can be classified through discontinuities in the properties of the order parameters and their fluctuations [1]. Specifically, first-order PTs manifest in a discontinuous change in the mean-field order parameters, second-order PTs show a discontinuity in the fluctuations of the system, whereas third-order PTs exhibit a discontinuity in the variance of the fluctuations. In the following sections, we show that the dissipation-induced PT connects between two steady states, which were disconnected in the closed system limit, but are characterized by the same mean-field order parameter. Moving between the two newly connected states, a discontinuity in the system's dynamical behavior appears but no discontinuity manifests-neither at the level of fluctuations nor for the variance of the fluctuations.
To show this, in the following, we calculate the covariance of the fluctuations using Eq. (8) for the time evolution of the expectation of correlation operators: Taking the lowest-order (mean-field) results for the covariance operators yields systems of coupled equations analogous to Eq. (9), which we solve for the steady-state.

C. Static observables of the harmonic oscillator
We are now ready to apply our aforementioned machinery onto the dissipative harmonic oscillator introduced in Eqs. (1) and (7). We will explore both the overdamped (ω 2 0 < κ 2 4 ) and underdamped (ω 2 0 > κ 2 4 ) regimes, as they allow us to introduce key signatures that will be useful for our discussion in Secs. III A and III B.

Overdamped oscillator
In the overdamped limit, dissipation can overwhelm the Hamiltonian potential of the system. Note that this limit involves a system-environment coupling that goes beyond the Lindblad limit discussed above [33,34]. Nevertheless, a sim-ilar mean-field, stability, and variance analysis as in Eqs. (9) and (10) holds here, see Appendix A. There is a single steadystate to the system, namely, an inert oscillator. In Fig. 4(a), we show the fluctuation eigenvalues and variance obtained on top of the empty state. We observe that at a certain point, the eigenvalues show an abrupt change (akin to an exceptional point [28]) in their derivative with respect to the control parameter (here the mode's eigenfrequency). Specifically, the imaginary parts become degenerate, while the real parts split in a region in parameter space. This marks the fact that the fluctuations become overdamped and do not exhibit an oscillating behavior. Instead, two characteristic decay times appear [dashed lines in Fig. 4(a)] and lead to squeezing of the fluctuations. Note that the sudden change in the fluctuations eigenvalues does not result in a discontinuous variance, ruling out the presence of a third-order PT.

Underdamped resonator
In open systems, we often deal with the presence of an external drive. As the system is expected to lock to the drive frequency in the long-time limit, a useful method to treat external drives is to move to a rotating frame with respect to the drive. If the oscillator is driven at a frequency ω drive , with the unitary transformation T = e −iω drive ta † a we move to a rotating frame at frequency ω drive , yielding the rotating frame Hamiltonian H rot = − a † a with detuning = ω drive − ω 0 . The sign of the detuning determines whether the drive stiffens or softens the harmonic potential and, additionally, whether the mean-field energy potential landscape,H rot = − |α| 2 , is a paraboliod with a minimum or a maximum around the phase-space origin, for negative or positive , respectively [see Figs. 3(a) and 3(b)]. Depending on the detuning, the empty cavity state appears in the rotating frame, either as the ground state (negative detuning) or as an effective excited state (positive detuning). Additionally, the change in curvature ofH rot as a function of the detuning manifests as the transformation of particle-to holelike excitations and vice versa through a change in the sign of the symplectic norm Eq. (3). The switch in the norm corresponds to clockwise (counterclockwise) rotations of the state of the system relative to the rotating reference frame (clock), see Fig. 3.
The effect of dissipation remains unchanged in the rotating frame picture and decreases the amplitude of oscillations, namely, it always creates a dragging force toward the center. The joint effect of coherent and incoherent forces leads to a spiraling evolution of the state of the system toward the minimum (maximum) of the potential in the case of negative (positive) detuning, see Fig. 3. The quasienergy potential inversion also manifests in the open system fluctuation spectrum in the form of a degeneracy point, where a pair of complex conjugated fluctuation eigenvalues, i , acquires a zero imaginary part while the real parts remain negative. The latter condition implies that the system's steady state remains stable at = 0 despite the fact that the effective confining Hamiltonian potential vanishes. At the same time, the variance as a function of detuning remains unchanged, marking the fact that the degeneracy point at = 0 does not imply a PT, see

D. Dynamical observables
In the closed system, the sign of the fluctuations' symplectic norm Eq. (3) helps us to distinguish between the ground state (positive norm) and an unstable excited state (negative norm) even if the two states are characterized by the same mean-field observable expectation values. In the open system scenario, we do not have the notion of a symplectic norm. We resort instead to study dynamical observables of the system to fully understand the nature of a steady state.
To study dynamical observables, we use a Keldysh action formulation, which for a bosonic theory takes the form [11,12] where the fields a c , a q are the standard classical and quantum bosonic fields acting on the Keldysh contour. This terminology signals that the former combination of fields can acquire a (classical) field expectation value, while the latter one cannot. The details of the derivation can be found in Ref. [12]. The strength of this formalism lies in its matrix structure, where the action Eq. (12) is directly connected to the Green's functions of the system, which are the retarded and advanced Green's functions, as well as the Keldysh component, respectively. It is useful to move to Fourier space and obtain a frequency-dependent action, where ω = dω 2π . Taking, for example, the harmonic oscillator Eq. (1), we can write the action Eq. (13) as where we assumed a Lindlad dissipator in the form of Eq. (7). We identify the inverse Green's functions and the Keldysh component with Inverting the matrix action Eq. (12), we obtain direct access to the frequency-resolved Green's functions of the system; for the harmonic oscillator Eq. (1) Using these Green's functions, we can calculate frequencydependent observable quantities (dynamical order parameters), such as the response and correlations of our system. Here, we provide a short list of those correlations that are most-commonly probed: (1) The spectral function (density) (A).
In linear response theory, the spectral function describes the excitation that the system undergoes when adding a single particle with frequency ω to it. It is defined as where the retarded Green's function, G R (ω), describes the linear response of a system to a weak external perturbation. For bosonic systems, A(ω) fulfills the sum rule: (2) The power spectrum (C). The power spectrum describes the occupation of the individual quantum excitation modes of the system. It is defined as the Fourier transform of the correlation function {a(t ), a † (0)} , yielding where G K (t, t ) = −i {a(t ), a † (t )} . Considering the steadystate solution, and taking the equal-time limit, we can find the mode occupation via (3) The fluorescence spectrum (S). The autocorrelation of the systems' excitations describes the probability of measuring an excitation of frequency ω leaving the system, dubbed fluorescence spectrum. It is defined as Each one of these functions presents peaks of different heights at the resonance frequencies of the system under study. In the simple case of the dissipative harmonic oscillator, we find that the spectral function exhibits the familiar line shape peaked at the resonance frequency of the oscillator: Interestingly, the mode occupation C(ω) in this system coincides with the spectral function as can be seen by comparing A(ω) and G K (ω). This will not be the case in more complicated systems.
To conclude this section, we consider again the harmonic oscillator in a rotating frame due to the presence of an external drive. The response of the oscillator can be found in the same fashion as before with the only difference that now the resonance frequency is the detuning from the drive: This implies that the oscillator response is peaked at positive frequencies for negative and at negative frequencies for positive , as shown in Fig. 4(c). In this example, the empty (α = 0) steady state has a different dynamical behavior depending on the sign of , where the two different regions are separated at = 0. Two key elements emerge from this simple discussion. First, the steady-state of the system does not change, i.e., it is always the empty state with no excitations in the system. Second, as we tune the external parameter , the dynamical response response shifts from positive to negative frequencies.

E. Limitations and applicability
The framework we outlined in this section involved certain assumptions, namely, (i) we treat bosonic systems or systems that are mappable to bosonic ones, e.g., via the Holstein-Primakoff transformation [35], (ii) we employ a mean-field (saddle-node) treatment of the order the parameter of the system on top of which quadratic fluctuations appear, and (iii) we truncate strong correlation between the system and environment (Lindblad approximation), where the environment is taken to be at zero temperature. Concerning (i), our analysis can be readily extended to fermionic or mixed systems, where the particles' commutation relations will impact the specific form of the derived expressions. The assumption (ii) is well justified whenever the system is well within a phase. Expanding beyond mean-field such that higher-order correlations are taken into account can play a crucial role whenever the system is close to criticality. Recently, several approaches have been developed to study higher-order correlations, including dynamical mean-field theory, diagrammatic expansions, functional renormalization group, exact diagonalization, and numerical or density matrix renormalization group analysis [36][37][38][39][40][41][42][43][44]. The approach presented in this work can be in principle extended to include higher order fluctuations. Going beyond (iii) can lead to interesting strongly-correlated states between system and environment, e.g., the Kondo effect [45]. Systematic diagrammatic expansion of the system environment coupling is challenging and depends crucially on the environment memory, size, and coherence, see e.g., Ref. [32] and discussion therein. We limit ourselves here to the simplest (and ubiquitous to light-matter systems) case to highlight the profound impact that even such a simple environment has on out-of-equilibrium PTs.

III. DISSIPATION-INDUCED PHASE TRANSITIONS
Following the pedagogical introduction, we are at a vantage point to better define dissipation-induced PTs. Generally, this includes any scenario where incoherent terms overwhelm coherent Hamiltonian terms and lead to a change in the steady-state topology in phase space, i.e., to bifurcations. Specifically, we now focus on the case where the Hamiltonian prompts a spontaneous Z 2 or U (1) symmetry breaking [cf. Figs. 1 and 2]. Dissipation can overwhelm the expected symmetry breaking and stabilize non-symmetry-broken states. A transition into this region manifests in the form of an exceptional pointlike scenario for the excitation spectrum and squeezing of the variance, as discussed in Fig. 4(a). More importantly, the dissipation-stabilized region bridges between two disconnected parts of the closed system phase diagram. As one traverses these regions, imaginary parts of some modes in the fluctuation spectrum invert, as discussed in Figs. 4(a) and 4(b). Correspondingly, the dynamical response functions display inversions in some of their peaks. In the following, we discuss two seemingly unrelated paradigmatic drivendissipative systems that exhibit such phenomenology. Thus, we highlight the ubiquitous nature and characteristics of such transitions.

A. Parametric Kerr oscillator
The first exemplary system we consider is the Kerr oscillator driven by a parametric (two-photon) pump, while subject to single-photon dissipation [21,[46][47][48][49][50][51]. As a function of the two-photon drive, this KPO exhibits a continuous, Z 2 time-translation symmetry-breaking transition from an empty cavity state, dubbed normal phase (NP), to a state with finite photon number, dubbed parametric phase state (PPS) [52]. The symmetry-broken PPSs appear in pairs of coherent states with equal amplitude but π -shifted phase [53]. At low photon numbers, their superpositions form Schrödinger cat states of opposite parities [20,54,55]. Possible applications of such phase states include annealing-based optimization algorithms in classical and quantum KPO networks [56][57][58], as well as universal quantum computation [53,59]. Furthermore, the competition between single-and two-photon drives leads to PTs that can be used for sensing [21,[60][61][62]. At weak two-photon drives, the dissipation stabilizes the KPO, while squeezing its fluctuations, with various applications in sensing [63][64][65]. We now turn to review how the single-photon dissipation changes the phase diagram of the KPO according to the dissipation-induced phenomenology described above, cf. Ref. [66].

Closed system
We first consider the (closed system, nondissipative) KPOs, with Hamiltonian where a annihilates a bosonic particle on the oscillator and U is the Kerr nonlinearity amplitude. A parametric two-photon coherent pump of strength G and frequency ω G is described We use the unitary transformation T = e −ia † aω G t/2 to move to a rotating frame with respect to half the parametric drive frequency. Neglecting fast oscillating terms, we obtain the time-independent KPO Hamiltonian in the rotating wave approximatioñ where = ω G /2 − ω c is the half-pump cavity detuning. In the following, we will focus on positive U , but similar physics is obtained for negative U , using the transformation → − and G → −G.
The time-translation symmetry of the Hamiltonian in the nonrotating frame maps to a Z 2 symmetry in the rotating frame (associated with a → −a) [67]. The symmetry is spontaneously broken as we tune either the pump strength G or the detuning . In Fig. 5(a), we draw the corresponding meanfield rotating energy potential landscape, obtained by substituting the mean-field ansatz a = α Re + iα Im into Eq. (21). The ground state is obtained by minimizing the Hamiltonian Eq. (21) and calculating the order parameter |α| 2 , see Appendix B for details. Considering the mean-field energy potential landscape, the phase diagram comprises three qualitatively distinct regions in parameter space with I, a paraboloid potential, where the Z 2 -symmetry is preserved and the NP is the ground state; II, a double-well potential separated by a saddle, where the Z 2 symmetry will be spontaneously broken and either of the PPSs, which are phase shifted by π , becomes the ground state; and III, a double-well potential separated by a potential hill (maximum), where the PPSs are still the ground states of the system. Crucially, in III the NP is a physically excited state of the system, while in II the NP is unphysical, i.e., its fluctuation spectrum presents complex frequencies, cf. Sec. II A. In Fig. 5(b), we plot the order parameter |α| 2 and show that it acquires a finite value within regions II and III, corresponding to the so-called parametric instability [19].
To corroborate this statement, we follow the procedure outlined in Sec. II A and find the excitation spectrum associated with the Hamiltonian Eq. (21). Here, we present the excitations on top of the NP, while the PPS case is detailed in Appendix B. The dynamical matrix associated with the NP reads [cf. Eq. (2)] In the NP, the cavity is empty and therefore there is no contribution from the Kerr nonlinearity. The eigenvalues of the dynamical matrix describing the excitations are In Fig. 5(c), we plot the excitation spectra on top of the NP and PPS, see also Appendix B for additional details. The NP spectrum becomes fully imaginary along the dashed line cut in Fig. 5(a) when entering region II (at c1 ), and at the same time the PPS spectrum becomes real, signaling the onset of the PT. Increasing the detuning further, a transition to region III occurs (at c2 ), and the NP excitation spectrum becomes real once more. This occurs in concomitance with a change in the sign of the symplectic norm of the excitations, cf. Eq. (3) and see Appendix B for the derivation of Eq. (26). A positive symplectic norm is a signature of particle-dominated fluctuations, whereas a negative one corresponds to holelike fluctuations. We therefore observe a swap in the norm of the two NP excitations between regions I and III of the phase diagram, i.e., the white and light-blue regions in Fig. 5(a). The two NPs are dynamically different, reflecting that in region III the NP is an excited state of the system. We dub it an excited normal phase (e-NP). As we show later, in the presence of dissipation, the e-NP becomes a stable steady state of the system. We here note that excitations on top of the PPS retain their respective norms throughout the entire parameter space, i.e., the positive frequency has a positive norm, whereas the negative one has a negative norm.

Open system
Following the methodology introduced in Sec. II B, we now study the impact of dissipation on the KPO Eq. (21). We introduce single-photon losses in the form of the same Lindblad dissipator as in Eq. (7) and study the master equation: We proceed and perform a mean-field approximation and characterize the steady state of the system via the expectation value of the operatorâ, i.e., via α = a . The equation of motion describing the evolution of α is given by We look for the steady state dα/dt = 0, and find the solutions for the order parameter α, see Appendix C for the analytic expressions of the solutions. In Fig. 6(a), we plot the order parameter behavior as a function of detuning along the orange cut in Fig. 6(c). From negative to positive detuning, we start from a parameter regime where only one stable solution is possible, i.e., the trivial NP. We then enter regions where up to five solutions for the steady-state mean-field equations are possible. Due to the Z 2 symmetry of the problem, the nonzero amplitude solutions have pairwise the same amplitude, but are π -shifted in phase. We therefore focus only on three representative solutions.
As before, we perform a stability analysis against linear fluctuations and study the eigenvalues of the possible solutions, cf. Sec. II B and Appendix C. As mentioned at the end of Sec. II B, a solution is unstable if at least one of the fluctuation eigenvalues has a positive real part. In Fig. 6(b), we plot the fluctuation eigenvalues associated with the various physically relevant steady-states. We observe exceptional pointlike sce-narios both for the NPs and PPSs. The NP becomes unstable at c3 where a pitchfork bifurcation takes place and stable PPSs appear, denoted PPS st . Increasing the detuning further, we reach the point c4 where the low-population solution becomes stable again (e-NP) alongside additional high-density unstable solutions that appear, denoted PPS un . The former high-density solutions stay stable, and we have a region of costability between the e-NP and the PPS st .
We extend our study to the whole parameter space and draw the open steady-state phase diagram in Fig. 6(c). Similar to the closed system case, the phase diagram comprises three different regions: (I) the white region presents only the NP as a stable steady state, (II) the blue region only the PPS st , and, finally, (III) the light-blue one is a region of costability between the e-NP and the PPS st . At the same time, the open system phase diagram is strikingly different with respect to the closed system one. Specifically, dissipation overwhelms the PPS tendency of the system, and shifts regions II and III upward toward stronger drives G, see Fig. 5(a). In turn, new regions appear: In II', the NP is stabilized despite the fact that it does not correspond to a physical state of the closed system and, in III', the NP remains the sole stable solution while being an excited state of the system. The appearance of region II' is what we dub a dissipation-induced PT, which in turn directly connects two distinct phases (regions I and III) in the closed system. As we shall see in the following, as the dissipation stabilizes the e-NP in a region of the parameter space where it is an excited state of the closed system, it is a dynamically different NP with respect to the NP at negative detuning.
We conclude this subsection by showing that the NP → e-NP transition, involving exceptional pointlike discontinuities in the fluctuation spectrum, is not a third-order PT. Specifically, we show that the variance of the fluctuations [cf. Eq. (11)] are differentiable. Without loss of generality, we consider the brown cut in Fig. 6(c), and study the fluctuation variance over the NPs. The results are shown in Figs. 6(d)-6(e). We identify the transition region, c5 < < c6 , as a region of overdamped excitations [cf. Sec. II C and Fig. 4(a)], where the real parts of the eigenvalues, Re( ), split but remain negative and the variances remain finite while changing continuously. The variance is continuous despite the discontinuity in the eigenvalues due to a compensation coming from their associated eigenvectors, see Appendix C for details. Note that there is a crucial difference between the overdamped oscillator in Fig. 4(a) and the KPO in Fig. 6. In the former case, the overdamped region appears as a consequence of strong dissipation, whereas in the latter it appears due to the interplay between a weak dissipation channel [Eq. (7)] and the effect of the parametric drive [Eq. (21)], leading to a locking of the motion to the drive, which manifests as overdamped fluctuations.

Keldysh KPO
We now employ the formalism introduced in Sec. II D and show that the two NPs, that are now connected via the dissipation-stabilized region II', are dynamically different. We first write the Keldysh action Eq. (12) for the KPO, where we did not resort to the matrix notation due to the presence of the quartic interaction term. As is customary in the Keldysh path integral formulation [11,12], we perform a saddle-point approximation minimizing the action with respect to the quantum field: We then set a c (t ) = √ 2α, a q = 0, and solve for the steady state, i.e., ∂ t a c = 0. We note here that such a solution is the same on the two branches of the Keldysh contour and coincides with the mean-field solution of Eq. (28) for the cavity field.
We now study the fluctuations on top of the mean-field solution and address the spectral properties of our system. We do so by expanding around the mean-field solutions where α is either one of the NPs or PPSs and δa cq are the fluctuation fields. We insert Eq. (31) into the action Eq. (29) and retain up to second order in the fluctuation fields to obtain the fluctuation action where we introduced the spinor notation φ cl,q = (δa cl,q , δa * cl,q ) T and defined the inverse Green's functions and Keldysh component as with I 2×2 the identity matrix of dimension 2. We now study the spectral properties of the NPs and of the PPSs.

Normal phase
In the NPs, we have a zero mean-field solution, α = 0. We invert Eq. (32) with Eqs. (33) and (34) to obtain the Green's functions We consider two distinct points in the parameter space, 1 , and 2 [cf. Fig. 6(c)], where the NPs are stable steady states. The former lies in the region of negative detuning, where the NP coincides with the ground state of the closed system. The second one instead is in the positive detuning region, where dissipation stabilizes the NP resulting in the e-NP. We calculate the spectral function Eq. (15), the power spectrum Eq. (16), and the fluorescence spectrum Eq. (18). We find that the first two observables present a peak inversion between the two different regions in the parameter space whereas the last is perfectly symmetric, see Fig. 7(a). The peak swap between the two region bears important physical significance and underlines the fundamental dynamical difference between the NP and the e-NP. It is reminiscent of the particle-hole inversion in the excitation spectrum of the NP already highlighted in the closed system. More specifically, it also marks the change in oscillation chirality relative to the reference frame [cf. Fig. 3] by the fact that a negative peak appears at positive response frequency in the spectral function, A(ω) [30]. Ultimately, dissipation stabilizes the NP in region II , where it is not a physical solution of the closed system and therefore does not exhibit proper excitations, as the excitations are overdamped. In region III , an otherwiseinaccessible excited state is reachable by adiabatic ground  Fig. 6(c). (b) The responses on top of the PPSs at points 2 , 3 in Fig. 6(c). The spectral response and power spectrum of the NPs feature a mode inversion when going from the NP to the e-NP, 1 → 2 . This is a signature of the dynamically different nature of the two normal phases in the different regions of parameter space. On the contrary, for the PPSs, no qualitative change in the spectra is observed-neither when crossing from negative to positive detuning nor when considering the costability region. state evolution from negative to positive detuning. At the same time, the system response retains features of an excited state, signaled by the holelike nature of the dynamical fluctuations (peak inversion).

Parametric phase state
The PPSs are characterized by a finite cavity field α, see Appendix C. We invert Eqs. (33) and (34) once more with the new mean field and obtain the Green's functions for the PPSs. As for the NPs, we consider two points in the parameter space, 2 , and 3 in Fig. 6(c). The point 3 lies in region II, where only the PPS st are stable steady states, whereas 2 is positioned in the coexistence region III. In Fig. 7(b), we plot the dynamical response functions on top of the PPS st at 2 , and 3 . There are no qualitative changes in the spectra as opposed to the NPs, which undergo a peak swap between the two regions in parameter space. This agrees with the discussion above on the differences that the closed system experiences between the NP and e-NP, whereas the PPSs do not present a particle-hole fluctuation inversion.

B. Interpolating Dicke-Tavis-Cummings
In the previous section, we studied the KPO as a ubiquitous model describing a broad class of systems, where one directly drives a bosonic resonator (be it vibrational, electric, photonic, or composed of more complex particles) with a two-particle drive. In this section, we consider instead a many-body lightmatter system, where PTs appear due to the coupling between the matter degrees of freedom and a bosonic resonator cavity [23,[68][69][70][71][72][73][74][75][76]. We specifically study a ubiquitous model describing light-matter systems, the IDTC model [22,30,77,78], which is a generalized version of the Dicke [68,79,80] and the Tavis-Cummings [81] models. The model system comprises of a single-mode bosonic cavity (light) coupled to the x and y components of N spinlike (two levels) degrees of freedom (matter). As a function of the coupling between the cavity and the spins, the IDTC model features transitions between an empty cavity state, the NP, to a highly occupied cavity state, dubbed superradiant phase (SP) [82,83].
In similitude to the KPO above, dissipation stabilizes and extends the NP into a new parameter regime, dubbed e-NP [23,30]. This dissipation-induced PT also allows for regions of coexistence between the e-NP and SP [23,84]. The analogy between the KPO and IDTC is even more explicit when considering the fluctuation spectrum, which shows a soft-mode invertion along the low-density dissipation-facilitated NP → e-NP transition. We highlight here these distinctive features of the dissipative IDTC and draw the universal parallels between this many-body light-matter system and the KPO while employing the methodology introduced in the previous sections.

Closed system
We start our discussion with the closed system IDTC Hamiltonian [22,23,30,85] H I =hω c a † a +hω z S z where a is the bosonic annihilation operator of the cavity field, and ω c is the cavity's resonance frequency. The ensemble of N two-level system is described in terms of collective spin operators S x,y,z = N i=1 σ {i} x,y,z , with Pauli matrices σ {i} j describing theith two-level system. Due to the tunable coupling between the spins and the cavity, this model interpolates between the Tavis-Cummings (λ x = λ y ) [81] and Dicke (λ x = λ y ) [68,79,80] models. The IDTC Hamiltonian Eq. (36) has a Z 2 × Z 2 symmetry corresponding to four different superradiant states with finite spin magnetization either along the x or y directions. Furthermore, tuning the couplings to be equal, λ x = λ y , boosts the symmetry to become U (1). In either case, increasing the coupling above a critical coupling λ c = √ ω c ω z /2 leads to spontaneous breaking of the symmetries, into a SP with finite in-plane magnetization [83,86]. This is shown in the phase diagram in Fig. 8(a). The mean-field energy potential landscape lies over the phase space of the resonator and the Bloch sphere of the collective spin. Focusing on the south pole of the Bloch sphere, where the NP lies, it suffices to draw the impact of the spins on the mean-field energy functional of the resonator, see Fig. 8(a). The phase diagram comprises of three qualitatively distinct regions in parameter space: I, a paraboloid potential, where the Z 2 × Z 2 /U (1) symmetry is preserved and the NP is the ground state; II, a double-well  Fig. 5), the system has Z 2 symmetries and is characterized by three different regions; with ω c = ω z and schematic illustration of the mean-field energy potential as a function of the real and imaginary parts of the cavity field, α Re and α Im , respectively. Quantum phase transitions occur between a normal phase (white) and a symmetry broken phase (blue). Analogous to the KPO, the light-blue region indicates a parameter regime where the normal phase is an excited state of the system and the SPs represent the ground state. (b) The order parameter, |α| 2 , and (c) the excitation spectrum, ω, on top of the NP (SPs), as a function of the coupling along the orange-dashed cut line in (a). In (c), real (solid) and imaginary (dashed) parts on top of the NP (green) and SPs green (blue), dark (light) hues encode the particle (hole) excitations, i.e., ds2 > 0 [ds2 < 0]. The excitation behavior is analogous to the KPO. At the boundary, λ c , between the white and blue region in (a), |α| 2 acquires a finite value in (b). Concomitantly, the eigenfrequencies on top of the NP ω NP I [cf. Eqs. (39)] in (c) become fully imaginary (dashed gray lines), thus signaling the onset of the symmetry-breaking transition. When λ > λ c1 the eigenvalues return to be real-valued, marking the transition to the coexisting region. Mirroring the KPO, the NP excitations exhibit a norm swap from λ < λ c to λ > λ c1 . potential separated by a saddle, where the Z 2 × Z 2 symmetry is spontaneously broken either along x or y directions and either of the positive or negative magnetization states becomes the ground state; and III, a double-well potential separated by a potential hill (maximum), where the SPs are still the ground states of the system. Crucially, as in the KPO, in III, the NP is a physical excited state of the system, while in II, the NP is unphysical (imaginary excitation frequencies).
Following the procedure outlined in Sec. II A, we find the ground state, order parameters α, X, Y , and spectrum of excitations of the IDTC model, Eq. (36). In the top panel of Fig. 8(b), we plot the order parameters |α| 2 and show that it acquires a finite value within regions II and III, corresponding to the superradiant PTs; this parallels the NP → PPS of the KPO, cf. Fig. 5. To show the existence of the e-NP phase and draw a comparison with the KPO studied above, see Sec. III A 1, we consider the IDTC's normal phase and analyze its excitation spectrum. Details regarding the SP can be found in Ref. [30]. The NP fluctuation Hamiltonian and associated dynamical matrix for the IDTC model are therefore, we find the closed system eigenvalues [30], cf. Sec. II A, Unlike the KPO, the IDTC exhibits two excitation eigenmodes. We focus on the soft-mode [with the negative sign in Eq. (39), plotted in Fig. 8(c)] since it exhibits an interesting behavior. In close analogy with the KPO, the NP soft-mode excitation eigenfrequencies become imaginary upon crossing of the superradiant boundary, λ c . They, then, become real again when crossing to the light-blue region III at λ c2 . As discussed in Ref. [30], the transition from region I to III occurs alongside with a norm flip [cf. Eq. (3)] in the soft excitation eigenmode. This flip signals a particle-hole inversion in the excitation spectrum describing the system, marking the e-NP as an excited state of the system. As in the KPO, this feature survives the introduction of dissipation and, in the next section, we show that the e-NP becomes a dissipation-stabilized steady state. The SP excitation spectrum does not show any norm sign changes [30].

Open system
In a previous work [23], we analyzed the impact of adding single-photon losses to the Hamiltonian (36), cf. Eq. (7). The dissipator drastically changes the phase diagram. We here reiterate, mirroring Sec. II B, the main steps that led to the modified phase diagram to stress the close similarities between the KPO and the IDTC.
We start from the master equation, where L[a]ρ. We define the mean-field order parameters according to a = √ Nα, S i = Ni with i = X, Y, Z, and study the associated equations of motion [cf. Eq.(9)]: It is possible to find an analytical expression for the steadystate solutions [23], and in Fig. 9(a) we plot the resulting stationary order parameter behavior along the orange cut in Fig. 9(c). As for the open KPO [cf. Fig. 6], initially only the trivial NP solution is a stable steady state. Then, in regions II and III, up to four additional nontrivial superradiant solutions are allowed. Given the Z 2 symmetries of the problem, the nonzero solutions appear pairwise opposite of each other in phase space. We study the stability of the solutions against linear fluctuations [cf. Eq.(10)] and identify the stable and unstable solutions. In Fig. 9(b), we plot the eigenvalues associated with the possible steady states. The NP becomes unstable at λ c3 , where the superradiant transition takes place. Increasing the coupling further, the costability region III appears, where the low-population solution (e-NP) becomes stable again at λ c4 , and an unstable SP appears, SP un . Stable (unstable) solutions are highlighted as solid (dashed) lines in Fig. 9(a).
The open steady-state phase diagram is plotted in Fig. 9(c). Mapping one-to-one to the KPO, the open IDTC phase diagram includes the variety of regions, I, II, III, II , and III . Specifically, due to dissipation, the superradiant regions in Fig. 9(c) are now separated by a NP sliver (II') that connects two regions (I and III'), where the NPs are the only stable steady state. Stressing the KPO ↔ IDTC analogy further, dissipation stabilizes the NP in region III' of the parameter space, where it is an excited state of the closed system, this justifies the name e-NP [30].
In Figs. 9(d) and 9(e), we study the eigenvalues and variance of the IDTC NPs along a cut line that does not encompass the symmetry broken phase, brown cut in Fig. 9(c), cf. Figs. 6(d) and 6(e). The fluctuations' eigenvalues in region II' show the same overdamped behavior as the KPO (cf. Sec. II C, Fig 6) with the imaginary parts coalescing to zero, and the real parts splitting but remaining negative. This is accompanied by a finite variance squeezing along the tran-sition between the NP and the e-NP. Here, too, the exception pointlike discontinuity in the fluctuations' spectrum does not lead to sharp transitions in the variance. In conclusion, the open IDTC model exemplifies another (many-body) family of systems, where dissipative effects are enhanced by the interactions among the components of the system, leading to dissipation-induced PTs (region II') and the plethora of associated corollary signatures. In the following, we turn to highlight that similar dynamical signatures appear in this model as well.

Keldysh IDTC
We conclude our treatment of the IDTC model using the formalism introduced in Sec. II D to show that the two NPs are dynamically different. To find a path integral representation for the Hamiltonian Eq. (36), we bosonize the spin around the noninteracting ground state. The ground state is characterized by the magnetization S z = −N/2. The bosonized spin is obtained using the Holstein-Primakoff transformation, where we introduced the eight-component cavity-spin Nambu spinor, = [ c (ω), q (ω)], defined as the concatenated classical and quantum four-spinors: ) T with i = c, q, respectively. For the explicit form of the Green's functions, G R 4×4 , G A 4×4 and the Keldysh self-energy, D K 4×4 , in Nambu space, see Ref. [30].
Varying the action Eq. (42) with respect to the quantum components of the fields and substituting a c (t ) = √ 2α, b c (t ) = √ 2β, a q = 0, b q = 0, we obtain the coupled equations where α, β ∈ C. We can split them into their real and imaginary parts and recast them in terms of the mean-field order parameters α, X, Y, Z using the inverse Holstein-Primakoff transformation. This procedure yields the mean-field Eqs. (41).
We turn now to study the fluctuations on top of the mean-field solutions and address the spectral properties of the open IDTC. The different stationary phases of the IDTC can be described by the field a (b) given by the mean-field solutions α (β), and a fluctuating term, i.e., a → α + δa (b → β + δb) [cf. similar treatment of the KPO in Eqs. (31) and (32)]. The resulting IDTC quadratic fluctuation action is with the Green's functions and the Keldysh self-energy The definition of the coefficients δω 1 , δω 2 ,λ 1 ,λ 2 can be found in Ref. [30].
To obtain the photon-only action, we integrate out the Holstein-Primakoff fluctuation field δb and obtain a Keldysh functional integral only over the photon fields. The resulting photon-only action reads where A 4 (ω) = (δa c (ω), δa * c (−ω), δa q (ω), δa * q (−ω)) T is the photon four-vector that collects the classical and quantum field components. The inverse retarded Green's function of the photon is while the Keldysh component of the photon action is We now study the spectral properties of the NPs and of the SP.

Normal phase
In the NPs, we have as stationary mean-field solution, α = 0, X = Y = 0, Z = −1/2. Inverting Eq. (48) and using Eqs. (49) and (50), we obtain the photon Green's functions. The analytical expressions are quite lengthy and we report here only the retarded Green's function: We focus on the points λ 1 and λ 2 marked in Fig. 9(c), both of which are in the parameter space where the NPs are a stable steady state. The point λ 1 lies in region I below criticality, where the NP coincides with the ground state of the closed system. The point λ 2 lies instead in the critical region III , where the e-NP is stabilized by dissipation. Bottom: The fluorescence spectrum, S(ω), in (a) on top of the NPs, at points λ 1 , λ 2 , and in (b) on top of the SPs, at points λ 2 , λ 3 , see Fig. 9(c). The spectral response and power spectrum of the NPs feature a soft-mode peak inversion when going from the NP to the e-NP, λ 1 → λ 2 . This is a signature of the dynamically different nature of the two normal phases in the different regions of parameter space. In comparison, the SPs show no qualitative change in the spectra, neither upon crossing the superradiant threshold nor when considering the costability region.
In Fig. 10(a), we plot the spectral function Eq. (15), the power spectrum Eq. (16), and the fluorescence spectrum Eq. (18). Similar to the KPO case [cf. Fig. 7], we find a soft-mode peak inversion between the two different NP regions. The soft-mode peak swap marks the fact that the e-NP is an excited state of the closed system, which is dynamically distinct from the NP [30]. This is a final dynamical feature that highlights the universal phenomenology of dissipation-induced PTs and connects between the KPO and IDTC models. In both systems, coherent interactions facilitate that a weak dissipation channel stabilizes an excited state of the closed system and signatures of the dissipation-induced stabilization manifest in the dynamical response of the system.

Superradiant phase
In Fig. 10(b), we plot the dynamical responses on top of the superradiant solutions for points λ 2 , and λ 3 in parameter space, see Fig. 9(c). The former lies in the coexistence region, whereas the latter lies in the region where only the SP st is the stable steady state. The SP, as the PPS in the KPO, does not present qualitative changes in the spectra, marking once more the close analogy between the two a priori distinct physical systems.

IV. DISCUSSION AND OUTLOOK
At equilibrium, different phases and transitions between them are well understood using statistical physics arguments based on equilibrium distribution functions. Thus, a coupling to a thermal bath affects the distribution function, and at T = 0 (quantum) PTs occur only between ground states of the system. Moving to an out-of-equilibrium setting, the presence of drives and excitations leaking out into the environment requires us to break away from the framework of equilibrium distribution functions. Crucially, in such open systems, as we also highlight in this paper, a T = 0 environment does not imply a quantum PT in the resulting phase diagram of the system, i.e., the ground state does not necessarily correspond to the new stationary state of the system coupled to the environment.
Over the years, a plethora of methods have been applied to analyze out-of-equilibrium PTs, including diagonalization of Liouvillians [17] or non-Hermitian Hamiltonians (cf. third quantizations [13]), Keldysh action approach [11,12], and equations of motion [2] to name a few. Various observables were devised to characterize scenarios appearing in the outof-equilibrium, e.g., exceptional points marking stabilized regions due to gain and loss [15,28] and negative peaks in response functions as signatures of stabilized excited states [30,87]. In this paper, we collected such observations under a unifying framework and related them to simple underlying physical implications, using a mean-field analysis on top of which dynamical fluctuation and responses were explored.
We thus showed how dissipation can profoundly change the closed system physics by lifting the boundaries of the closed system phases. We, specifically, highlighted two scenarios that arise when dissipation overwhelms a symmetrybreaking term in the system: both cases involved the stabilization of a phase with higher energy than that of the ground state of the system, either corresponding to an excited physical state or to an unphysical state of the closed system. The latter exhibited an exceptional pointlike behavior, where the system is locked to the environment and its fluctuations become overdamped. On the other hand, when dissipation stabilizes a maxima of the closed system Ginzburg-Landau potential (i.e., an excited state), we observed an inversion in the dynamical response function, in accord with the inverted curvature of the potential. Crucially, the appearance of both types of scenarios is correlated to one another, where the overdamped region bridges between areas in phase space where the same state can be a ground or excited state.
We demonstrated the general phenomenology through two prominent example of open systems, namely, the KPO and the IDTC model. These models mark a wide family of experimental out-of-equilibrium domains, where nonlinearities, dissipation, and time-dependent drives interplay with one another [17,20,88,89]. Examples thereof range from nonlinear mechanical, electric, and photonic resonators, to driven light-matter systems in a similarly broad range of frequencies, including collective PTs in cold-atom experiments [90][91][92][93] and light-induced modifications of material properties [94,95]. As such, we believe that our work paves the way for a systematic study of dissipative stabilization of high-energy states, where key future directions can involve beyond mean-field analysis of the steady state and the introduction of more complex dissipation channels with their own symmetry-breaking tendencies, e.g., cat-state stabilization via two-photon losses and superradiant suppression via spin dephasing [20,96].
To study the stability of the various solutions, we split Eq. (27) into its real and imaginary parts and expand around the steady-state solutions as outlined in Sec. II B. Here we generally write α = α i for the solutions with i = 0, 1, 2, 3, 4. We obtain d dt to then find the steady-state fluctuation eigenvalues: Studying the real parts of ± KPO for the different α i and assuming U > 0, we identify the α 3,4 pair as the stable PPSs, and the α 1,2 one as the unstable PPSs, dubbed PPS st and PPS un in the main text, respectively. Note that for U < 0, the α 3,4 pair is unstable and the α 1,2 one is stable.
The variance of fluctuations is calculated from the correlation functions [cf. Eq. (11)] as outlined in Sec. II B. We solve the resulting coupled equations for the NP solution in the linear regime (U ≈ 0) and obtain Var(α Re ) = Var(α Im ) = Crucially, the result discontinuous only at the phase boundary = ± √ G 2 Re + G 2 Im − κ 2 . Thus, the abrupt splitting in the fluctuation eigenvalues (exceptional pointlike scenario) does not lead to any abrupt changes in the variance.