The nonlinear semiclassical dynamics of the unbalanced, open Dicke model

In recent years there have been significant advances in the study of many-body interactions between atoms and light confined to optical cavities. One model which has received widespread attention of late is the Dicke model, which under certain conditions exhibits a quantum phase transition to a state in which the atoms collectively emit light into the cavity mode, known as superradiance. We consider a generalization of this model that features independently controllable strengths of the co- and counter-rotating terms of the interaction Hamiltonian. We study this system in the semiclassical (mean field) limit, i.e., neglecting the role of quantum fluctuations. Under this approximation, the model is described by a set of nonlinear differential equations, which determine the system's semiclassical evolution. By taking a dynamical systems approach, we perform a comprehensive analysis of these equations to reveal an abundance of novel and complex dynamics. Examples of the novel phenomena that we observe are the emergence of superradiant oscillations arising due to Hopf bifurcations, and the appearance of a pair of chaotic attractors arising from period-doubling cascades, followed by their collision to form a single, larger chaotic attractor via a sequence of infinitely many homoclinic bifurcations. Moreover, we find that a flip of the collective spin can result in the sudden emergence of chaotic dynamics. Overall, we provide a comprehensive roadmap of the possible dynamics that arise in the unbalanced, open Dicke model in the form of a phase diagram in the plane of the two interaction strengths. Hence, we lay out the foundations to make further advances in the study of the fingerprint of semiclassical chaos when considering the master equation of the unbalanced Dicke model, that is, the possibility of studying a manifestation of quantum chaos in a specific, experimentally realizable system.


I. INTRODUCTION
Since the advent of chaos theory in the 1960's, the question of how classical chaotic behavior arises from quantum dynamics has been a topic of significant interest.After all, the linearity of the formulation of quantum mechanics disallows the existence of exponential sensitivity to initial conditions that is characteristic of classical chaotic dynamics [1].In the late 20th century, the study of the onset of chaotic behavior in quantum mechanical models found progress in the field of many-body quantum electrodynamics (QED).More recently, quantum chaos has been extremely topical due to its connection with information scrambling [2,3] and black hole physics [4,5]; it has even been suggested to provide routes to experimentally testing the holographic principle of quantum gravity [6,7].
A quantum mechanical model that has recently brought attention to a plethora of uniquely many-body effects is the Dicke model.This paradigmatic model of quantum optics describes the interaction of an ensemble of N two-level atoms with a light field confined to an optical cavity [8][9][10].The Dicke model has been of particular interest because of a quantum phase transition that the model can exhibit due to the onset of significant cooperative interactions of the atoms with the light field [10][11][12].If the strength of the atom-field coupling is suf- * ksti263@aucklanduni.ac.nz ficiently high, the system can undergo a quantum phase transition from the so-called normal phase into a superradiant phase, which features a non-zero photon number at equilibrium [13][14][15][16].
The Dicke model has been used to study chaos in many-body QED in the semiclassical regime [17][18][19][20] and in relation to the superradiant phase transition in a seminal pair of papers by Emary and Brandes [21,22].However, these treatments depend critically upon the Hamiltonian formulation of quantum mechanics and are, thus, limited to closed quantum systems.A more realistic approach that includes dissipative effects introduces a quantum master equation [10]; such a model is referred to as the open Dicke model.Recently, chaotic behavior was uncovered in a driven version of this model [23].
The superradiant quantum phase transition described by the Dicke model has proved difficult to realize experimentally.However, technological advancements over the past decade have led to a number of successful realizations.Baumann et al. [24,25] and Klinder et al. [16] have produced an experimental realization of the Dicke model with a Bose-Einstein condensate coupled to an optical cavity.With this setup, the Dicke superradiant phase transition is produced when the Bose-Einstein condensate undergoes a symmetry breaking self-organization into a supersolid phase.Other realizations have been achieved by Hamner et al. with a spin-orbit-coupled Bose-Einstein condensate confined in a harmonic trap [26] and by Safavi-Naini et al. with a trapped-ion array [27].
On the optical side, Dimer et al. [28] proposed engi-  [28] in a cavity QED system.An optical cavity with a resonant frequency ω and decay rate κ is coupled to N atoms with frequency splitting ω0.The atoms are also driven by two independent lasers, which produce two independent Raman coupling strengths λ±.Each atom is a four-level system with two stable ground states |↓ and |↑ .The laser and cavity fields are sufficiently detuned from the atomic transitions that the two excited states |+ and |− can be adiabatically eliminated in the theoretical model, producing an effective two-level system with comparable rotating and counter-rotating atom-cavity interaction terms.
neering an effective Dicke model in a cavity QED experiment with atom-photon interactions mediated by Raman transitions between ground atomic spin states of a four-level system.This scheme is illustrated in Fig. 1.The proposal considers a four-level atom with two ground spin states, |↓ and |↑ , and two excited states |− and |+ , highly detuned from the cavity and laser frequencies.Two independent lasers drive the transitions |↑ ↔ |− and |↓ ↔ |+ , respectively, while the cavity mode couples simultaneously to the transitions |↑ ↔ |+ and |↓ ↔ |− .Together, the lasers and cavity mode drive two, independent Raman transitions between the ground spin states and, if the detunings of the fields from the excited states are sufficiently large, these states can be adiabatically eliminated to yield an effective twostate-atom model.Significantly, in this effective model counter-rotating and co-rotating interactions feature on an equal footing; i.e., the excitation of the atom from |↓ to |↑ may occur with either absorption of a photon from, or emission of a photon into, the cavity mode, as a result of the two possible Raman transition paths between |↓ and |↑ .The generalization of the Dicke model to feature two independent coupling strengths for the co-and counter-rotating interactions, which we refer to as the unbalanced, open Dicke model, is the focus of this paper.It was recently realized by Zhiqiang et al., whose experimental and theoretical results have shown that this generalization yields interesting nonequilibrium phase behavior and, in particular, oscillatory phases [29,30].
In the work presented here, we take a dynamical systems approach to analyze the set of nonlinear differential equations that arise from the semiclassical model.In doing so, we demonstrate that, in contrast to equal coupling, the option of unbalanced coupling reveals a wide variety of dynamical behavior.In particular, we analyze in detail the semiclassical description of a superradiant phase transition that occurs in two stages, corresponding to the onset of superradiance and the subsequent disappearance of the normal phase.In between the two stages, the system features a phase coexistence of the normal and superradiant phases, as noted in [30].We are also able to use our approach to explain the results of Zhiqiang et al. [29], showing the appearance of superradiant oscillations as arising from Hopf bifurcations.In addition to this, we describe a new type of transition that occurs between the two types of normal phases; spin-down, where the atomic population is entirely de-excited, or spin-up, where it is entirely excited.We then show that, at the transition point, an infinite set of energy-conserving periodic orbits emerge on the Bloch sphere for non-zero dissipation rates.We also analyze the onset of chaos in the system, demonstrating period-doubling cascades and collisions of chaotic attractors.To summarize all of the behaviors described above, we present an overview of the possible dynamics of the model in the form of a phase diagram in the plane of the two coupling strengths that appear due to the generalization of unbalanced coupling.Our exposition provides an original and comprehensive, nonlinear dynamics perspective on the unbalanced Dicke model, complementing other recent interest in this quite fascinating system [30,31].
The paper is organized as follows.In Sect.II, we outline the unbalanced, open Dicke model and its semiclassical approximation as a set of ordinary differential equations.In Sect.III, we describe the possible dynamical phases and their mathematical description as solutions to said differential equations in the semiclassical regime.In Sect.IV, we present a bifurcation analysis of the steadystates and their multistability, where we demonstrate the two-stage superradiant phase transition.In Sect.V, we discuss multistability between phases in more detail and show the basins of attraction for the normal and superradiant phases on the Bloch sphere for an initially empty cavity.In Sect.VI, we analyze the emergence of superradiant oscillations, as demonstrated experimentally in [29].In Sect.VII, we explore the transition that arises between the two types of normal phases, before moving onto an analysis of the emergence of chaotic dynamics in Sect.VIII.Conclusions and an outlook are presented in Sect.IX.

