Dissipation-induced bistability in the two-photon Dicke model

The Dicke model is a paradigmatic quantum-optical model describing the interaction of a collection of two-level systems with a single bosonic mode. Effective implementations of this model made it possible to observe the emergence of superradiance, i.e., cooperative phenomena arising from the collective nature of light-matter interactions. Via reservoir engineering and analogue quantum simulation techniques, current experimental platforms allow us not only to implement the Dicke model but also to design more exotic interactions, such as the two-photon Dicke model. In the Hamiltonian case, this model presents an interesting phase diagram characterized by two quantum criticalities: a superradiant phase transition and a spectral collapse, that is, the coalescence of discrete energy levels into a continuous band. Here, we investigate the effects of both qubit and photon dissipation on the phase transition and on the instability induced by the spectral collapse. Using a mean-field decoupling approximation, we analytically obtain the steady-state expectation values of the observables signaling a symmetry breaking, identifying a first-order phase transition from the normal to the superradiant phase. Our stability analysis unveils a very rich phase diagram, which features stable, bistable, and unstable phases depending on the dissipation rate.

One-and two-photon Dicke models. The standard Dicke model was originally used to describe the behaviour of a collection of atoms with the electromagnetic field inside a high-quality-factor cavity. It can be derived by making several assumptions about the system. For instance, the atomic size must be small compared to the field wavelength, making the atoms insensitive to the field modulation 1,65 . The atoms must couple to a single mode of the field. Finally, the atomic energy level structure must be highly anharmonic, so that only one transition is resonant with the field, allowing us to approximate the atoms by two-level systems (or qubits). This so-called two-level approximation, when handled improperly, can lead to a gauge ambiguity in the USC regime. Solutions to this problem have been proposed only recently 8,9 . When these assumptions hold, the system is described by the one-photon Dicke Hamiltonian (here = 1), where â ( â † ) is the annihilation (creation) operator of the bosonic mode, N is the number of qubits, and σ j x,y,z are the Pauli matrices describing the j-th qubit. This Hamiltonian exhibits a Z 2 symmetry, corresponding to the simultaneous exchange For low values of the coupling, the ground state of this Hamiltonian is given by the product state of the field vacuum and the atoms individual ground states. When the coupling constant enters the USC regime, however, the system experiences a second-order phase transition in the thermodynamic limit which breaks the Z 2 symmetry. The system enters the so-called superradiant phase, in which the bosonic field is described by a coherent state, while the qubits are collectively rotated 3 . The possibilities offered by quantum simulations have brought this model far beyond the atom-cavity setting. For instance, the cavity electromagnetic field may be replaced by microwave resonators in superconducting circuits, or by vibrational motion in atomic platforms. These effective implementations made it possible to circumvent the problems raised by gauge ambiguities and to observe the superradiant phase transition [10][11][12]14 . These platforms also pave the way to the experimental exploration of novel forms of light-matter interactions. In particular, quantum simulation schemes make it possible to implement two-photon interaction models both in the SC and in the USC regime 42 www.nature.com/scientificreports/ can be used to couple the internal state of the ions to their motional degrees of freedom. Let us assume that the properties of the trap allow us to single out a single vibrational mode with frequency ν . If the detuning between the laser and the internal transition is close to −2ν (red-detuned laser), then the laser can excite a process in which two phonons are destroyed and one qubit excitation is created. If the detuning is close to 2ν (blue-detuned laser), then the energy brought by the laser can be used to simultaneously create one qubit excitation and two phonons. Therefore, by using both a red-detuned and a blue-detuned lasers, one can engineer a qubit-boson coupling similar to the Dicke model, but where the standard one-boson interaction term is replaced by a twoboson term, which is generically called in the literature two-photon or two-phonon coupling term. For simplicity, we will only use the term "two-photon" and the cavity-QED terminology in the following. Furthermore, the modulation of photonic states by the laser pump permits to effectively renormalize both the bosonic frequency and the coupling constant, thus allowing to bring the two-photon coupling to the USC regime.
Similarly, it has recently been shown 44,45 that two-photon interactions can also be implemented in superconducting circuits, engineering an intrinsic nondipolar coupling between a superconducting artificial atom and superconducting quantum interference device (SQUID). In this case, the standard linear coupling is suppressed, while the two-photon coupling terms emerge as the natural light-matter interaction in an undriven system and not as the result of a quantum simulation scheme.
These various possibilities of implementing the two-photon coupling term motivates the study of the twophoton Dicke model, whose Hamiltonian reads (setting = 1), This Hamiltonian exhibits a four-folded symmetry, stemming from the simultaneous exchange of In the USC limit, this model exhibits an instability known as spectral collapse 42 , where the discrete spectrum collapses into a continuous band for a critical value of g. Some intuition about this effect can be gained through the following reasoning. When the coupling constant g in Eq. (3) becomes large, the interaction term dominates the physics. Since this term commutes with the σ i x , we can study the qubits domains σ i x = − 1 2 and σ i x = + 1 2 independently. Let us consider σ i x = − 1 2 for all i. Then we have an effective boson dynamics described by this Hamiltonian which is a quadratic potential for the field quadratures x =â † +â and p = i(â † −â) . When g is large enough, this potential becomes almost flat, shrinking the gap between the different energy levels. Ultimately, these levels coalesce into a continuous band, causing the so-called spectral collapse 42 . When g is increased even further, the potential becomes an upside-down harmonic well, and is unbounded from below for x → ∞ . Therefore, the dynamics of the system will become unstable, signaling the breaking down of the model. By contrast, in the one-photon Dicke model (1), the interaction term adds only a linear correction, meaning that the Hamiltonian can never be unbounded from below.
Very recently it has been shown 60 that the two-photon Dicke model can also display a second-order quantum phase transition very similar to the superradiant transition of the one-photon Dicke model. Instead of a coherent state, however, the bosonic field here will be described by a squeezed state for high values of the coupling. The ground-state phase diagram has been analyzed with different numerical and analytical techniques also for other two-photon coupling models 63,64 . However, two-photon light-matter interaction models have so far never been considered from an open quantum system perspective.

