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.

Introduction

The Dicke model[1] describes the interaction of a collection of two-level systems with a single bosonic mode. In the thermodynamic limit, this model exhibits a superradiant phase transition at zero temperature[1-5]. Namely, the ground-state number of photons changes non analytically from zero to finite values as the light-matter coupling strength is increased across a critical value. The relevance of the Dicke model to capture the physics of light-matter coupling near the critical point is the object of an ongoing debate. In particular, an obstacle to the observation of a superradiant phase transition arises due to the presence of the so-called diamagnetic term. In this regard, resolutions of gauge ambiguities have been recently proposed, such as employing modified unitary transformations or going beyond the two-level system description[6-9].

It is possible, however, to circumvent this controversy entirely by using driven systems and bath engineering to simulate effective Hamiltonians. Thanks to this approach, it was possible to observe the superradiant phase transition in several platforms, in particular atomic systems in cavity[10-12], and trapped ions[13,14]. Other proposals have been put forward, using NV centers array[15], or superconducting circuits[16]. In general, the implementation of analogue quantum simulations[17,18] provides an ideal playground to test driven-dissipative physics in a controlled setting. Their experimental feasibility has motivated increasing research efforts devoted to the study of driven-dissipative quantum optical models, as it is known that noise and dissipation can drastically change the properties of the steady-state phase diagrams and the emergence of phase transitions[19,20]. Analytical studies[21-27] are extremely challenging, as the intrinsic non-equilibrium nature of driven-dissipative systems does not allow a determination of the stationary state of the system via a free energy analysis[28]; for Dicke-like models, it is especially true when local spin dissipation is considered[25,26,29].

Effective implementations offer others significant advantages. First, in a cavity-QED setting, driving the system allows turning virtual excitations inside the cavity into real photons which can exit the system[30], thus giving immediate access to the intra-cavity dynamics. Second, by manipulating the exchanges between a system and its environment, reservoir-engineering techniques allow us to realise previously inaccessible quantum phases of matter[31-37]. For instance, it is possible to stabilize phases without an equilibrium counterpart[38, 39], and reservoir engineering methods for complex many-body phases have been thoroughly explored in different contexts[32,40]. Several experiments and theoretical proposals have applied these ideas to study generalized Dicke models[11,15,41-45] and the ultrastrong coupling regime (USC)[46,47], i.e., the regime of parameters where the coupling constant becomes a sizable fraction of both the qubit and bosonic frequencies.

Among the variety of phenomena that are made accessible by analogue quantum simulations[17,18], a particularly interesting one is the possibility to engineer a coupling involving the simultaneous exchange of several photons. In superconducting circuits[48], this is possible with nonlinear virtual processes in the USC regime[49-53], and experimentally a photon-pair driving mechanism has been realized[37], leading to the generation of so-called photonic Schrödinger cats[54-56]. The possibility to control and protect such states is promising for the implementation of quantum computation protocols[57-59].

In the two-photon Dicke model[42,44,60], it is the light-matter interaction that creates or annihilates one pair of bosonic excitations per qubit flip. This exotic exchange leads to several unusual properties in the USC regime[46,47]. In particular, it was shown very early on[61,62] that there exists a critical value of the coupling strength for which the discrete spectrum collapses into a continuous band[42]. For higher values of the coupling, the two-photon Dicke Hamiltonian is no longer bounded from below, indicating the breaking down of the model itself. It was recently shown that, in spite of the spectral collapse, a superradiant-like phase transition can take place also in the two-photon Dicke model[60]. This transition has been characterized also for other two-photon interaction models[63,64] but, so far, only the Hamiltonian case has been considered.

In this paper, we investigate the fate of the spectral collapse in a driven-dissipative case. Specifically, we consider the N-body two-photon Dicke model connected to an engineered Markovian bath. The Lindblad formalism is used to introduce different incoherent processes, such as photon loss , individual qubit decay , and local qubit dephasing . To analyze the system dynamical properties, we resort to a mean-field decoupling of the equation of motion, allowing us to determine (semi-)analytically the steady-state of the system. We show the emergence of a first-order phase transition from a normal to a superradiant phase, for a critical value of the light-matter coupling g. A numerical study of the stability of the different phases unveils a very rich phase diagram: depending on the strength of the atomic dissipation, the system may be stable, bistable, or unstable. In particular, the spectral collapse can disappear altogether. Interestingly, we find that both atomic decay and atomic dephasing can be beneficial to the stabilization of the superradiant phase. We found no evidence of supression and restoration of the transition, contrary to the phenomenology of the standard Dicke model[25].