II. UNBALANCED DICKE MODEL
We consider an ensemble of N identical two-level atoms interacting cooperatively with a single mode of the radiation field.The unbalanced Dicke model features two different coupling strengths, λ ± , for the co-and counterrotating terms in the Hamiltonian (with = 1), where â is the annihilation operator of the cavity field mode, ω is the frequency of the field mode, and ω 0 is the frequency splitting between the atomic levels.Ĵ±,z are the collective angular momentum operators, given by where σ(ν) ±,z are the (spin-1/2) angular momentum operators for the νth atom.The collective operators form an angular momentum algebra such that [ Ĵ+ , Ĵ− ] = 2 Ĵz , [ Ĵz , Ĵ± ] = ± Ĵ± , and they have a linear degenerate subspace spanning the 2 N -dimensional Hilbert space [9].A basis is {|J, M }, where with The states |J, M are referred to as the Dicke states.Here we consider the case of maximal angular momentum, J = N/2.This has the effect of limiting the dimensionality of the collective atomic state to N + 1 as only states with the same angular momentum are coupled.In particular, the J = N/2 subspace is referred to as the Dicke manifold [10].The Hilbert space of the entire system is then spanned by the basis {|n ⊗ |J, M }, where |n is a Fock state of the cavity mode defined by â † â |n = n |n .
In the experiments of [29], the dominant source of dissipation is cavity loss.We model this with the standard quantum optical master equation in Lindblad form [28,29] where κ is the cavity field decay rate and ρ is the reduced density operator of the atom-cavity system.The Hamiltonian (1) reduces to the usual, balanced Dicke Hamiltonian when λ − = λ + and reduces to the Tavis-Cummings Hamiltonian when λ + = 0.

A. Semiclassical model
In the thermodynamic limit of N → ∞, we consider the time-evolution of the expectation values of the operators â, Ĵ− , and Ĵz , scaled as where γ ∈ R because Ĵz is a Hermitian operator.We use ( 1) and ( 7) to derive a system of nonlinear differential equations for α, β, and γ, given by These equations are the semiclassical approximation of the unbalanced, open Dicke model.They form a fivedimensional dynamical system whose behavior we wish to explore.
In the case of maximal angular momentum, the model exhibits spin conservation of the form which confines the β-and γ-components of trajectories onto the Bloch sphere.Hence, the conservation law (10) eliminates a degree of freedom.Therefore, the system is four-dimensional with state space R 2 × S 2 , which is a four-dimensional hypercylinder.
In the semiclassical regime, the parity symmetry manifests itself as a Z 2 -equivariance of system (9), given by the linear transformation In particular, this symmetry produces trajectories with a symmetric counterpart on the Bloch sphere.However, symmetry-related trajectories represent the same physically observable states, since the observables |α| 2 , |β| 2 , and γ are invariant under T Z2 .
In addition to the parity symmetry (12), the unbalanced, open Dicke model also features a parameter exchange symmetry, such that the exchange leaves the dynamics invariant subject to a reversal of the roles of the North and South poles of the Bloch sphere.Hence, the dynamics that occur near the South pole take place on the North pole under the parameter exchange (13).The physical interpretation of the parameter exchange symmetry is that the rotating frame about which the co-and counter-rotating terms oscillate has changed direction.Thus, the "co-rotating terms" are now rotating against the frame and vice versa for the counter-rotating terms.Then for the case ω 0 < 0, ω > 0 the meaning of the co-and counter-rotating terms is reversed.Conversely, the associated exchange of λ − and λ + also implies a change in sign of the effective atomic frequency ω 0 .

C. Continuation and Stereographic Transformation
All computations are performed with the numerical continuation software package auto-07p [33,34], which is used to find and continue equilibria, periodic orbits, and their bifurcations as system parameters are varied.
However, the spherical part of the state space R 2 × S 2 of system (9) can introduce issues when attempting to continue solutions in the presence of conserved quantities.To remedy this, we transform the Bloch sphere to a plane by stereographic projection, such that the conservation law (10) is always implicitly satisfied.
More specifically, we map every point (β, γ) on the Bloch sphere S 2 , except for the chosen projection point, to the plane; hence, the transformed variables can take on any value in R 2 .
With the real-valued variables we choose the "North pole," of the Bloch sphere to be the projection point, and define the stereographic map Its inverse is The equations of motion for the transformed system in the variables {a 1 , a 2 , x, y} are obtained by multiplication of the Jacobian matrix of the coordinate change g by the vector containing the equations of motion for {a 1 , a 2 , b 1 , b 2 , γ}, then substituting the expressions for {b 1 , b 2 , γ} by {x, y} (for more information see [35]).The resulting equations of motion for the stereographically projected system on R 4 are System (20) is the central object of study of this paper; all bifurcation analysis is performed with these equations.Images of the Bloch sphere are obtained with the inverse transformation g −1 .Note that the Z 2 -equivariance now manifests itself as the symmetry transformation which is the point reflection in the origin of R 4 .

III. PHYSICAL ATTRACTORS
A key result of the experiments performed by Zhiqiang et al. [29] is the existence of three phases of behavior: normal, superradiant, and oscillatory.The normal phase features a zero photon number in the steady-state, â † â → 0 as t → ∞, and trivial atomic states, either |N/2, N/2 for spin-up or |N/2, −N/2 for spin-down.
That is, all atoms are either excited or in their ground states.The superradiant phase, on the other hand, features two asymmetric steady-states, with non-zero photon number.The oscillatory phase is identified by persistent oscillations in the photon number about an unstable superradiant state.
In the semiclassical regime, the normal phase is characterized by a stable equilibrium point with zero photon number, whereas the superradiant phase corresponds to a pair of stable equilibria with non-zero photon number.The oscillatory phase is marked by the existence of a pair of stable periodic orbits oscillating about the two unstable superradiant equilibria.Thus, we classify equilibria of the system as normal phase equilibria or superradiant equilibria.Figure 2 shows a selection of characteristic trajectories on the Bloch sphere (left column) and temporal traces of the photon number (right column) for each of the three cases.In Figure 2(a1), we show a single trajectory from the North pole to the South pole (i.e., γ → −1/2 as t → ∞) in the normal phase, during which the cavity emits a pulse of photons, shown in Fig. 2(a2).In the superradiant phase, the symmetry of the steadystate is broken and the steady-state on the Bloch sphere is located off of the poles as shown in panel (b1).Here the photon number, shown in panel (b2), tends to a non-zero value.In the oscillatory phase, a pair of periodic orbits develop on the Bloch sphere, one of which is shown in panel (c1), and the photon number oscillates about a superradiant state, shown in panel (c2).