Effect of dissipation.
The physics of the one-photon Dicke model changes drastically once dissipation is taken into account. It was shown in Refs. 25,69 that in the presence of qubit decay and dephasing, the transition of the Dicke model could be modified, suppressed, or restored. Similarly, the presence of dissipative processes in the two-photon Dicke model raises intriguing questions.
Assuming a Markovian environment and performing the Born approximation, the dissipation may be described by a Lindblad master equation 28 .
In our analysis, we will assume that the Hamiltonian part of the evolution remains that of Eq. (3), and we will include three dissipation channels, that is, individual qubit decay and dephasing, and photon loss. We obtain the following Lindblad equation (we recall = 1), where ρ(t) is the density matrix of the system at time t, σ   (4). In this regard, also the Lindblad master equation presents a four-folded symmetry similar to that of the Hamiltonian case. Let us note that this equation is a purely phenomenological one and, while it is the most appropriate in a quantum simulation framework, it would fail to describe the true evolution for genuine implementations of the model in the USC regime . Indeed, in the presence of bare local dephasing and qubit decay processes, the system would not tend towards the dressed ground state. In fact, these processes would effectively pump energy into the system, forcing it away from the true polaritonic ground state and toward a different steady-state; while considering dressedoperator incoherent processes would lead to the polaritonic ground state 29 . Moreover, a microscopic theory needs to be developed for arbitrary strengths of the dissipative couplings, which has shown that USC effects can be robust in loss-dominated systems 71 . However, our analysis is meant to describe effective implementations of the model, where the considered decoherence and dissipation processes can themselves be implemented via bath-engineering techniques 67 . For instance, this Lindbladian dynamics could be observed in a strongly driven atomic cloud, where Raman processes effectively engineer the wanted processes as demonstrated in several experiments [10][11][12] . Note also that, since quantum simulation allows to renormalize both the effective frequencies ω c and ω 0 and the coupling constant g 42 , the dissipation constant may be large compared to ω c , ω 0 and g, while remaining small compared to the actual frequencies of the system. We expect that this model will have very different behavior from its Hamiltonian counterpart. On the one hand, the critical behaviour and the properties of the superradiant phase can drastically change, as for the onephoton model 5,25 . On the other hand, the spectral collapse may be modified or avoided due to the presence of dissipation. Indeed, the photon loss term in Eq. (6) acts like a stabilizing quadratic term which can balance the effect of the Hamiltonian unstable potential.