This paper is organized as follows: first, we introduce the one- and two-photon Dicke model, and briefly describe some of their properties using heuristic arguments. Secondly, we discuss dissipation in the two-photon model. Thirdly, we present the phase diagram of the system obtained via a decoupling mean-field approximation, and show the existence of two different regimes of dissipation. Finally, we summarize our conclusions and present some perspectives for future work.

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 ),where () is the annihilation (creation) operator of the bosonic mode, N is the number of qubits, and are the Pauli matrices describing the j-th qubit. This Hamiltonian exhibits a symmetry, corresponding to the simultaneous exchangeFor 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 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-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,66-68]. For instance, in trapped-ions experiments, laser-induced interactions 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 (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 (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 two-boson 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 two-photon Dicke model, whose Hamiltonian reads (setting ),This Hamiltonian exhibits a four-folded symmetry, stemming from the simultaneous exchange ofIn 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 , we can study the qubits domains and independently. Let us consider for all i. Then we have an effective boson dynamics described by this Hamiltonianwhich is a quadratic potential for the field quadratures and . 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 . 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 ),where is the density matrix of the system at time t, , indicates the commutator, and are the Lindblad dissipation superoperators defined as[28]Even if, in general, the symmetries of the Lindblad equation cannot be directly obtained from those of the Hamiltonian[70], Eq. (6) remains identical under the transformation given in Eq. (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 dressed-operator 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-12]. Note also that, since quantum simulation allows to renormalize both the effective frequencies and and the coupling constant g[42], the dissipation constant may be large compared to , 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 one-photon 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: , , , . The 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 , 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 and . The other two phases have , , (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 , 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, , and . For now, we have set , and . For small values of g, the superradiant phase yields nonphysical complex values for , showing that the system can only reach the normal phase 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.

Figure 1

Photon number in the superradiant phase versus normalised light-matter coupling, for several qubit numbers N. We have set and . For small g, the superradiant phase becomes unphysical. The horizontal black line is a guide to the eye indicating the value , which is the value the field adopts in the normal phase. The vertical dashed line indicates the point where the normal phase becomes unstable (this point is independent of N). In all cases, the superradiant phase becomes stable before the normal phase becomes unstable, indicating bistability.

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 steady-state, 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-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]. In the nonequilibrium case, a first-order phase transition takes place in the resonant regime , 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- plane, for two values of : and , and for various number of qubits N. For and the smaller value 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 , 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 , 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 limit.

Figure 2

Phase diagram of the model. a) qubits, , b) qubits, , c) qubits, , d) qubits, , e) qubits, , and f) qubits, . For all plots, . N: normal phase, S: superradiant phase, B: bistable phase, I: instability. For , the region of stability for the superradiant phase shrinks when the number of qubits increases. For , a bistable behavior is observed, and the phase diagram is almost invariant when the number of qubits is increased beyond a few dozens.

Values of 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 , 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 , the instability disappears and the superradiant phase becomes stable for most values of and g. Hence, the phase diagram as a whole changes drastically when goes across this threshold.

Figure 3

Stability of the superradiant phase versus coupling and dissipation. a), b), c). For all plots, and qubits. S,B: superradiant or bistable phase; N,I: normal phase or instability. Here we do not distinguish between the normal and the unstable phase or between superradiant and bistable phase because we want to focus on the stability of the superradiant phase only. Except for small values of or g, the value of for which the transition occurs is almost the same for all parameters, around 1.6. This means the phase diagram as a whole changes drastically when goes across this value.

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.

Figure 4

Stability of the superradiant phase versus qubit dephasing and decay. a), b), c). For all plots, , , and qubits. S,B: stable or bistable phase, N,I: normal phase or instability. Here, 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 , we have increased g towards higher values, up to (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 [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 , which becomes vanishingly small in the thermodynamic limit. As a consequence, a transition can take place only if the qubit frequency is also allowed to scale like [60]. By contrast here, the instability, when present, only emerges for , which allows a phase transition to take place for all values of . 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 , 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].

Methods

From the Lindblad master equation (6), we can obtain the field equation governing the dynamics of any operator . For the operators , , , , , that we have studied, this gives:where we have defined . 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 , while . 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 . 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 . 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 mean-field approximation, Eq. (8) becomesOne solution to this equation is , while . 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 and .

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 havewhile for the superradiant phaseOnly if all the eigenvalues of () are negative, the normal (superradiant) phase is stable. We have obtained these eigenvalues numerically to study the phase diagram of the system.

Note that in the thermodynamic limit , the solution (10) can be simplified as:The stability matrix can be correspondingly simplified. However, even in this case, the dependence in N does not disappear; furthermore, exact solutions are still out of reach. We have been able to check, however, that these exact results, for large N, showed consistency with the thermodynamic limit expressions (in particular, in both cases, the superradiant phase does disappear for small dissipation).