IV. STEADY-STATE BIFURCATIONS
We first consider the effect of the generalization to unbalanced coupling on the superradiant phase transition.Here, we let the coupling strengths λ ± vary while keeping κ, ω, and ω 0 fixed.This phase transition marks the 1.5 onset of significant cooperative interactions between the ensemble of atoms and the radiation field.In this section, we discuss the superradiant phase transition in terms of the steady-state bifurcations of the semiclassical model (20).
We show in Fig. 3  In Fig. 3(a), we show the bifurcation diagram for λ − = 1.The bifurcations of the South pole and superradiant equilibria are given in Fig. 3(a1), and the bifurcations of the North pole are shown in Fig. 3(a2).The system starts, for λ + small, on the stable branch N S 0 which then undergoes a supercritical pitchfork bifurcation P spr to become unstable along the branch N S 1 .Here, two branches of stable superradiant equilibria SR 0 emerge.Each of them then becomes unstable at (a pair of) Hopf bifurcations H, and regains stability at the same time as the North pole equilibrium point at λ + ≈ 2.236; see Fig. 3(a2).This transition marks a special point, labelled F, where the North pole becomes stable.This transition is discussed in more detail in Sect.VII.Slightly before this transition, at λ + ≈ 1.932, the unstable branch N S 1 undergoes a subcritical pitchfork bifurcation P sub leading to the creation of two unstable branches SR 1 , which then meet the stable branches SR 0 and vanish in the saddlenode bifurcations SN + .
When λ − = 2, shown in Fig. 3(b), the transition to superradiance emerges in two stages.Again starting for λ + small, the system starts along the (stable) branch N S 0 .As λ + increases, a pair of saddle-node bifurcations SN − create two superradiant branches each, the unstable branch SR 1 and the stable branch SR 0 .The unstable branch SR 1 moves down to the South pole until it collides with the branch N S 0 in a subcritical pitchfork bifurcation P sub , which destroys the superradiant equilibria and turns the South pole equilibrium unstable along branch N S 1 .The following steady-state bifurcations as λ + increases feature the same structure as those in Fig. 3(a).
In physical terms, the key difference between the two cases, λ − = 1 and λ − = 2, is the manner in which the superradiant phase comes into existence.Fig. 3(a) describes the superradiant phase transition when λ − = 1.Here, the superradiant phase emerges with the creation of two stable superradiant equilibria, after the supercritical pitchfork bifurcation P spr ; this is the semiclassical description of the classic Dicke phase transition as in the case of balanced coupling λ − = λ + .However, the transition described by Fig. 3(b) consists of two stages, where the emergence of the superradiant phase with definite large amplitude and the disappearance of the normal phase now occur separately.In between the two stages, the normal and superradiant phases coexist as a multistable configuration, the structure of which is discussed in Sect.V.

Superradiant phase diagram
Figure 4 shows the superradiant phase diagram in the (λ − , λ + )-parameter plane, which outlines regions where the number of equilibria is the same.It is constructed by determining the locations of the pitchfork and saddlenode bifurcations.A necessary condition for the pitchfork and saddle-node bifurcations to occur is that the determinant of the Jacobian matrix of system (20), evaluated at the respective equilibrium point, is zero, i.e., det(J) = 0, where J is the Jacobian matrix of (20).Evaluating at the South pole equilibrium point (α, β, γ) = (0, 0, −1/2) yields the quartic equation which is quadratic in λ 2 + .Its (positive) solution is which we identify as the pitchfork bifurcation curve in Fig. 4 when taken as a function of λ − .In particular, this curve gives the well known value λ * of the superradiant phase transition for the case of balanced coupling λ * = λ − = λ + [28], namely The asymptotic expression for the pitchfork curve in the limit of large coupling is λ + = λ − ± √ ωω 0 , which we note is independent of the cavity dissipation rate κ.Following a similar procedure, but evaluating the Jacobian matrix J at the superradiant equilibrium points, we find another quartic equation that we can identify with the locations of the saddle-node bifurcations: The solutions to this equation give the two saddle-node bifurcation curves in the (λ − , λ + )-parameter plane as the lines The pitchfork bifurcation changes criticality when the saddle-node and pitchfork bifurcation curves intersect to form codimension-two pitchfork-saddle-node bifurcations (PSN).These are, in fact, identical to the "dissipationinduced tricritical points" studied by Soriente et al. in [30].However, Soriente et al. use a different parameterization of the model (see Appendix); in our version, the PSN points are Multistability and a single tricritical point in the unbalanced model have also been noticed in a theoretical study of Keeling, Bhaseen, and Simons [36] in a different parameter regime.
In Fig. 4, there are two disconnected saddle-node curves SN ± , defined for λ − ≥ λ PSN − , corresponding to either the plus or minus sign in (26).The pitchfork bifurcation curve ( 23) is cut by the points PSN into three distinct sections where the pitchfork bifurcation is either subcritical or supercritical.
We also see that the steady-state bifurcations are symmetric about the diagonal line λ − = λ + , i.e., invariant under the parameter exchange λ − ↔ λ + .This is due to the fact that ( 22) and ( 25) contain only even powers of λ ± .Furthermore, in the limit of large coupling, the pitchfork bifurcation curves become parallel to the diagonal, λ − = λ + .Hence, the system cannot exhibit multistability phenomena between normal and superradiant phases in the case of balanced coupling, because the diagonal intersects the pitchfork curve only once, at the point λ * given in (24).

V. MULTISTABILITY
As we have seen in Sect.IV, the generalization to unbalanced coupling in the open Dicke model allows for a case of the superradiant phase transition where the emergence of the superradiant phase and the disappearance of the normal phase occur separately.In between the two stages of this transition, the two phases coexist due to the existence of multiple stable equilibria of the system, which are "selected" as the final state under time evolution by the initial conditions.