Results
Symmetry breaking. When spin dissipation is absent 3,21,23,27,69 or acts collectively on all qubits at once, one can significantly simplify the problem by treating the qubits as a single, collective spin, which allows to obtain the quantum state of the qubits. This is not possible here, since the dissipation acts on each qubit individually 25,26,29 . However, we can still obtain meaningful results by focusing on some specific, relevant observables. We have focused our analysis on the following quantities: Hamiltonian treatment 60 predicts a symmetry breaking during which the bosonic field becomes squeezed, which is captured by the second-order moment of the bosonic field. Therefore, we expect that the observables X , Ŷ and â †â will be valid order parameters also in the presence of dissipation. We have studied the evolution of these quantities using a mean-field decoupling approximation (see the Methods section). We found that the dynamics of these quantities has three possible solutions. The first one corresponds to �â †â � = �X� = �Ŷ � = �Ĵ x � = �Ĵ y � = 0 and �Ĵ z � = −1 . The other two phases have �â †â � , �X� = ±X s , �Ŷ � = ±Y s � = 0 (complete expressions in the Methods section, Eq. (10)). In accordance with previous results in the one-and two-photon Dicke model 3,5,25,60 , we can identify the first solution as the "normal phase", as it corresponds to the product state of the individual ground states of the field and the atoms. The other two solutions correspond to the "superradiant phase"; that is, they contain a macroscopic number of atomic and photonic excitations. Since the "superradiant" solutions have �X� � = 0 , the four-folded symmetry of the model is at least partially broken. The stability of each phase has been analyzed numerically (see Methods). Let us now illustrate the properties of the superradiant phase and the complex driven-dissipative phase diagram of the model considered.
Nature of the phase transition. In Fig. 1 we show the value of the steady-state photon number in the superradiant phase as a function of the coupling strength g, ω 0 = ω c , and κ = ω c . For now, we have set Ŵ ↓ = Ŵ φ = Ŵ , and Ŵ = 3ω c . For small values of g, the superradiant phase yields nonphysical complex values for �â †â � , showing that the system can only reach the normal phase �â †â � = 0 until the critical value of the coupling strength is achieved. When g is increased, the superradiant phase becomes physical, and the stability analysis reveals it is stable as well.
Therefore, a nonzero number of bosonic excitations can appear in the system. Analyzing the average photon number in the system steady-state, we can already identify two qualitative differences compared to the ground state in the Hamiltonian case. First, in the driven-dissipative case, the number of photons in the superradiant phase does not go to zero when one approaches the limit of stability from above. Second, the point at which the normal phase becomes unstable and the superradiant phase becomes stable do not coincide. Therefore, the driven-dissipative two-photon Dicke model exhibits bistability at the mean-field level.
The emergence of bistability in mean-field models is well-known in open quantum systems. A typical example is that of the Kerr resonator, where the semiclassical solution obtained via the Gross-Pitaevskii mean-field has three different solutions: two which are stable and one unstable. As soon as one considers the quantum steadystate, however, only one solution is found 72 . This apparent contradiction can be solved by considering the full Liouvillian spectrum, where the onset of bistability is in close relation to the emergence of criticality 19,73,74 . Indeed, several models presenting bistable behaviour at the mean-field level proved to display a genuine first-order phase transition in the thermodynamic limit of a full quantum model 56,[75][76][77][78][79][80][81][82] .
These results show that the Hamiltonian and dissipative versions of this model are strikingly different. In the equilibrium case, a second-order phase transition is predicted to occur, and only in the far-detuned regime 60 www.nature.com/scientificreports/ ω 0 ≪ ω c . In the nonequilibrium case, a first-order phase transition takes place in the resonant regime ω 0 = 2ω c , a condition that strongly simplifies possible experimental implementations.

Phase diagram.
Having established the existence of a phase transition, we can produce the phase diagram of the model by studying the stability of both phases for a broad range of parameters. The analysis of these diagrams revealed the existence of two regimes of dissipation. In Fig. 2, we display the phase diagram in the g-ω 0 plane,  N). In all cases, the superradiant phase becomes stable before the normal phase becomes unstable, indicating bistability. www.nature.com/scientificreports/ for two values of Ŵ : Ŵ = 1.5ω c and Ŵ = 3ω c , and for various number of qubits N. For Ŵ = 1.5ω c and the smaller value N = 10 qubits, we observe that the mean-field equations predict the existence of a zone where the superradiant phase is stable. However, the size of this zone shrinks when N increases. Since the mean-field description becomes correct only for N → ∞ , no phase transition can happen in the mean-field limit for this value of dissipation. The system will either reach the normal steady-state or be unstable. For Ŵ = 3ω c , however, we observe that bistability becomes possible, similarly to what happens in the driven-dissipative one-photon Dicke model 83 . In the thermodynamic limit, the region of stability becomes independent of the number of qubits, meaning that a phase transition can take place in the N → ∞ limit.
Values of Ŵ/ω c lower that 1.5 or higher than 3 yield qualitatively similar results, which allows us to conclude that there are two regimes of dissipation: a large dissipation regime in which a phase transition is possible, and a low dissipation regime in which only the normal phase is stable in the thermodynamic limit.
Interestingly, the transition between these two regimes of parameters when Ŵ is increased is quite sharp, especially in the thermodynamic limit. To visualize this, we study the stability of the superradiant phase versus both g and Ŵ , for 100 qubits, and for various values of ω 0 , the other parameters being the same (this amounts to taking horizontal slices in Fig. 2 and study their evolution when Ŵ changes).
The results are displayed in Fig. 3: for Ŵ/ω c ≈ 1.6 , the instability disappears and the superradiant phase becomes stable for most values of ω 0 and g. Hence, the phase diagram as a whole changes drastically when Ŵ/ω c goes across this threshold.
Hence, we have established that the presence of dissipation is instrumental in stabilizing the superradiant phase. If we compare this with the results obtained in the one-photon version of the driven-dissipative model 25 , an instructive analogy can be made. Adding enough qubit dissipation appears to preserve the superradiant phase transition, which normally would be destroyed in the presence of noise. In the one-photon case, however, decay and dephasing can play antagonistic roles: adding an infinitesimal amount of qubit dephasing without decay destroys the transition, while adding both dephasing and decay stabilizes it. To see if such effect is also present in the two-photon model, we study the stability of the superradiant phase with respect to both Ŵ φ and Ŵ ↓ . The results are displayed in Fig. 4. We see no evidence of suppression and restoration of the phase transition. Rather, these plots indicate that both dephasing and decay contribute positively to the stabilization of the superradiant phase.

