Liouvillian spectral collapse in the Scully-Lamb laser model

Phase transitions of thermal systems and the laser threshold were first connected more than forty years ago. Despite the nonequilibrium nature of the laser, the Landau theory of thermal phase transitions, applied directly to the Scully-Lamb laser model (SLLM), indicates that the laser threshold is a second-order phase transition, associated with a $U(1)$ spontaneous symmetry breaking (SSB). To capture the genuine nonequilibrium phase transition of the SLLM (i.e., a single-mode laser without a saturable absorber), here we employ a quantum theory of dissipative phase transitions. Our results confirm that the $U(1)$ SSB can occur at the lasing threshold but, in contrast to the Landau theory and semiclassical approximation, they signal that the SLLM"fundamental"transition is a different phenomenon, which we call Liouvillian spectral collapse; that is, the emergence of diabolic points of infinite degeneracy. By considering a generalized SLLM with additional dephasing, we witness a second-order phase transition, with a Liouvillian spectral collapse, but in the absence of symmetry breaking. Most surprisingly, the phase transition corresponds to the emergence of dynamical multistability even without SSB. Normally, bistability is suppressed by quantum fluctuations, while in this case, the very presence of quantum fluctuations enables bistability. This rather anomalous bistability, characterizing the truly dissipative and quantum origin of lasing, can be an experimental signature of our predictions, and we show that it is associated with an emergent dynamical hysteresis.