Basins of attraction
The nature of the multistability can be examined by considering the basins of attraction of the normal and superradiant phase equilibria.That is, the sets of initial conditions that converge to particular stable equilibria as t → ∞.In the multistable case, there exist six equilibria in total: the stable normal phase equilibrium at the South pole N S 0 , the unstable normal phase equilibrium at the North pole N N 2 , two stable superradiant equilibria SR 0 , and two saddle superradiant equilibria SR The basins of attraction of SR 0 and N S 0 are separated by the stable manifolds W s (SR 1 ) of the saddle superradiant equilibria, which are the set of points in phase space that converge to these equilibria as t → ∞.These equilibria have one eigenvalue of the Jacobian matrix of system (20) with positive real part, and three with negative real part.Thus, by the Stable Manifold Theorem [37], the stable manifolds are three-dimensional.To illustrate the basin, we consider initial conditions for α = 0, i.e., an initially empty cavity.In this case the boundaries of the basins of attraction are given by the onedimensional intersection of the three-dimensional stable manifolds and the two-dimensional surface defined by α = 0.The boundaries of the basins of attraction can then be found numerically by computing the manifold W s (SR 0 ) with numerical continuation techniques [38].
In Fig. 5, the basins of attraction are shown on the Bloch sphere for a parameter set in the region of multistability.The two basins of attraction of the superradiant equilibria SR 0 appear as teardrop-shaped lobes emanating from the North pole.The region outside of the lobes is the basin of attraction for the normal phase equilibrium N S 0 .As λ + is increased, the superradiant basins of attraction grow from the North pole.They continue to grow until they reach a point where the system is "halfway" between the superradiant and normal phases; that is, the basins of attraction occupy equal areas.The superradiant basins then continue to grow until they entirely envelop the Bloch sphere at the second pitchfork bifurcations, after which the system is in the superradiant phase.Persistent oscillations in the semiclassical model arise due to a pair of Hopf bifurcations, where a stable periodic orbit bifurcates from a superradiant equilibrium point, after which the superradiant equilibrium point becomes unstable.After the Hopf bifurcation, two asymmetric periodic orbits Γ emerge after the transition due to Z 2 -symmetry.Figure 6 shows this on the Bloch sphere.Before the bifurcation, in Fig. 6(a), the system is in the superradiant phase and trajectories converge to the superradiant equilibria SR 0 .After the Hopf bifurcation, shown in Fig. 6(b), the stable periodic orbit Γ emerges from the (now unstable) superradiant equilibrium point; this corresponds to the oscillatory phase as observed experimentally in [29].Note that a self-intersection on the Bloch sphere in Fig. 6 is not a violation of the uniqueness theorem as the Bloch sphere is a projection of the state space.
The bifurcation diagram for λ − = 2 in Fig. 7(a) shows the branch of the bifurcating periodic orbits as repre-0 1 max(a 1 ) The inset of Fig. 8 gives details of two additional bifurcations that arise due to oscillations in system (20).The first is a codimension-two Bogdanov-Takens bifurcation (BT) where the Hopf bifurcation curve ends on the saddle-node curve SN + [37,39].Here, the Jacobian matrix J features two zero-eigenvalues and satisfies det(J) = Tr(J) = 0.The second codimension-two point, which we label FP, is the first intersection of the pole-flip transition curve F with the pitchfork curve P spr .

VII. POLE-FLIP TRANSITION
Figure 8 shows the pole-flip transition curve F that forms the boundary of the spin-down normal phase (γ → −1/2), where the South pole of the Bloch sphere is stable, and the spin-up normal phase (γ → 1/2), where the North pole is stable.This transition therefore describes quite a dramatic change to the atomic dynamics, as the steady-state spontaneously changes from all atoms occupying their ground states, to all atoms becoming excited.
Along the pole-flip transition curve F, the Bloch sphere is foliated by an infinite number of periodic orbits, see Fig. 9.The North pole equilibrium point also features two eigenvalues with zero real part.To solve for the curve F, we use the bialternate product [39].It has the property that for a matrix A, the matrix A 2I has eigenvalues that are pairwise sums of the eigenvalues of A [40], where I is the identity matrix.Parameter values where the Jacobian matrix J has eigenvalues with zero real part can then be found when the matrix J 2I has a zero eigenvalue, i.e., solving det(J 2I) = 0. Evaluating this determinant at the North pole equilibrium point gives us the analytical expression for the curve F Similarly, the determinant could be used to compute Hopf bifurcations of the superradiant equilibria; however, this requires analytic expressions for these equilibria, which are quite complicated and not particularly enlightening.The parameter values for the points where the pitchfork bifurcation and pole-flip transition occur simultaneously, on the other hand, are readily found from ( 28) and ( 23) as There are two cases corresponding to two different foliations of the Bloch sphere.The first case consists of periodic orbits that all oscillate about the poles, which is illustrated in Fig. 9(a).This is the boundary case of the two configurations of the normal phase, spinup or spin-down, and is located along the curve F for λ − < min(λ FP1 − , λ FP2 − ).We refer to this case as the normal pole-flip transition.The second case occurs after the pitchfork bifurcation creates a pair of superradiant equilibria.Here oscillations develop around these equilibria, confined by a pair of orbits homoclinic to the South pole, forming a figure eight, as illustrated in Fig. 9(b).This transition, which we refer to as the superradiant pole-flip transition, occurs along the curve F when At the point of the pole-flip transition, there is a twodimensional surface of periodic orbits which features a continuous symmetry, and thus has an associated conserved quantity.We find that this conserved quantity is the semiclassical energy per atom, given by E = Ĥ /N , that is, ) is conserved along the two-dimensional surface foliated by periodic orbits.This property was confirmed by careful numerical investigation by calculating the energy E along a selection of periodic orbits computed via numerical continuation.
Moreover, this surface of periodic orbits is attracting, meaning orbits with initial conditions that lie outside of the surface will converge to a periodic orbit with constant energy.Hence, generically, orbits will initialize outside the surface of periodic orbits and during their convergence their energy will change until the energy approaches the constant energy of a particular periodic orbit, not necessarily with the same initial energy.
The presence of energy-conserving periodic orbits in an open system such as this is surprising.These periodic orbits represent the balance of processes that preferentially drive the spins up or down, where there is also a balance of the energy exchange between the atoms and the cavity field, even in the presence of dissipation.However, as shown in a recent article [31], this behavior is modified significantly with the addition of spin dissipation and pumping.

VIII. EMERGENCE OF CHAOTIC DYNAMICS
We now turn our attention to describing and analyzing the bifurcations of periodic orbits, in particular, those that lead to the onset of chaos.We observe perioddoubling cascades leading to the formation of two chaotic attractors, related by mirror symmetry, followed by their collision to form a single chaotic attractor through a sequence of global bifurcations.