Discussion
In this paper, we present the first analysis of the steady-state phase diagram of a two-photon interaction model in the driven-dissipative case. In particular, we have explored numerically the mean-field behavior of the N-body two-photon Dicke model. We have identified a rich behavior, including a superradiant phase transition of first order, a bistable phase, and an instability that is removed by dissipation. Although one may be tempted to interpret this instability as the dissipative counterpart of the spectral collapse, there are important differences between the two phenomena.
In the Hamiltonian case, the spectral collapse is expected to occur when g increases. In our case, for Ŵ = 3ω c , we have increased g towards higher values, up to 10 3 (not shown). We have found only a stable superradiant phase. To summarize, the spectral collapse is entirely controlled by the coupling, while the instability is controlled by a non-trivial interplay between dissipation and coupling. Another, perhaps even more important, difference is the scaling of the parameters with N. In the Hamiltonian case, the collapse occurs for g ∼ ω c √ N 42,60 (note that the definition of the coupling constant is not the same in our work and these references). The stable dynamics can only take place in the interval 0 ≤ g ≤ ω c √ N , which becomes vanishingly small in the thermodynamic limit. As a consequence, a transition can take place only if the qubit frequency ω 0 is also allowed to scale like 1/ √ N 60 . www.nature.com/scientificreports/ By contrast here, the instability, when present, only emerges for g ∼ ω c , which allows a phase transition to take place for all values of ω 0 . Indeed, the phase diagram we have obtained here are quite different from the one of the Hamiltonian model published in 60 . The fact that the spectral collapse and the instability scale differently with the number of atoms suggests that they could be qualitatively different processes. Hence, the two main properties of the two-photon Dicke Hamiltonian, second-order phase transition and spectral collapse, are qualitatively modified in the presence of dissipation. This illustrates how the critical behavior of a given phase transition can change radically when one goes from the equilibrium case to the non-equilibrium one.
Furthermore, we have pointed out conceptual differences between the behavior of the one-and two-photon Dicke models. For the latter, both local dephasing and decay appear to help stabilize the transition, in contrast to what was found in the dissipative one-photon Dicke model. In the thermodynamic limit of the two-photon Dicke model, the entire phase diagram changes very abruptly when dissipation is moved across a very narrow range of parameters.
We note that while a mean-field approach can predict several solutions for the steady-state equation in the bistable region, and a unique solution in the mono-stable regions, in a finite-size quantum system there can only be a unique steady-state in each region. Indeed, the introduction of both quantum and classical fluctuations prevents the fields from remaining stationary around their mean-field values. Several works 73,84 have illustrated these behaviours in the similar case of the two-photon Kerr resonator. In turn, the presence of multiple solutions at the mean-field level is translated into an observable bistable behavior of the critical parameters of the full quantum model in a quantum trajectory approach 85 .
Indeed, for future perspectives, a full quantum treatment of the driven-dissipative two-photon Dicke model would be an interesting and yet challenging task. Extracting information on the thermodynamic limit from direct numerical simulations of the full dynamics is far from straightforward, due to the exponentially increasing Liouvillian space. Exploiting the permutational invariance of the Liouvillian 29 can exponentially reduce the computational overhead with regard to the qubit degrees of freedom, but the photonic subspace, approximated by a cut-off photon excitation number n ph , needs to be larger than in the case of the single-photon Dicke model to avoid spurious results induced by the finite approximation of the otherwise unbounded Hilbert space. With this regard, two possible solutions could be considered, also simultaneously: first, the inclusion a two-photon dissipation process, which would effectively reduce the highest excitation number explored, for appropriately large values of the two-photon decay rate; second, the use of quantum trajectories, which reduces the computational overhead from being that of the Liouvillian space to just that of an effective Hilbert space, at the cost of averaging over many runs. We point out though that the intermittent dynamics characterizing a bi-stable phase can be grasped even by single quantum trajectory simulations. Alternatively to these approaches, a qubit-only description of the system could be obtained at the cost of abandoning the Lindbladian formalism for a full Redfield theory 86 .
From an experimental perspective, we have shown that in the driven-dissipative case the superradiant phase transition of the two-photon Dicke model can be observed also for resonant interactions, and does not require large detuning as in the static case.
Hence, this quantum phase transition could be observed in a more accessible regime of parameters than previously thought. The driven-dissipative two-photon Dicke model could be implemented with trapped ions 42 or in a cold atoms setup, similar to what has been done already for the driven Rabi model 87,88 . Superconducting circuits 48 provide another platform to simulate this dynamics, in which the two-photon interaction can be engineered between a flux qubit and a SQUID resonator 44,45 . www.nature.com/scientificreports/