emission and absorption coefficients [1], paved the way to the creation of devices and lasers exploiting the stimulated emission of optical radiation [2]. The Scully-Lamb laser model (SSLM) is a remarkably predictive quantum description of a laser oscillator [3]. Recently, the SLLM has attracted much interest in different fields, ranging from open quantum systems phenomena (e.g., exceptional points [4,5] and quantum thermodynamics [6,7]) to more exotic ones (the characterization of entropy in Bose-Einstein condensates and even in black holes [7]). As pointed out in, e.g., Ref. [8], the laser population inversion can never occur in thermal systems, since it would imply the presence of negative temperatures [9]. As such, the laser is an out-of-equilibrium system [10]. Despite this fact, analogies between thermal phase transitions and lasing nonequilibrium transitions were drawn more than forty years ago [11,12] using the Landau theory of phase transitions [8,[13][14][15]. Accordingly, the SLLM threshold is a second-order phase transition associated with the spontaneous breaking of the U (1) symmetry [12], and these results are supported by a semiclassical analysis.
Notwithstanding the success of the Landau and semiclassical models in providing a qualitative description of lasing, still they lack a rigorous description of its transition. Indeed, the semiclassical description includes dissipative processes but neglects quantum fluctuations. The Landau theory takes into account thermal-like fluctuations but, by assuming thermodynamic equilibrium, neglects the out-of-equilibrium fluctuations [16,17].
As such, it is natural to ask if there is some unveiled physics behind the lasing process resulting from a genuine theory of quantum dissipative phase transitions, i.e., based on the quantum Liouvillian framework [18]. Here, using the Lindblad master equation and its associated Liouvillian superoperator [18], we analyze the SLLM explicitly taking into consideration the U (1) symmetry of the model. We divide the Liouvillian into its symmetry sectors, i.e., we separate the dynamics into families of states which evolve independently one from another due to the presence of the U (1) symmetry [19][20][21]. We demonstrate that the SLLM second-order phase transition is characterized by a Liouvillian spectral collapse. That is, a high degeneracy of the Liouvillian spectrum emerges as diabolic points of infinite order, one for each symmetry sectors, triggering dynamical hysteresis and other critical properties which cannot be explained by the U (1) SSB alone. We borrow the term spectral collapse from the two-photon Rabi model-describing a parametric exchange of excitations between a bosonic field and an ensemble of spins-where high-order degeneracies of the Hamiltonian emerge at the critical coupling (see the discussion in, e.g., Refs. [22,23]). Although similar, this Liouvillian spectral collapse is not a trivial extension of its Hamiltonian counterpart, but rather a novel type of criticality. Using the theory recently developed in Ref. [24], we give a more insightful description of the nature of the spectral collapse by preventing spontaneous symmetry breaking (SSB) from the SSLM. To do that, we consider a generalized SLLM, i.e., the standard SLLM model with constant additional dephasing. While the U (1) symmetry is maintained, long-lived phase coherences are destroyed, thus preventing the emergence of a symmetrybroken phase and highlighting some properties which distinctively characterize the spectral collapse.
It is worth mentioning that, according to (semi)classical equations of motion, a bifurcation phenomenon in equilibrium systems can also be accompanied by a second-(third-) order phase transition with (without) symmetry breaking [25][26][27]. Moreover, in thermal systems, a second-order phase transition without symmetry breaking may also occur, as can be described in a topological-transition framework [28]. However, the SSLM with the removed SSB always has a unique steady state, as shown in Ref. [24]. This property prohibits the occurrence of semiclassical bifurcations and multistabilities, which, in turn, demonstrates that its nature is profoundly different from that of the known effects.
Notably, some pioneering works grasped hints of the emergence of a critical timescale in the laser not associated with SSB. Indeed, in Ref. [29] the presence of a slow timescale was argued using a potential-well approximation to compute the tunneling time in a Fokker-Planck formalism. Moreover, in Ref. [30], it is explicitly shown that there is a closure of the spectral gap (to be associated with the emergence of a critical timescale in the photon-number evolution) in the limit of a vanishing saturation rate. These results corroborate the validity of the Liouvillian analysis of the SLLM transition.
An analogous effect of a mixed-order phase transition was described in Ref. [31], where it was shown that a spin model, also characterized by a U (1) symmetry, can undergo a dissipative phase transition without SSB. Furthermore, a Liouvillian spectral collapse, similar to the one studied here, has also been observed in more exotic lasing models, such as the squeezed laser in Ref. [32], or in a U (1) symmetric superradiant model [21]. In this regard, the model under consideration exhibits behavior similar to that in these other U (1)-symmetric systems. With the help of our analytical description via the P -function, we can precisely estimate the nature of quantum fluctuations, which induce this peculiar character of the transition. Note also that the SSB in our model can be arbitrarily removed according to the theory developed in Ref. [24]. This approach enables us to describe more deeply the interplay between the U (1) symmetry and quantum fluctuations in triggering critical timescales.
The article is organized as follows. In Sec. II, we recall the SLLM, and we introduce its generalization in the Liouvillian description. We recall the results of the semiclassical theory, and we derive the quantum solution in the weak-gain-saturation regime using the Glauber-Sudarshan P -representation. We discuss the connection between dephasing and the convergence rate to the steady state, and connect the P -function to a thermodynamic potential within the Landau theory. In Sec III, we demonstrate the one of the main results of this article, namely the presence of a Liouvillian spectral collapse in the SLLM. To do that, we consider the standard SLLM (i.e., without additional dephasing). Exploiting the Liouvillian symmetries, we confirm the presence of SSB within the Liouvillian theory. Nevertheless, we show the presence of another critical phenomenon, not captured by the Landau theory, i.e., a Liouvillian spectral collapse at the critical point. In Sec. IV we investigate the generalized SLLM, i.e., the standard SLLM in the presence of additional dephasing. As proved in Ref. [24], this leads to a second-order phase transition in the absence of SSB. Notably, both the semiclassical theory and the Landau approach fail to describe this model. Hence, we validate the necessity to consider a Liouvillian theory to properly characterize the lasing transition. We discuss the emergence of multistability at the critical point, witnessed by a dynamical hysteresis. In Sec. IV A, we detail the novelty of the spectral collapse, as well as similarities and differences with respect to other phase transitions. We characterize the nature of the fluctuations at the critical point, showing that the critical slowing down of the spectral collapse can be interpreted as the time required for a random process with constant diffusion to cover a larger and larger system size. In Sec. V, we discuss the results obtained and provide conclusions and perspectives. Appendix A details the procedure to obtain the Liouvillian form of the SLLM equation of motion, as well as its limits of validity. Moreover, it provides additional details on the properties of the SLLM Liouvillian. Finally, in Appendix B, we provide additional details concerning the lack of SSB in the SLLM using quantum trajectories.

II. THE SCULLY-LAMB MODEL AND ITS GENERALIZATION
The time evolution of a quantum system under the Born and Markov approximation [16] is captured by the Lindblad master equation ( = 1): whereρ(t) is the system reduced density matrix at time t, H is the Hamiltonian describing the coherent part of the system evolution, L is the Liouvillian superoperator [57], and D[L j ] are the so-called Lindblad dissipators, whose action is The operatorsL j are the jump operators, which describe how the environment acts on the system inducing loss and gain of particles, energy, and information. A central role in the following discussion is played by the steady stateρ ss , i.e., the state which does not evolve any more under the action of the Liouvillian: We indicate the expectation values of operators in the steady state as ô ss = Tr[ρ ssô ].

A. The model
We consider here a quantum description of the Scully-Lamb laser model, describing a photonic mode pumped by the injection of inverted atoms into the laser cavity [58]. By tracing out the atomic degrees of freedom [59], the system Hamiltonian becomes that of an optical cavity (a harmonic oscillator) whereâ (â † ) is the bosonic annihilation (creation) operator. The evolution of the photonic field is captured by the three jump operators (see the discussion in Appendix A and Refs. [5,30,[59][60][61][62][63][64]): whereL 1 describes the laser gain,L 2 captures the field dephasing, andL 3 represents the particle loss. The jump operators are characterized by rates: A for the unsaturated gain of the active medium [10] (we are considering here a laser without an absorber), B for the gain saturation, Γ for the dissipation rate, which corresponds to the inverse of the photon lifetime. The rate β represents the overall dephasing rate of the system; accordingly, η represents an additional dephasing rate beyond the rate induced by the gain saturation. If η = 0 (η = 0) the system is a standard (generalized) SLLM.
The Lindblad form of the SLLM master equation is well-defined in the weak-gain saturation (WGS) regime, valid under the condition The WGS conditions are satisfied around the lasing threshold [59]. Far away from the threshold, the jump operatorL 1 in Eq. (5) can be nonphysical, and a different master equation needs to be employed [59,62]. Since the main aim of this paper is to investigate the nonequilibrium phase transitions of the laser at the threshold, we can safely useL 1, 2, 3 in Eq. (5) (see also Appendix A 2).

B. The lasing transition within the semiclassical approximation
As discussed in Ref. [56], for this U (1) model we can set ω = 0 via a reference frame change: from the laboratory one to the one rotating at the cavity frequency. This leaves the overall properties of the system unchanged. Within such a rotating frame, the semiclassical approximation (often-but not always-justified in the limit of intense fields) amounts to consider that the correlation functions â † mân = â † m â n → α * m α n . Equivalently, the density matrixρ(t) of the system is, at all times, in the coherent stateρ(t) = |α α|. Under this approximation, and working in the WGS regime of Eq. (6), one can compute the mean-photon number in the steady state [60,64], obtaining either a trivial solution â †â ss = n ss = 0 or Stability analysis reveals that, for A > Γ+η, the solution n ss = 0 becomes unstable, and the system passes nonan-alytically from zero to a nonzero number of photons. Thus, the system undergoes a second-order phase transition, whose order parameter α = â ss = 0 indicates the U (1) SSB of the model (see the discussion in Sec. III B). Let us remark that the semiclassical equations of motion for the fields obtained from Eq. (1) do not include any contribution from quantum fluctuations. We also note that a rigorous quantum solution for the steady state reveals that at the threshold â †â ss = 0 [59,62].

C. The full quantum solution
Let us analytically proof that the semiclassical analysis is wrong in the case η = 0, thus demonstrating the necessity to go beyond the semiclassical approximation to correctly describe the lasing transition of the SLLM.
The Glauber-Sudarshan P -representation is a mapping of the density matrix from the bosonic Hilbert space to the complex phase space C via a real-valued function P (α), which readŝ where |α is a coherent state [58,65,66]. In this formalism, the matrix equation forρ(t) is transformed into a differential equation for the quasiprobability function P (α, t) [67,68] by repeatedly applying the rules: Under the weak-gain-saturation condition in Eq. (6), Eq. (5) becomes In the steady state ∂ t P ss (α) = 0, the P function depends only on the field intensity |α| 2 . By rewriting Eq. (10) in the polar coordinates (r, φ), one can separate P (α, t) into its radial [P (r, t)] and angular [P (φ, t)] components. This is also true for the steady state, and P ss (r, φ) = P ss (r)P ss (φ). At the steady state, P ss (φ) = 1/(2π), while By solving Eq. (11), one obtains P ss (α) = P ss (φ)P ss (r) = where N is a normalization constant.
In striking contrast to the semiclassical approach, the parameter η does not appear in the steady-state solution (7). As such, the evident contradiction between the photon number obtained using the semiclassical approximation in Eq. (7) (depending on η) and the P -representation in the weak-gain-saturation regime (which does not de-pend on η) cannot be reconciled, even in the thermodynamic limit. In other words, the effect of the photon dephasing is overrated in determining the intensity of the field in the semiclassical approach, and a (superposition of) coherent state(s), with the amplitude in Eq. (7), can never be the steady state of the system for η = 0. We provide a detailed description of this semiclassical approximation failure in Appendix B.

Diffusion coefficient and the loss of coherence
Although η plays no role in the determination of the steady state, it strongly affects the dynamics of the system. To prove the last affirmation, we notice that the loss of coherence is encoded in the evolution of P (φ, t).
We can estimate such time by considering an initial state P (r, φ, t = 0) = P ss (r)δ(φ − φ 0 ). Accordingly, we have [c.f. the right-hand side of the Eq. (10)]: This is just a diffusion equation with a diffusion coefficient The diffusion coefficient D determines the spectral line of a laser. In the standard Scully-Lamb laser theory, the parameter η = 0, B/A 1, and the laser line tends to zero with increasing photon number A/ n ss → 0. Thus, as expected, a state will maintain its coherence even for infinitely large times t → ∞.
On the contrary, in the generalized SLLM, the laser linewidth remains finite, since D = β/2 η/2 in the large-gain limit. Thus the P function already hints at the fact that the system cannot break the U (1) symmetry, since it cannot retain a phase. That is, the generalized SLLM, in the presence of η = 0, cannot be characterized by a U (1) SSB.
We can already see the profound and nontrivial interplay between the dynamics of the system, the symmetry breaking, and the emergence of a phase transition.

D. The Landau theory
One can try to go beyond the semiclassical analysis of the phase transition by using the well-known Landau theory of phase transitions, assuming an analogy with thermal phase transitions in equilibrium systems. This approach seems particularly effective in the standard SLLM (for which η → 0). Indeed, the exponent of P ss in Eq. (12) can be directly associated with the "free energy" G of the system [12]. Remarkably, the relatioṅ E ≡ − ∂G ∂E can be derived, where E ≡ α ss + α * ss is the electric field obtained via the semiclassical approximation in Eq. (7). From the Landau theory one deduces the presence of a second-order phase transition, signalled by a nonanalytical change in the intensity of the electric field (i.e., the photon number). This is associated with the U (1) SSB in the system. Within the generalized SSLM (η = 0), no conclusions about such a correspondence between the P -function and the corresponding G can be made. Even if blindly assuming that the exponent of P ss still accounts for the free energy, one would erroneously conclude that still the phase transition, associated with the SSB of U (1), is taking place [12].

III. DISSIPATIVE PHASE TRANSITION OF THE SCULLY-LAMB MODEL
A dissipative phase transition is a discontinuous change in the steady stateρ ss as a function of a single parameter [18,33]. In the following analysis, the gain A in Eq. (5) plays such a role (see, e.g., Fig. 1), and the definition for the second-order dissipative phase transition reads where A c is the critical point. This nonanalyticity occurs in the thermodynamic limit. However, the effects of an emerging criticality can also be witnessed in finite-size systems. Indeed, there is a profound connection between the dynamical properties of a dissipative system and the emergence of a phase transition, as proved in Refs. [18,33] and experimentally demonstrated in Refs. [41,45,69]. Criticality is accompanied by the emergence of the so-called critical slowing down, i.e., the appearance of infinitely long timescales in the system dynamics.
To properly characterize the emergence of a critical timescale, we resort to the Liouvillian spectrum. Mathematically, we can define the eigenvalues λ i (representing the time scale of the problem) and eigenvectorsρ i (encoding the states explored along the dynamics) of L via While Re(λ i ) indicates a decay time towards the steady state, Im(λ i ) encodes the frequency of the associated oscillations. The Liouvillian is a superoperator, which acts on operators to generate new operators, in the same way in which an operator acts on a vector to generate a new vector. We use the symbol • to indicate the density matrix placeholder. For example, Since the Liouvillian is a linear superoperator [i.e., L(aÂ+bB) = aLÂ+bLB], L can be represented in a matrix form. A detailed discussion of superoperator properties can be found in, e.g., Refs [17,57,70].

A. The thermodynamic limit
Phase transitions and nonanaliticity can only emerge in the thermodynamic limit. Thus, we need to introduce a well-defined thermodynamic limit to properly compare the Liouvillian results with those of the Landau theory and the semiclassical approach discussed above. In a lattice system with L sites, this corresponds to increasing the size L to infinity. In bosonic systems, one can exploit the infinite dimension of the Hilbert space to observe a nonanalytical change in the steady state, as discussed, e.g., in Refs. [18,47,56,[71][72][73]. To do so, one considers an appropriate rescaling of the parameters, so that the size of the Hilbert space increases, but the rescaled observables merge far from the critical point. We introduce an effective scaling parameter N transforming the photon number n → N n = N (A − Γ)/B in Eq. (7). Hence, any transformation is valid for any real number µ, and the thermodynamic limit is reached by increasing N . While the steady state is unchanged by different choices of µ, the choice µ = 0 ensures that the dissipation rate Γ is kept constant (providing a natural timescale) for the problem. Considering the laser model of Ref. [59] and Appendix A 1, where the gain is caused by the injection of inverted atoms, the scaling parameter N in Eq. (17) corresponds to decreasing the light-matter coupling g, while increasing the pumping rate (i.e., without increasing the gain A). Thus, the parameter N can be seen as the number of constantly pumped atoms required to compensate the decreasing light-matter interaction. The thermodynamic limit N → ∞ is well defined, since we are increasing the number of particles, while keeping constant the interaction energy density.
In Fig. 1(a), we plot the photon number â †â /N , obtained by numerically solving Lρ ss = 0, for different values of N . The photon number becomes sharper and sharper with increasing N , converging towards the semiclassical result of Eq. (7). This confirms the presence of criticality in the thermodynamic limit, and corroborates the analysis of the scaling parameter N in Eq. (17).

B. Symmetry breaking
We can refine the spectral analysis in Eq. (16) by introducing the symmetries of the Liouvillian. The SLLM is characterized by a U (1) symmetry of the Lindblad master equation, because the transformationâ →â · exp(iφ) leaves the Liouvillian unchanged for any real number φ. Thus, the superoperator  quantity. This is not always the case for Liouvillian symmetries [19,20,74]. However, all eigenmatrices of L must be eigenmatrices of U [18,19,74,75]; that is, where u k is one of the eigenvalues of U, and k is an integer (see Appendix A 3 for more details). We refer to the span of all the eigenmatrices characterized by the same u k as a symmetry sector. Each symmetry sector is a part of the Liouvillian space which is not connected to its other parts (sectors) by the Liouvillian dynamics [21,56,76]. In other words, the Liouvillian can be decomposed into a direct sum of superoperators L k acting on different symmetry sectors as As such, we can relabel the eigenvalues and eigenmatrices 7 as λ i , respectively, where k tracks the symmetry sector: We order the eigenspectrum of each sector according to Using this convention, the steady state is proportional toρ We define the Liouvillian gap λ as the nonzero eigenvalue of the Liouvillian whose real part is the closest to zero, i.e., The inverse of λ indicates the slowest relaxation time of the system. The closure of the gap (λ → 0) indicates a diverging timescale and, thus, criticality.
Within this representation, we can write any initial density matrix asρ for an appropriate choice of coefficients c j , and, thus, we haveρ 1. Spontaneous breaking of the U (1) symmetry In a Liouvillian framework, the breaking of a Z N symmetry takes place when (N − 1) eigenvalues, each belonging to a different symmetry sector, become zero: λ [18]). These eigenvalues must become zero not only at the critical point, but also in the entire broken-symmetry regime, ensuring that there exist several steady states which are not eigenstates of the symmetry operator. Since U = lim N →∞ Z N , a SSB in the SLLM means that in each symmetry sector (except u 0 ) there is a vanishing eigenvalue. As suggested by the semiclassical analysis, we should observe this fea- In Fig. 1(b), we demonstrate that, converging towards the thermodynamic limit N → ∞, the Liouvillian gap [defined in Eq. (23)] tends to zero for A > A c . This indicates a symmetry breaking, as also confirmed for similar models in Refs. [32,50]. To prove that the SSB is that of U (1), in Fig. 1(c) we plot λ (k) 0 for k = 1, 2, 3 for fixed N = 100. In each sector we observe a similar "closure" of the gap, confirming the presence of U (1) symmetry breaking. We verified that in all sectors (up to the numerical cutoff) the closure of the Liouvillian gap becomes more evident as we increase N .
Technically speaking, the SLLM lasing transition is a time crystal. The imaginary part of the plotted eigenvalues is never zero (Im[λ (k) 0 ] = kω), and therefore undamped oscillations take place. However, as proved in Ref. [56], the imaginary part of λ (k) 0 can be set to zero via an appropriate change of the reference frame, resulting in a genuine symmetry-broken phase.
C. Spectral collapse in the standard Scully-Lamb laser model For this standard SSLM (η = 0) the results of the Landau theory and those of the semiclassical approximation agree with the Liouvillian theory in predicting the U (1) SSB [12]. However, we just analyzed λ (k) 0 , i.e., the slowest relaxation rate in each symmetry sector, ignoring the remainder of the spectrum. By analyzing λ (k) j for j > 0, we show that the SLLM transition is characterized by a different kind of criticality, which we call a Liouvillian spectral collapse. For this standard SLLM, this does not mean that the SSB does not take place, but rather that the lasing threshold has a richer structure.
In Fig. 2(a), we plot the smallest nonzero eigenvalue in the sector u 0 , namely, λ describes the slowest relaxation rate of all those observables of the form (â † ) mâm . Surprisingly, we observe that λ (0) 1 becomes zero at the critical point. Note that the sector u 0 comprises the Liouvillian eigenmatrices which are phase-independent (i.e., diagonal such as the steady state). Consequently, the SSB, which is associated with the development of a well-defined phase, is not directly associated with the criticality in u 0 . The absence of a critical timescale in u 0 has been observed in models with a Z 2 SSB, as can be inferred in, e.g., Refs. [18,52]. We stress that, to numerically distinguish λ (0) 1 from all the other eigenvalues, which become zero at the transition, one needs to properly divide the system in its symmetry sectors, as detailed in Sec. III B.
Remarkably, this point-like closure of the gap occurs for all the symmetry sectors for the same critical value of A c = Γ [we plot those of λ . We numerically tested that, by increasing N , more and more eigenvalues collapse, in their real part, towards zero.
One might think that a possible (although involved) explanation of this phenomenon is that the real part of each eigenvalue converges to zero independently of the others, so that some eigenvalues are "more critical" than the others. However, this conjecture is false. In Fig. 3, we plot the minimum of the eigenvalues λ   Fig. 2 we demonstrated that the position of the minimum tends to the critical point Ac = Γ, here we can clearly see that all the eigenvalues have the same convergence rate, further demonstrating the existence of a spectral collapse. Parameters are the same as in Fig. 1.
transition for those eigenvalues which present a spectral collapse (i.e., λ (k) j for j > 0). Obviously, the convergence rate of all these eigenvalues is identical, meaning that they are all resulting from the same phenomenon.
In other words, we demonstrated that the lasing threshold according to the SLLM corresponds to an infinite degeneracy of nondecaying states at the critical point. Mathematically speaking, there is a diabolic point of an infinite order in each of the infinitely many symmetry sectors of the Liouvillian. We are sure that the criticality is not associated with exceptional points since the degeneracies emerging in the thermodynamic limit, where Re[λ (k) j ] = 0 prevents the occurrence of exceptional points [70].
We refer to the emergence of infinitely degenerate manifolds as a Liouvillian spectral collapse, in analogy to the definition applied to the two-photon Rabi model [22]. Notice that, while the real part of λ (k) j → 0, the imaginary part is always fixed by the symmetry sector, since Im λ Let us remark that this degeneracy does not imply that all the processes are infinitely long lived. Indeed, there is a part of the spectrum which never collapses to zero. For instance, at the critical point, a state initialized in a Fock state |n n| rapidly changes its parity, defined bŷ P = exp (iπâ †â ), on a timescale which is always of the order of 1/(Γn).
This nonanalytical change occurs at the same point as the U (1) SSB. As we discuss in more detail in the next section, the spectral collapse is a fundamental ingredient of criticality, and the SSB is rather a byproduct of this process.

IV. WITNESSING AND CHARACTERIZING THE LIOUVILLIAN SPECTRAL COLLAPSE
As indicated by the P -representation in Sec. II C, the steady state is unaffected by η. On the contrary, the prediction of the semiclassical analysis in Eq. (7) is that the phase transition point is shifted by η. Indeed, η cannot affect the dynamics of the steady-state symmetry sector, as it was recently demonstrated in Ref. [24]. For this reason, by keeping η constant even in the thermodynamic limit, it is possible to witness a second-order phase transition, even if we prevent SSB.
In Fig. 4(a), we show the validity of Eq. (12) and confirm the validity of the results in Ref. [24] because, despite an additional dephasing, the photon number is identical to that of the SLLM for η = 0 [c.f. Fig. 1(a)]. Interestingly, the semiclassical coherent-state approach completely misses the onset of criticality. A detailed discussion of why the semiclassical approximation fails is given in Appendix B. We also stress that a similar "wrong" prediction of SSB by a semiclassical analysis has been pointed out in Ref. [31] for a spin model, where the quantum solution shows a second-order dissipative phase transition without SSB.
As discussed in Sec. II D, a thermodynamic free energy, based on the P -function, would instead predict the pres-ence of SSB. We conclude that such a prediction is wrong, and we cannot use the P -function as a thermodynamic potential within the Landau theory.
Notice that, in agreement with Ref. [24], and as one can see from Fig. 4(b), the nonzero η precludes any SSB (after the transition, there are no multiple steady states). However, the Liouvillian gap shows all the signs of the gap closure associated with the spectral collapse, as shown in Fig. 2. Indeed, the spectral collapse takes place in the u 0 symmetry sector, i.e., those states which are unaffected by η, because they are symmetric under U (1) action (see Appendix A 3). For this reason, using the generalized SSLM, we can distinguish the effects induced by the spectral collapse from those depending on the SSB.
We also confirm the predictions of Eq. (14). Indeed, after the phase transition occurs, η prevents the emergence of a symmetry-broken phase, and the Liouvillian gap for a large-enough N is bounded by η/2. We conclude that the dephasing destroys the coherence of a laser state and no SSB is taking place because the retaining of a phase is the characteristic of the U (1) SSB. Further numerical simulations confirm that, for all values of A, an initial state with a nonzero phase rapidly looses its phase (on a timescale of the order of η/2).

A. Anomalous multistability and dynamical hysteresis
The presence of the Liouvillian spectral collapse and the degeneracy of eigenvalues is not just a theoretical artefact, because the existence of nondecaying processes at the critical point has profound consequences in the physics of the SLLM which can be experimentally observed.
The presence of symmetry sectors of the Liouvillian affects also the dynamics of certain classes of operators. For instance, the evolution ofâ †â depends only (explicitly or implicitly) on that of â † mâ m . This correspondence between operators and density matrices becomes particularly useful when using the spectral decomposition in Eq. (24). For instance, since U(â †â ) =â †â one can easily prove that That is, only the eigenmatrices of the symmetry sector u 0 have nonzero photon number expectation value. Therefore, the spectral collapse can be revealed by the photon-number evolution, since one is monitoring only the u 0 sector, i.e., the symmetric one which cannot experience criticality due to the U (1) SSB [see Appendix A 3 and Fig. 9(a)]. At the spectral collapse, the mean photon number â †â (t) must experience a critical slowing down [42,79,80]. Although the spectral collapse can be witnessed only for N → ∞, a dynamical hysteresis characterizes finite-size systems and can reveal an emerging criticality [41].  To witness a dynamical hysteresis, one can slowly (with respect to Γ) change the value of A [i.e., ∂ t A(t) Γ]. Far from the critical point, we expect the evolution to be adiabatic, since λ We analyze the emergent critical slowing down and the dynamical hysteresis of the SLLM lasing transition in Fig. 5. In Figs. 5(a) and 5(b), we initialize the system to be close to the vacuum, and we observe its evolution at the critical point. Figure 5(a) shows that the time needed to reach the steady-state photon number is much longer that Γ, which represents the timescale far from the critical point in the unbroken symmetry phase. Figure 5(b) shows that, apart from the initial transient dynamics, the decay towards the steady-state is purely exponential. An exponential fit of the curve retrieves λ In Fig. 5(c), we show the presence of the dynamical hysteresis for N = 100. The red curve represents the master equation evolution for a linear ramping of the parameter A, reading for t f = 200Γ. Similarly, the blue curve represents the inverse descending process for the same t f = 200Γ. The thin superimposed curves represent the results of a single quantum trajectory (i.e., reproducing the result of a single ideal experiment [81], see also Sec. IV B). Far from the critical point, i.e. A Γ or A Γ, the evolution is adiabatic, and the two curves overlap with that for the steady state (black dotted curve). Instead, around the critical value, the two dynamics are significantly separated, proving the presence of a dynamical hysteresis.
In Figs. 5(d) and 5(e), we show the results of the same simulation but for N = 1 and N = 10. While for N = 1 there are no signs of the hysteresis, for N = 10 we observe some minimal differences between the ascending and descending curves around the critical regime. Comparing Figs. 5(c), 5(d), and 5(e) we conclude that the hysteresis becomes more and more pronounced with increasing N , confirming the results of Fig. 2. Hence, Figs. 2(a) and 2(b), describe the emergence of a nonsymmetry breaking critical timescale.
Importantly, these results are independent of the additional dephasing in Eq. (5) (all the plots in Fig. 5 perfectly overlap). Indeed, one can demonstrate that the u 0 symmetry sector is unaffected by the value of β. Such a dynamical hysteresis, together with the nonanaliticity of the photon number in the cavity, could constitute an experimentally realizable test of the spectral collapse and of the presence of a second-order phase transition with or without U (1) symmetry breaking.
We also stress that the presence of a dynamical hysteresis is not an artefact of the Liouvillian formulation of the SLLM, but similar results, based on the original model introduced by Scully and Lamb in Ref. [3], have been obtained in Ref. [30]. Nevertheless, using the Liouvillian framework, we can assign a clear physical meaning to these phenomena.

B. Quantum fluctuations at the spectral collapse
To understand the cause of multistability at the spectral collapse, let us show that correlations play a fundamental role in the SLLM criticality. First, let us consider the second-order single-time correlation function defined as g (2) We plot the results for g (2) (0) in Fig. 6(a). We see that, at the critical point, where the phase transition takes place, the value of g (2) (0) abruptly passes from 2 to 1.
Remarkably, at the critical point, g (2) (0) = 1, indicating that the system is not a coherent state.

P -Representation
The role of fluctuations in determining the emergent slow timescale can be proved by analyzing the Fokker-Planck equation in Eq. (10). As we previously discussed and demonstrated, all the Liouvillian eigenmatrices belonging to the u 0 sector (i.e., those which are invariant under any rotation, c.f. Appendix A 3) are unaffected by the parameter β. Thus, to investigate the dynamics of any initial state which is only a combination of matrices belonging to the u 0 sector, we can set β = 0. The Fokker-Planck equation can be separated into radial and spatial components, i.e., P (α, t) = P (θ, t)P (r, t), and P (θ, t) = P ss (θ) = 1/2π. As such, the P -function for the field intensity r acquires a similar form to that in the Eq. (11) and In this case, the Fokker-Planck equation is an Ornstein-Uhlenbeck process, where the drift and diffusion terms correspond to the first and second terms in the right-hand side of Eq. (31), respectively. Moreover, the drift is not constant, but depends on the field intensity, i.e., the drift sign flips depending on whether the current intensity is larger or smaller than its mean value n ss in the process, and which allows the system to evolve to its equilibrium.
We can identify three regimes: (1) well below the critical point (A < Γ), the mean photon number is small, thus the field intensity (photon number) just oscillates near its mean n ss ≈ 0 [see also Fig. 7(a)  coefficient A is small, resulting in a decreasing deviation from the mean with a decreasing gain.
(2) At the critical point (A = Γ), the mean photon number diverges in the thermodynamic limit n ss → ∞. However, the rescaled photon number vanishes because n ss /N → 0. We conclude that also the drift coefficient vanishes because r(A − Γ − Br 2 /N ) → 0 for sufficiently small r in the thermodynamic limit. Consequently, possible fluctuations of the photon number around their mean induced by the drift are suppressed. The typical timescale of the system is, thus, determined solely by the diffusion coefficient (proportional to the gain A). Notably, the diffusion remains constant irrespective of the thermodynamic parameter N . As such, the typical timescale, which characterizes a system near the origin, can be seen as the time required from the system to explore the entire available configuration space, which becomes divergently large in the thermodynamic limit. For the finite-size systems, this diffusion-induced slowing down can be observed by trajectories initialized near the origin [as confirmed by Fig. 7(b)].
(3) Finally, above the threshold (A > Γ), the drift becomes large again. Hence, the drift quickly drags the photon number to the proximity of its mean value, around which it remains to oscillate in the long-time limit [see Fig. 7(c)]. The size of these fluctuations, however, becomes small with respect to the mean photon number.

Quantum trajectories
To visualize the effects of fluctuations in the full quantum simulation and a "chaotic-like" nature of the system at the critical point, we resort to quantum trajectories. Quantum trajectories representρ(t) as the average over several stochastic evolutions of a pure state wave function. An in-depth discussion of the quantum trajectory formalism goes beyond the purpose of this article, and we refer the interested reader to, e.g., Refs. [81,82]. Suppose that we can monitor each time a quantum jumpL 1 , L 2 , orL 3 takes place. Between two quantum jumps, the infinitesimal evolution of the wave function is dictated by a non-Hermitian Hamiltonian [77,83], and This nonunitary evolution is interrupted by a random quantum jump of one of the operatorsL j , each one occurring at each interval (t, t + dt) with probability The occurrence of a quantum jump can be seen as a "click" of a perfect detector [84]. In this case, the wave function instantaneously changes as where dt is an infinitesimal time.
In Fig. 7 we plot the photon number along single quantum trajectories for different values of the gain A. Before the critical point, A < Γ in Fig. 7 (a), the signal is noisy and randomly passes between the vacuum and states with several photons. Despite this fact, we see that the quantum trajectories change rapidly reaches its steady state value, around which it oscillates. At the critical point, A = Γ in Fig. 7 (b), we still see a very noisy signal, but in all the trajectories a long time is required before the system reaches states with a high number of photons. Finally, in Fig. 7 (c) where A > Γ, we observe a very different behavior: the system rapidly reaches a large photon number, and then oscillates around this value.
The intuition which we gain from this picture is that the critical point, similarly to the region A < Γ, is characterized by large oscillations. However, on a short timescale, the size of this oscillation is quite small [c.f. Fig. 7 (b)]. This implicates that the critical slowing down associated with the spectral collapse is due to both the chaotic nature of light (requiring to have large oscillations) and the fact that the size of the oscillations is comparatively small. Thus, it requires a large amount of time for any quantum trajectory to explore the entire space of accessible configurations. This picture confirms our predictions based on the P -representation analysis.
In this article, we analyzed the SLLM transition in the WGS regime, within the open quantum system framework provided by the spectral properties of the corresponding Liouvillian superoperator. Our analysis takes into account the nonequilibrium character of lasing and its quantum fluctuations, making this analysis more rigorous than previous studies, where the Landau theory and semiclassical analysis have been used. Within the Liouvillian formalism, we demonstrated that the lasing transition in the SLLM is associated with a Liouvillian spectral collapse.
The spectral collapse indicates that, in one or more symmetry sectors [i.e., the Liouvillian eigenspaces whose dynamics is not interdependent due to the presence of the U (1) symmetry], infinitely many eigenvalues retain their imaginary part, but collapse to zero in their real part. As such, in the presence of symmetry breaking, the spectral collapse can be seen as the presence of infinitely many diabolical points [one for each symmetry sector of the U (1) symmetry] of an infinite degeneracy. Surprisingly, this spectral collapse does not signal the presence of an unphysical regime of the parameter space at larger gain rate A. Instead, it indicates a complete change in all the characteristics of the system. This criticality, even in the absence of SSB, can be witnessed by the presence of a dynamical hysteresis. Remarkably, a semiclassical analysis, which neglects field fluctuations, completely misses this phenomenon. Contrary to "usual" multistability, which is suppressed by the presence of quantum fluctuations, this multistability is enabled by fluctuations. For example, in the seminal treatment of the Kerr resonator of Ref. [85], the semiclassical bistable solutions become a single density matrix, once quantum fluctuations are taken into account. This can be explained by the presence of quantum tunneling between the different semiclassical solutions [67,86], and a bistability emerges in the thermodynamic limit around the critical point of a first-order phase transition of the Kerr resonator [18,42,43,79,[87][88][89][90][91].
Our findings also prompt further investigations to analyze the Liouvillian spectral collapse within the framework of dynamical phase transitions and even quantum chaos [92][93][94][95][96]. Indeed, the presence of highly degenerate diabolical points describing nondecaying processes at different frequencies open perspectives for the study of the Loschmidt echoes and other indicators of dynamical criticalities [97,98].
The Landau theory and semiclassical analysis have also been used to study first-order lasing phenomena. For example, early studies of lasers, with a low-saturationintensity saturable absorber [99,100], sparked a number of theoretical [101][102][103][104] and experimental [105,106] works in characterizing first-order lasing transitions. Similarly, first-order transitions in dye lasers were studied [107][108][109], and the analogy between the first-order phase transition and laser threshold for multimode lasers led to the prediction of multicritical points [110][111][112][113]. We plan to use the Liouvillian theory to investigate these classes of models.
Here, we consider the minimal model to obtain the Scully-Lamb laser master equation, inspired by the derivation of Yamamoto and Imamoǧlu in Ref. [59].
We consider a beam of four-level atoms whose states, in order of ascending energy, are {|0 , |g , |e , |f }. The operatorsσ p,q ≡ |p q| indicate the atomic transitions from |q to |p . Spontaneous decays enable the system to decrease its energy via the transitions |q → |p , if |q has a higher energy than |p , described by a set of Lindblad dissipators γ p,q D[σ p,q ] [see Table I and Fig. 8(a)].
A cavity of frequency ω c is at resonance with the |g ↔ |e transition, while all the other transitions are far-off resonance. Under this condition, light and matter can exchange a single excitation, and the coherent part of the system is described by the Jaynes-Cummings Hamiltonian [114] H JC = ω câ †â + ω q 2σ z + g (âσ e,g + h.c.) ,  8. Pictorial representation of the Scully-Lamb laser model obtained via a driven four-level atoms system. (a) A beam of four-level atoms interacts with an optical cavity. An external pump drives the atoms from the lowest energy level |0 to the highest one |f . (b) In the limit in which the decay from |f is rapid, the emerging physics is that of an incoherently driven three-level atom interacting with a dissipative cavity, as detailed in the non-Lindbladian master equation in Eq. (A2). The rates are explained in Table I, and g is the Jaynes-Cummings coupling strength between the atom and cavity. The broadening of |e schematically visualizes the effect of dephasing with a rate γϕ.
whereσ z =σ e,e −σ g,g . The cavity photons are lost via ΓD[â]; the term γ ϕ D[σ z ] describes dephasing due to, e.g., collisions between different species of atoms in the laser, which destroy the coherence between |e and |g induced by the cavity field.
An external coherent field continuously pumps (excites) the atoms from their ground level |0 to the uppermost level |f and vice versa. Let us assume that the drive is strong and that the decay rate γ e,f , describing the transition |f → |e , is faster than all the other rates of the system. Then, one can trace out the |f degree of freedom, obtaining a non-Hermitian effective pumping γ ↑ D[|e g|] [c.f. Fig. 8(a)].
In both cases, by tracing out the atomic degrees of freedom, a non-Lindblad-form for the reduced density matrix ρ(t) ≡ρ(t) cav of the cavity field can be obtained: where A, B 1 , and B 2 are the gain, dephasing, and gain saturation rates. In the fast dephasing case, they are: where γ tot = γ ϕ + γ 0e + γ ge + γ 0g . In the lifetime broadened laser, instead, one has The superoperator K is not a Lindblad dissipator, and it makes Eq. (A3) non-Lindbladian, with where {Â,B} =ÂB +BÂ. In the following, and for the sake of simplicity, we focus on the lifetime broadened laser (the case 2) by choosing B = B 1 = B 2 . Similar results can be obtained for case (1) upon appropriate rescaling of the parameters.

Lindblad form in the weak-gain saturation limit
The master equation in Eq. (A3) can be obtained in several ways [62,64]. In particular, Eq. (A3) is an approximation of the original Scully-Lamb equation in Ref. [3] valid up to third-order corrections in the field amplitude.
Although trace preserving, Eq. (A3) does not conserve the complete positivity of the density matrix (i.e., it can lead to a negative probability distribution of the eigenstates ofρ(t) [61]). Since all Lindblad master equations are completely positive trace-preserving (CPTP) maps, Eq. (A3) cannot be brought to a Lindblad form. However, by ignoring terms of second order inââ † B/(2A) [which corresponds to the WGS limit in Eq. (6)] Eq. (A3) becomes [5,60] whereL j are exactly the jump operators obtained in Eq. (5) (see also Table I). As demonstrated in Ref. [61] for a micromaser model, such a Lindblad form can also by obtained by an expansion in the light-matter coupling strength that explicitly preserves the Lindblad form of the master equation. Equations (A3) and (A7) must be taken with a grain of salt. Since they are valid in the WGS limit, they cannot be used to describe processes involving a high number of excitations. Consider, e.g., an initial Fock stateρ (t = 0) = |m m|. If m is small enough and B A, the leading term inL 1 is √ Aâ † , and the dynamics rapidly converges towards its steady state, where A and B compete. If, however, m 2 A/B, the dynamics is dominated by the B term which is unbounded, and the system cannot converge anymore to the "true" steady state, and the number of photons diverges towards infinity. In other words, there exist values of m for which the states become nonphysical. From a Liouvillian point of view, the SLLM is unbounded and, increasing the cutoff, spurious eigenvalues from the high-excited eigenmatrices of the Liouvillian (which do not respect the WGS conditions) produce quasi-zeros of the spectrum.
In both simulating the dynamics and studying the Liouvillian spectrum, this problem plays a very minor role, since spurious eigenvalues can be easily eliminated. Indeed, when discussing the dynamics of the system, we verified that all the solutions are dynamically stable (i.e., the presence of a diverging dynamics does not affect the steady state). Moreover, in the same way in which the low-lying part of the dynamics is unaffected by the presence of the diverging dynamics for a high-photon number, the Liouvillian eigenstates converge to a well-defined value, even when spurious eigenvalues are present. Thus the unphysical eigenvalues and eigenmatrices can be easily discarded a posteriori.
A more detailed analysis of the problem and a deeper discussion of the validity of this approximation can also be found in Ref. [30]. We recall that the eigenmatrices have, by themselves, no physical meaning but they need to be combined witĥ ρ ss in order to represent a physical process [18]. Whilê ρ ss + j c jρ (0) j is well-defined state, for a generic k one needs to consider whereD α = exp (αâ † − α * â ) is the displacement operator. Notice that, for a generic density matrix, the Wigner representation is the standard Wigner function describing a quasiprobability. For the eigenmatrices under consideration, this is not the case since, except the steadystate, they are traceless. Moreover, the negativities of these Wigner representations of the eigenmatrices are not indicators of the "quantumness" of the state.
In Fig. 9, we plot different eigenmatrices of the Liouvillian, corresponding to different symmetry sectors. Each sector is characterized by a different structure. A matrix corresponding to u 0 must be invariant under any rotation and, therefore, it has a circular symmetry in its Wigner representation [c.f. Fig. 9(a)]. As it can be seen in Fig. 9(b), a state belonging to u 1 must be invariant upon a π rotation (up to the minus sign). Similarly, for Figs. 9(c) and 9(d), the rotation becomes π/2 and π/3, respectively.
We conclude that u 0 describes all those phase-independent exponentially decaying processes. The sector u 1 captures all those processes that have a single lobe in the phase space. This means that they capture the dynamics towards the steady state of those processes rotating at a frequency ω. Similarly, u k is characterized by k identical lobes in the phase space, and it rotates at the frequency kω. To understand why the semiclassical approximation (which often correctly captures criticality in other singlecavity models [43,47,48]) fails for η = 0, we consider here the fluctuations induced by dephasing. The semiclassical approximation assumes the trivial relation n = |α| 2 . However, for symmetry reason, in the full model we have that â ss = exp(iφ) â ss = 0. As such, the value â can be different from zero only in the dynamics towards the steady state or in the thermodynamic limit, where the U (1) SSB takes place. Equivalently, the semiclassical approximation is mixing dynamical quantities with steady-state ones, and it is therefore valid only if η = 0 and coherence can be maintained for an "infinite" time.
Using the previously introduced counting-trajectory formalism, we can characterize the effect of η on a coherent state. If no quantum jump happens, the systemĤ eff contains the term −η/2(â † ) 2â2 , which acts as an imaginary Kerr nonlinearity, inducing dephasing and changing the shape of a coherent state. Similarly, a quantum jumpL 2 destroys the coherence of a state (i.e., the coherent state is not an eigenstate ofâ †â ). These are not the fluctuations of the photon number discussed above.
We plot the results of such a quantum trajectory in Fig. 10 for an initial coherent state. In Fig 10(a), we show the expectation value ofx = (â +â † )/2 both for η = 0 (corresponding to the standard SLLM, red curve) and η = 0 (i.e., the generalized SLLM, blue curve). In both cases, x oscillates at a frequency ω. While for η = 0 the oscillations are long-lived, for η = 0 they rapidly approach zero. This is corroborated by the Wigner functions at time t = 20 Γ in Figs. 10(d) and 10(e). While in the standard Scully-Lamb the system retains a welldefined phase, in the presence of dephasing η = 0 the state rapidly loses any coherent-like feature. Finally, in Fig. 10(c) we plot p( â †â ), i.e., the probability distribution that the photon number attain a certain value along a single very long quantum trajectory, both with or without additional dephasing. Also, this distribution is independent of the value of η.
The results of Fig. 10 seem to indicate that, for a very long trajectory, the system would completely lose its coherence, thus resulting in a state which does not oscillate anymore. One would be tempted to interpret the Pfunction solution as a superposition of states which have completely lost their coherence. This picture, however, is due to the particular assumptions leading to Eq. (32). Indeed, as Eqs. (13) and (14) imply, the dephasing η determines a random-walk phase diffusion of the field. To illustrate this, we can consider a different kind of unraveling (e.g., a diffusive homodyne-like quantum trajectory [82]), the result of the statistical unraveling then would be completely different. In this case, a reference field of intensity β is mixed with the output field of each jump operator. Consequently, the jump operators be-comeL j (β) =L j + β j , while the Hamiltonian readŝ In the ideal "homodyne" trajectory limit β j → ∞, the detector continuously reads a signal. The effect of the measurement on the system, however, is minimal, since each quantum jump is largely due to the presence of the local oscillator β j . We plot these results in Fig. 11. This time, the effect of decoherence can be interpreted as random displacement (walk) which a coherent state undergo [ Fig. 11(a)]. As such, the dynamics of a single quantum trajectory always retain a nonzero phase [c.f. Figs. 11(b) and 11(c)]. In this case, the P -function would result from the superposition of several coherent states with different phases and amplitudes.
We conclude that the photon number and its fluctuation depend only on the gain A, saturation B, and dissipation Γ. We confirm that η plays no role in determining them, but induces decoherence, which can be interpreted in different ways: (i) in a counting trajectory approach, the decoherence changes the shape of an other-wise coherent-like state; (ii) in a homodyne trajectory, η randomly displaces the phase of the wave function, preventing long-term coherence. Neither dephased states nor random fluctuations can be described as a single co-herent state. Thus, the semiclassical approach based on Eq. (7) assigns to η the same role as Γ, wrongly predicting the presence of SSB in the generalized SLLM.