A. Period-doubling cascades
In a period-doubling bifurcation, a periodic orbit bifurcates to form a new periodic orbit with twice the period of the original, effectively modulating the amplitude of oscillation with a period of two.
A very common route to chaos is through a sequence or cascade of period-doubling bifurcations leading to the creation of chaotic attractors [37].Note that due to the system's Z 2 -symmetry, most periodic orbits and chaotic attractors have symmetric counterparts; hence, perioddoubling bifurcations also occur in pairs.
A period-doubling route to chaos, shown in Fig. 10, can be observed in system (20).Here, we show trajectories for three values of λ + (with λ − held constant) on the Bloch sphere, temporal traces and Fourier spectra log |F{|α| 2 }| of the photon number (where F{•} denotes the Fourier transform).The system is in the oscillatory phase in Fig. 10(a): here the two periodic orbits (only one of which is shown) feature a single loop [Fig.10(a1)] and the photon number |α| 2 temporal trace is shown in Fig. 10(a2).Furthermore, the Fourier spectrum has a single prominent peak at the frequency of oscillation [Fig.10(a3)], and a much smaller first harmonic.Upon increase of λ + (as in Fig. 10(b)), a period-doubling bifurcation has occurred; the original periodic orbit on the Bloch sphere is now unstable (and not shown) and a periodic orbit with two loops is now the attractor [Fig.10(b1)].The temporal trace of the photon number shows the amplitude modulation of the oscillation cycle [Fig.10(b2)], and the Fourier spectrum shows the emergence of peaks in between the previous peaks [Fig.10(b3)], corresponding to additional contributions from other oscillation frequencies.As λ + is increased further (as in Fig. 10(c)), the periodic orbits bifurcate many more times as the sys-tem undergoes an infinite sequence of period-doubling bifurcations and two chaotic attractors emerge on the Bloch sphere (Fig. 10(c1) shows one of them).The temporal trace of the photon number shows a more complicated, non-repeating oscillation cycle [Fig.10(c2)], and the Fourier spectrum is now broad with many peaks due to contributions from infinitely many oscillation frequencies [Fig.10(c3)], which is characteristic of chaotic motion [41].Figure 11 shows a one-parameter bifurcation diagram for λ − = 2 and varying λ + ; here, we represent periodic solutions by the maximum of the real part of the cavity field amplitude in panels (a) and (b).We observe a cascade of period-doubling bifurcations (PD) bifurcating off of the primary branch of periodic orbits Γ.Every period-doubled branch of periodic orbits then has another period-doubling bifurcation along it, where another branch of periodic orbits emerges.The accumulation point of the period-doubling cascade is approximated near PD 6 , which occurs at λ + ≈ 2.624.

B. Global bifurcations to chaos
Numerical continuation of the periodic orbits emanating from Hopf bifurcations (see Fig. 7) reveals an interesting spiral structure of the primary branch of periodic orbits Γ, shown in Fig. 11(a), with an enlarged version in Fig. 11(b).The turning points (with respect to λ + ) of the spiral are saddle-node bifurcations of periodic orbits, denoted SNP [37], where a pair of stable and unstable periodic orbits collide and vanish, or are created.In the case observed in Fig. 11(b), there are sequences of saddlenode bifurcations of periodic orbits accumulating onto a bifurcation point Hom.This oscillating behavior is more apparent in Fig. 11(c), where the bifurcation diagram of the period T is shown.Here, saddle-node bifurcations of periodic orbits SNP occur at the turning points or folds of the curve.Notice that the period T of the periodic orbit increases without bound, and accumulates to a value of λ + where system (9) exhibits a homoclinic bifurcation Hom.Here, system (9) exhibits a pair of trajectories that converge forward and backward in time to a pair of saddle equilibria, forming two homoclinic orbits, in this case, to the superradiant equilibria SR 1 .The following eigenvalue inequality at the moment of the homoclinic bifurcation Hom is satisfied where v s ≈ −0.469 ± 2.729i and v u ≈ 4.318 are the (stable) complex and (unstable) leading eigenvalues of the superradiant equilibria.Hence, the homoclinic bifurcation point (Hom) that is exhibited in system ( 9) corresponds to a wild Shilnikov bifurcation [42]; and it is responsible for the existence of infinitely many saddle-node bifurcations and the oscillatory behavior observed in Fig. 11.The Shilnikov theorem [43,44] then asserts the existence of infinitely many m-homoclinic bifurcations, where the associated homoclinic orbit makes m loops around the primary homoclinic orbit.Figure 12(a1) shows the Shilnikov homoclinic orbit in (a 1 , x, y)-space, as computed by numerical continuation of the primary periodic orbits up to a sufficiently high period.The temporal trace of a 1 along the orbit is shown in Fig. 12(a2); note the Shilnikov orbit's characteristic tail of decaying oscillations.
Additionally, we show in Fig. 12(b-c) 2-homoclinic and 4-homoclinic orbits, respectively, out of the infinitely many that exist nearby the point Hom.These homoclinic orbits are computed, also by continuation, as perioddoubled orbits of sufficiently high period.The temporal traces of a 1 show the formation of additional peaks corresponding to each loop of the homoclinic orbits at the moment of the corresponding 2-homoclinic and 4homoclinic bifurcation, shown in Fig. 12(b2) and (c2).
Figure 13(a) illustrates in the (a 1 , y)-plane the Shilnikov homoclinic orbits that were shown in Fig. 12(a).Here, we show both asymmetric homoclinic orbits; these are homoclinic to superradiant equilibria, so are labelled Hom SR   1   with the subscript referring to the number of their loops.In addition to Shilnikov homoclinic orbits to asymmetric equilibria, we also find (at different parameter values) a pair of homoclinic orbits to the South pole of the Bloch sphere N S 1 , which is the origin of the stereographically projected system (20); these homoclinic orbits are denoted Hom S i and Fig. 13(b) shows the pair Hom S 5 that make five loops and connect to the saddle equilibrium N S 1 on the symmetry subspace of the parity transformation (21).Furthermore, N S 1 has real eigenvalues and, hence, the pair of homoclinic orbits feature exponentially decaying tails without oscillations.
Nearby (in parameter space) the homoclinic bifurcation to the South pole equilibrium N S 1 , we find an interesting phenomenon related to the two chaotic attractors formed via period-doubling cascades and isolated in phase space, shown before in Fig. 10(c).We observe in numerical simulations the formation of a single, larger chaotic attractor as parameters are varied.Figure 14(a1) shows the two isolated chaotic attractors in (a 1 , x, y)space, with the temporal trace of the photon number given in Fig. 14 FIG.13.Two types of homoclinic orbits found in the system (20).Panel (a) shows a pair of homoclinic orbits to the two (saddle-focus) superradiant equilibria SR1 for (λ−, λ+) = (5.75,7.05).Panel (b) shows a pair of homoclinic orbits to the (real saddle) normal phase equilibrium N S 1 for (λ−, λ+) = (0.8, 1.627).Here κ = ω = ω0 = 1.
tors then merge into a single chaotic attractor, illustrated in Fig. 14(b1).This transition connects two previously isolated regions of phase space, for positive and negative a 1 , and produces a dramatic change in the dynamics of the system.The most obvious change is noticeable in the temporal trace of the photon number, shown in Fig. 14(b2), which now displays large jumps corresponding to when the state switches from one region of phase space to the other.This difference in the photon number dynamics should, in principle, be distinguishable experimentally, provided the timescale on which the experiments take place can be resolved.The Fourier spectra before and after the collision of chaotic attractors are given in Fig. 14(a3-b3), respectively.Notice that there is no appreciable difference in the spectrum before and after the collision.
The collision of asymmetric chaotic attractors shown in Fig. 14 must occur through interaction with the symmetric subspace of the parity transformation (21).We find that this global transition involves homoclinic orbits of the symmetric equilibrium N S 1 , that is, the South pole of the Bloch sphere.Each panel of Fig. 15 shows the one-dimensional unstable manifold W u (N S 1 ) of N S 1 .When the two chaotic attractors are separated as in Fig. 15(a), each side of the manifold W u (N S 1 ) accumulates onto each respective localized chaotic attractor.Figure 15(b) shows the homoclinic orbit Hom S 5 , which exists for λ + ≈ 1.627, connects the two regions of phase space.
Collision of a pair of asymmetric chaotic attractors to form a single symmetric chaotic attractor.For λ+ = 2.21 in (a1), there are two isolated chaotic attractors, shown in light and dark blue.The temporal trace of the photon number is given in (a2), and the Fourier spectrum of the photon number is shown in (a3).For λ+ = 2.26, the merged, symmetric chaotic attractor is shown in (b1), with the temporal trace of the photon number in (b2), and the Fourier spectrum in (b3).
In Fig. 15(c), where λ + has been increased further, the unstable manifold W u (N S 1 ) accumulates onto the merged chaotic attractor; note in particular, that each side of the unstable manifold W u (N S 1 ) accumulates on the entire, non-localized attractor.
As λ + is increased further, the merged chaotic attractor contracts [see Fig. 15(d)] as it approaches the pair of homoclinic orbits associated with the superradiant poleflip transition F, shown in Fig. 15(e); compare with the image on the Bloch sphere in Fig. 9(b).At this transition the merged chaotic attractor disappears.When viewed for decreasing λ + , we conclude that this flip transition generates a large symmetric chaotic attractor and the associated chaotic temporal trace of the photon number.