Methods
From the Lindblad master equation (6), we can obtain the field equation governing the dynamics of any operator ∂ t �Â� . For the operators X , Ŷ , â †â , Ĵ x , Ĵ y , Ĵ z that we have studied, this gives: The solution of Eq. (8) is, in general, a formidable task. If one is interested in the properties of the steady-state, however, the time derivatives can be set to zero. This approximation is not sufficient to solve Eq. (8), since some operators are a function of higher-order correlation functions, thus resulting in an infinite hierarchy of coupled equations. In the normal phase (i.e., when no symmetry is broken), one can reduce the complexity of the problem by considering κ ≫ Ŵ ↓ ≃ Ŵ φ , i.e., that the bosonic field reaches a steady-state long before the qubits do. Indeed, in this region the Liouvillian gap must be opened 19 , and in the absence of critical-slowing down the typical timescale is dictated by the dissipation rates. Using adiabatic elimination, one can easily find the behavior of the system close to the normal phase characterised by �X� = �Ŷ � = �â †â � = �Ĵ x � = �Ĵ y � = 0 , while �Ĵ z � = −1 . The results of adiabatic elimination, however, fail to capture the superradiant phase: in this regime, the diverging timescale coming from the closure of the Liouvillian gap makes the photonic timescale comparable to the qubit one.
In order to truncate the hierarchy of equations stemming from a Liouvillian problem, in many-body quantum physics one often resorts to a Gutzwiller mean-field approximation. In this case, one assumes that the system density matrix can be factorised as a tensor product between the qubit and photonic part, decoupling all the high-order correlation in Eq. (8). For instance, we will assume that �Ĵ xâ †â � = �Ĵ x ��â †â � . Note that in general, the mean-field approximation can lead to incorrect predictions in the presence of strong quantum correlations. It is expected to be true for high-dimensional models (i.e., when the number of nearest neighbors is elevated) and in the thermodynamic limit N → ∞ . In Dicke-like models, all spins are coupled to the mode, which induces an effective all-to-all coupling. In the thermodynamic limit, the model is effectively infinite-dimensional, and the qubits act as a collective, classical spin. Furthermore, in the superradiant phase, the fluctuations of both the spins and the field are negligible compared to their mean-field values 89 . As a consequence, it was shown that the mean-field approximation gives indeed correct results for Dicke-like models with large N 26,65 . Under the meanfield approximation, Eq. (8) becomes One solution to this equation is �X� = �Ŷ � = �â †â � = �Ĵ x � = �Ĵ y � = 0 , while �Ĵ z � = −1 . That is, the normal phase with no photons is always a solution to Eq. (9). However, there are other two solutions to this equation, where the bosonic field is populated. Namely: where we have introduced β = ω c Ŵ ′ 2ω 0 NŴ ↓ and g t = 2ω c + κ 2 2ω c 2ω 0 + Ŵ ′2 2ω 0 /8. Having identified the three possible solutions, we study their stability by considering a linear perturbation of the steady-state value: For the normal phase, we have − → A = �X�, �Ŷ �, �â †â �, �Ĵ x �, �Ĵ y �, �Ĵ z � www.nature.com/scientificreports/