C. Overall Phase diagram
A two-parameter bifurcation diagram in the (λ − , λ + )plane outlining all of the dynamics discussed thus far is presented in Fig. 16.Specifically, the phase diagram shown in Fig. 16 shows the locations of saddle-node bifurcations SN ± , super-and subcritical pitchfork bifurcations P spr/sub , and Hopf bifurcations H; compare with Fig. 8.It also shows the curves of period doubling bifurcations PD i , and homoclinic bifurcations to superradiant equilibria Hom SR i , which are found with numerical con-  tinuation; note that we use the same notation Hom i to refer to the homoclinic orbits and their associated bifurcation curves.The different bifurcation curves divide the (λ − , λ + )-parameter plane into several distinct regions of different dynamics.These regions are the spin-up or spin-down normal phases N ↑/↓ where the photon number tends to zero, the superradiant phase SR where the atoms collectively emit into the cavity mode, multistable regions N ↑/↓ + SR where there are simultaneous normal and superradiant phases, and the oscillatory phase Osc where the photon number features simple persistent oscillations.Moreover, there is a region Chaos which is more complicated.It is bounded by the accumulation of the curves PD i and the curve F. Inside the region Chaos we find localized and non-localized chaotic attractors.The transition between the two cases is indicated by graded shading, with the lighter shade indicating two isolated chaotic attractors and darker shade indicating a merged chaotic attractor.As we discussed, the transition between the two cases involves the curves Hom S i near the curve P sub .Figure 16 also features a number of codimension-two points that act as organizing centers for the creation of codimension-one bifurcation curves; they will be discussed later in the section.

Step-through the phase diagram
To illustrate typical examples of the dynamics as λ − and λ + are varied, we follow a relevant closed path through the two-parameter bifurcation diagram Fig. 16 as we traverse across the parameter plane, beginning in the bottom right-hand region of the spin-down normal phase N ↓ ; see [45] for an interactive web version of the bifurcation diagram.
Consider the path illustrated in Fig. 17 A number of codimension-two points are part of the phase diagram in Fig. 16. Figure 18 is an enlargement that shows these near the intersection of the curves P sub and F, with a further enlargement as an inset.Codimension-two points act as organizing centers for the codimension-one bifurcation curves of the phase diagram.Apart from the points PSN, FP, and BT already seen in Fig. 8, we now also show in Fig. 18 the points FP c and BV, as well as the associated bifurcation curves P sub , F, and the homoclinic bifurcation curves Hom SR i , whose associated orbits are homoclinic to SR 1 , and Hom S , whose orbits are homoclinic to N S 1 .The point FP c in Fig. 16 and Fig. 18 is the main organizing center of global bifurcations and chaotic dynamics.Here, the (subcritical) pitchfork bifurcation curve P sub crosses the pole-flip transition curve F. This point is responsible for the creation 0.4 0.9 1 of the period-doubling bifurcation curves PD i and of infinitely many homoclinic bifurcations to both superradiant and normal phase equilibria.Out of the Bogdanov-Takens point BT emerges a single (superradiant) homoclinic bifurcation curve Hom SR which reaches the organizing center FP c .Coming out of the point FP c we expect an infinite number of m-homoclinic bifurcations to superradiant equilibria Hom SR m and to the South pole of the Bloch sphere, one of which (the same homoclinic orbit as shown in Fig. 13(b)) is shown in the inset of Fig. 18 as Hom S .
The inset of Fig. 18 also shows the curve labelled CC where the South pole equilibrium N S 1 transitions from having real to having complex conjugate eigenvalues.At the intersection point BV of the curves Hom S and CC, the homoclinic bifurcation transitions from a homoclinic bifurcation to a real saddle to a Shilnikov bifurcation (to a saddle focus with associated decaying oscillations).More specifically, the point BV is a Belyakov transition point (see [42] for more information), with the eigenvalue condition where v s = −0.5 and v u ≈ 1.0985 are the stable and unstable leading eigenvalues of the South pole equilibrium point.Theory [42] tells us that there are infinitely many Shilnikov bifurcation curves to the South pole equilibrium emanating from the Belyakov point BV.We expect that this transition is an accumulation of m-homoclinic orbits Hom S m emerging from the organizing center FP c .One such curve, for the homoclinic orbits shown in Fig. 15(b), is shown in the inset of Fig. 18 (the Shilnikov bifurcation curves emerging on the right-hand side of the CC curve are not shown).This transition should occur for each homoclinic bifurcation curve emerging from the organizing center FP c , so we expect infinitely many Belyakov points BV: one on each homoclinic curve as it crosses CC.

D. Effect of detuning
So far we have considered the case ω 0 = ω of zero detuning between the effective cavity frequency ω and the effective atomic frequency ω 0 .As we show now, this is not as restrictive as it initially may seem.Namely, the bifurcations we have studied here are generic in the sense that, for sufficiently small changes of the parameters κ, ω, and ω 0 , the phase diagram in Fig. 16 remains qualitatively the same.We now address the natural question of how much of a detuning between ω and ω 0 is still allowed before the phase diagram changes qualitatively.Note that in the experiments of [29] ω and ω 0 depend on the sum and difference, respectively, of the frequencies of the driving lasers, so either can be changed quite easily in an experiment.To be definite, we consider changing the effective atomic frequency ω 0 , while holding the effective cavity frequency ω constant.
We find that the relative positions of organizing centers (codimension-two points) of the phase diagram in Fig. 16 do not change qualitatively when ω 0 is decreased over the entire interval 0 < ω 0 ≤ ω.More specifically, it follows from the expression (28) for its slope that, as ω 0 approaches zero, the pole-flip transition curve F approaches the diagonal in the (λ − , λ + )-plane where λ − = λ + .Similarly, according to expression (23), the curve of pitchfork bifurcations approaches this diagonal as well as ω 0 → 0. In other words, the curves F and P collapse onto the diagonal in this limit.Nevertheless, we find that F and P do not change their relative positions until the limit ω 0 = 0 is reached, and the same is true for the curves emerging from their intersection point, the organizing center FP c .
When ω 0 is increased from ω 0 = ω the relative positions of the organizing centers of phase diagram in Fig. 16 do not change qualitatively over the interval ω < ω 0 ≤ ω * 0 , where ω * 0 = √ ω 2 + κ 2 .At this critical frequency ω * 0 there is a qualitative change of the phase diagram because then the expression (28) for the slope of pole-flip transition curve F becomes which, according to (26), is exactly the (ω 0 -independent) slope of the saddle-node curve SN + .Thus, at the critical frequency ω * 0 the curve SN + coincides with F, which corresponds to the codimension-three situation where the points FP, FP c , PSN, and BT all coincide as a unique point of tangency between the pitchfork bifurcation curve P sub and the pole-flip transition curve F. For ω 0 > ω * 0 the relative position of the points FP and FP c along the curve P sub has changed.This means that these two codimension-two points effectively move through each other at ω 0 = ω * 0 , leading to a qualitative change of the phase diagram.A detailed analysis of the phase diagram in the (λ − , λ + )-plane for ω 0 > ω * 0 is beyond the scope of this work, but we conjecture that it does not change qualitatively until the limit ω 0 → ∞ is reached, when the pole-flip transition curve F also reaches the diagonal where λ − = λ + .
We conclude that the phase diagram as shown in Fig. 16 is representative of the dynamics if 0 < ω 0 < ω * 0 , which constitutes a considerable range of detunings.Because of the parameter exchange symmetry (13), the same is true when −ω * 0 < ω 0 < 0, subject to an exchange of the coupling strengths λ − ↔ λ + ; this also implies a reversal of the roles of the North and South pole of the Bloch sphere, i.e., the pitchfork bifurcation to superradiance takes place at the North pole in this case.
We finally remark that the critical frequency ω * 0 is characterized by a specific relationship between the cavity decay rate κ and the effective atomic and cavity frequencies ω 0 and ω.That is, if one sees the light field amplitude α and the atomic polarization β as isolated oscillators (i.e., zero coupling) then ω * 0 separates the situation where the frequency of the atoms is faster or slower, respectively, than the combined decay and frequency of the light field.It is an interesting observation that this relation has consequences for the overall organization of the phase diagram of the unbalanced, open Dicke model, namely by generating a codimension-three bifurcation of the nonlinear dynamics.

IX. CONCLUSION AND OUTLOOK
We have analyzed phase transitions in the unbalanced, open Dicke model in the semiclassical regime with a dynamical systems approach.In doing so, we examined how unbalanced coupling leads to a splitting of the emergence of superradiance and the destruction of the normal phase in the semiclassical description, giving rise to the phase coexistence studied in [30].We have also investigated the emergence of superradiant oscillations in terms of Hopf bifurcations that the system undergoes.
During the analysis, we discovered a transition between the two types of normal phase configurations, spin-up and spin-down, which we have called the pole-flip transition.We found that at the point of transition, an infinite number of energy-conserving periodic orbits foliate the Bloch sphere.This can occur in two ways: the "normal" pole-flip transition where the periodic orbits oscillate about the poles, or the "superradiant" case where periodic orbits develop also about superradiant equilibria.Explicit expressions for the pole-flip transition and the conserved energy were determined.
We have also demonstrated the onset of chaos in the system arising from period-doubling cascades of periodic orbits created at the Hopf bifurcations.These create two symmetry-related localized chaotic attractors, which can then collide via a sequence of homoclinic bifurcations to form a single symmetric and non-localized chaotic attractor as parameters are varied.We also showed that this type of symmetric chaotic attractor can arise when the foliation of periodic orbits that arise from the superradiant pole-flip transition is perturbed, which constitutes a completely different mechanism for the creation of chaotic behavior.
The different regions of possible behaviours and transition curves between them have been represented as a phase diagram in the (λ − , λ + )-plane, which clearly identifies where superradiance, multistability, oscillations, and chaos can be found.We have shown how the overall phase diagram is organized by a number of codimensiontwo bifurcation points, from where different transition curves meet and/or emerge.In fact, the point FP c , where the curves of symmetry breaking and the pole-flip transition cross, acts as the main organizing center that generates the chaotic dynamics that can be observed in the system.We showed that the relative positions of the different codimension-two points do not change when the effective cavity and atomic frequencies are detuned over a very large range, up to a critical value that we determined explicitly.This implies that the phase diagram presented here is valid well beyond the case of zero detuning for which it was computed.
The generalization of unbalanced coupling, which may seem like a simple generalization, reveals such rich and complicated structure of additional dynamics.It is quite surprising given that considering unbalanced coupling does not add any new, explicitly nonlinear terms to the Dicke Hamiltonian; it simply splits the contribution of the interaction term of the standard Dicke model to two independently varying terms.Indeed, the dynamics of the standard Dicke model can be found along the diagonal of the (λ − , λ + )-plane, and all the additional behaviours we found require a quantifiable amount of unbalancing.
There are clearly some promising directions for further investigation of the semiclassical approximation of the unbalanced, open Dicke model.First of all, ongoing work is devoted to a more detailed study of the bifurcation sequence involved in the collision of two localized chaotic attractors to form a non-localized chaotic attractor.Another task is to determine what the phase diagram looks like beyond the critical value of the detuning between the effective cavity and atomic frequencies.We have already presented some ingredients in terms of the relative positions of organizing centers, but a proper answer will require a substantial bifurcation analysis in the spirit of the one presented here.Finally, it is an interest-ing mathematical challenge to study the exact nature of the main organizing center FP c , which involves the existence of a two-dimensional surface that connects both poles of the Bloch sphere, and which is foliated by periodic orbits and a pair of homoclinic orbits.Since such a foliation is rather special, in the sense that there is a conserved quantity, an important related question is how this situation, which occurs along an entire line in the (λ − , λ + )-plane, can be unfolded.This will require the consideration of additional decay channels, such as (either individual or collective) dissipation of the spins as considered in [31,46].In the semiclassical regime, the addition of atomic decay means that angular momentum is no longer conserved so the atomic dynamics are freed from the Bloch sphere, adding an additional degree of freedom to the dynamics.
Our study has also shown that regions in the (λ − , λ + )plane corresponding to the different behaviors described in this work are sufficiently large that it should be possible to reach them experimentally.This will obviously require the implementation of sufficiently unbalanced coupling strengths, among other experimental challenges.The possibility of future experiments naturally raises many questions regarding the correspondence of the semiclassical results presented here with the fully quantum mechanical model for finite N .Specifically, the region of multistability, where there is coexistence of the normal and superradiant phases, could provide an avenue for studying the correspondence of quantum to semiclassical dynamics since basins of different attractors are typically not well defined in a quantum mechanical context.In addition, any oscillations in the quantum model will typically be washed out of the master equation description by its statistical nature.A possible approach to analyzing oscillations and chaos in the fully quantum mechanical model is to utilize a quantum trajectory method [47][48][49].In particular, our demonstration that its semiclassical approximation of the unbalanced, open Dicke model possesses rich chaotic dynamics implies that studying its quantum mechanical counterpart could provide an interesting context in which to study the onset of quantum chaos.

FIG. 1 .
FIG. 1. Schematic of the unbalanced, open Dicke model as proposed by Dimer et al.[28] in a cavity QED system.An optical cavity with a resonant frequency ω and decay rate κ is coupled to N atoms with frequency splitting ω0.The atoms are also driven by two independent lasers, which produce two independent Raman coupling strengths λ±.Each atom is a four-level system with two stable ground states |↓ and |↑ .The laser and cavity fields are sufficiently detuned from the atomic transitions that the two excited states |+ and |− can be adiabatically eliminated in the theoretical model, producing an effective two-level system with comparable rotating and counter-rotating atom-cavity interaction terms.

FIG. 3 .
FIG. 3. Bifurcation diagrams for the two scenarios of the superradiant phase transition.Panel (a) shows the one-stage onset of superradiance, and (b) shows the two-stage case.Solid curves indicate stable equilibria, dashed curves unstable equilibria and saddles.For (a1-a2) λ− = 1, for (b1-b2) λ− = 2; other parameters are κ = ω = ω0 = 1.Normal phase equilibria are labelled N, and superscripts indicate whether the equilibrium is the North or South pole with N or S. Superradiant equilibria are labelled SR.Subscripts indicate the number of unstable directions.
the equilibria of system(20) in a bifurcation diagram for λ − = 1 in panels (a) and λ − = 2 in panels (b).Stable equilibria are shown as solid curves, whilst unstable equilibria are shown dashed.Normal phase equilibria are labelled N. The normal phase al-ways has two equilibria which are distinguished with a superscript N for the North pole or S for the South pole.Superradiant equilibria are labelled SR.Subscripts indicate the number of unstable directions.

FIG. 4 .
FIG. 4. Two-parameter bifurcation diagram describing the superradiant phase transition in λ− and λ+.Shown are the locations of the saddle-node bifurcations SN (red curves) and the locations of the pitchfork bifurcations P (blue curve), light-blue indicates the pitchfork bifurcation is subcritical and darker blue indicates supercritical.Codimension-two pitchfork-saddle-node bifurcations are shown as black dots and labelled PSN.The vertical gray lines indicate the slices of the parameter plane that produce Figs.3(a1) and (b1), respectively.Other parameters are set at κ = ω = ω0 = 1.

FIG. 8 .
FIG. 8.The two-parameter bifurcation diagram describes the oscillatory phase transition in λ− and λ+.Hopf bifurcations are located on the black curve H.The oscillatory region is bounded by the curves H and F; compare with Fig. 4.Here κ = ω = ω0 = 1.

FIG. 11 .
FIG. 11.One-parameter bifurcation diagram (a)showing the primary (black curve) and period-doubled (orange curves) branches of periodic orbits; the inset shows the period T of the first period-doubled orbit.Panel (b) shows an enlargement of the primary branch as a homoclinic bifurcation denoted Hom is approached.The spiral structure of the primary branch is shown in panel (c) in terms of the period T of the periodic orbit.Here κ = ω = ω0 = 1, λ− = 2.

FIG. 16 .
FIG. 16.Phase diagram in the (λ−, λ+)-plane of the unbalanced, open Dicke model for κ = ω = ω0 = 1, showing the bifurcation curves from Fig. 8 and a selection of bifurcation curves associated with chaotic behavior.Color indicates regions of different behavior: the two normal phase regions, spin-down (white) and spin-up (light gray), are denoted N ↓ and N ↑ , respectively.The superradiant region (blue) is denoted SR, the multistable regions (red) are labelled SR + N ↓ and SR + N ↑ , the oscillatory region (gray) is denoted Osc, and the chaotic region (gradient orange to purple) is denoted Chaos, and its shading indicates the transition from isolated chaotic attractors (orange) near the period-doubling curves PDi to the merged chaotic attractor (purple) near the poleflip transition curve F.
: from N ↓ it enters the superradiant region SR either directly or via the multistable region SR + N ↓ .The path then moves counter-clockwise through Osc and Chaos to N ↑ and back to N ↓ .The associated sequence of behavior and bifurcations is as follows: N ↓ The system is in the normal phase, with all spins down.N ↓ Pspr − −− → SR Moving across the supercritical pitchfork bifurcation P spr the system undergoes the single-stage transition into superradiance.N ↓ SN− − −− → SR + N ↓ Alternatively, moving across the P sub − −− → SR saddle-node bifurcation SN, the system becomes multistable, with coexistent spin-down normal and superradiant phases, and then becomes entirely superradiant when the subcritical pitchfork bifurcation P sub is crossed.SR H − → Osc As the system moves across the Hopf bifurcation curve H, the superradiant equilibria bifurcate into stable periodic orbits, and the system enters the oscillatory phase; see Fig. 6.Osc PDi −−→ Chaos The system undergoes a perioddoubling cascade as it moves across the period-doubling curves PD i , forming a pair of asymmetric, localized chaotic attractors; see Fig. 10.Chaos Inside of this region the pair of chaotic attractors collide via a sequence of homoclinic bifurcations to the origin Hom S i to form a merged chaotic attractor, indicated qualitatively in Fig. 16 by shading.Chaos F − → SR + N ↑ Transitioning across the curve F, the system undergoes the superradiant pole-flip transition, where the chaotic attractor disappears by contracting onto the pair of homoclinic orbits at F; compare with Fig. 15(d-e).As a result, the system is now multistable and simultaneously in the spin-up normal and superradiant phases.SR + N ↑ SN+ − −− → N ↑ At the saddle-node bifurcation curve SN + , the superradiant equilibria disappear and the system remains in the spin-up normal phase.N ↑ F − → N ↓ Moving again across the pole-flip transition curve F, the normal phase transitions from spin-up to spindown; see Fig. 9(a).

FIG. 18 .
FIG.18.An enlargement of Fig.16accentuating the codimension-two bifurcations.The inset shows the curve (CC) that delineates regions where the equilibrium N S 1 has real or complex conjugate (green shaded region) eigenvalues.Also shown in the inset is a homoclinic bifurcation curve to the South pole equilibrium point (Hom S ).