Nonlinear dynamics of the dissipative anisotropic two-photon Dicke model

We study the semiclassical limit of the anisotropic two-photon Dicke model with a dissipative bosonic field and describe its rich nonlinear dynamics. Besides normal and 'superradiant'-like phases, the presence of localized fixed points reflects the spectral collapse of the closed-system Hamiltonian. Through Hopf bifurcations of superradiant and normal fixed points, limit cycles are formed in certain regions of parameters. We also identify a pole-flip transition induced by anisotropy and a region of chaotic dynamics, which appears from a cascade of period-doubling bifurcations. In the chaotic region, collision and fragmentation of symmetric attractors take place. Throughout the phase diagram we find several examples of phase coexistence, leading to the segmentation of phase space into distinct basins of attraction.

An interesting variation of the Dicke model considers a collective two-photon coupling and is motivated by recent progress in achieving strong and ultra-strong light-matter interactions.
In these regimes, multi-photon interaction processes (previously suppressed by the weak coupling strength) become more prominent [36][37][38][39][40][41][42]. Correspondingly, the two-photon Dicke model has gained considerable interest. Several studies highlight the collapse of discrete energy levels into a continuous band, occurring at a certain threshold coupling [43][44][45][46]. Furthermore, before the spectral collapse takes place, there is a phase transition to a 'superradiant'-like phase, where the fourfold discrete symmetry is spontaneously broken [47][48][49][50][51][52]. On the other hand, while stationary states of this model are relatively well understood, its dissipative nonlinear dynamics has not been investigated in detail.
Thanks to the collective nature of the coupling, Dicke models can be studied through the semiclassical approximation, valid in the limit of large N. For a regular single-photon interaction, nonlinear dynamics induced by counter-rotating terms leads to classical chaos in the strong-coupling regime. Away from the thermodynamic limit, the system undergoes a transition from quasi-integrability to quantum chaos, caused by the precursors of the quantum phase transition [53][54][55][56][57][58][59][60][61]. Significant efforts have been devoted to studying chaos in semiclassical and quantum regimes of closed Rabi and Dicke models [62][63][64][65][66]. Furthermore, including the effect of dissipative channels, nonlinear dynamics and chaotic behavior in a driven-dissipative setting [67] and considering anisotropic couplings [68] have been recently discussed. The rich nonlinear dynamics of the one-photon Dicke model suggests that similar interesting phenomena can be found in the two-photon model as well.
With these motivations in mind, we present in this article a study of nonlinear dynamics in the two-photon Dicke model, including the effects of bosonic field dissipation and anisotropic couplings (i.e., where the rotating and counter-rotating terms are unbalanced). After introducing the model in Sec. 2, together with a detailed justification of the mean-field approximation (see Appendix A), we discuss three main dynamical behaviors allowed by the system: (i) Various types of stable fixed-point, see Sec. 3; (ii) Limit cycles arising from Hopf bifurcations of normal and superradiant fixed points, see Sec. 4; (iii) Chaotic motion, see Sec. 5. The interplay of these dynamical regimes determines a complex behavior in parameter space, summarized by various phase diagrams. In particular, chaos emerges beyond a pole-flip transition point from a cascade of period-doubling bifurcations. We also find several types of phase coexistence, leading to fragmentation of phase space and sensitive dependence of the asymptotic dynamics on the initial state. The coexistence between different types of stable fixed points is discussed in more detail in Appendix B. Finally, Sec. 6 contains our concluding remarks.

The model
The Hamiltonian of the two-photon anisotropic Dicke model can be expressed as follows (setting = 1) whereâ (â † ) is the annihilation (creation) operator of the bosonic field with frequency ω 0 , σ (j) are the Pauli operators of two-level system j (qubit or atom), whileσ y ). The identical qubits have energy transition frequency ω q and two-photon interaction strength g with the bosonic field. The parameter λ (we assume λ > 0) models an imbalance of rotating and counter-rotating couplings, which may be altered by the intensity of the electric and magnetic fields in circuit QED implementations [39] or the power of lasers in trapped ions setups [37,38]. Differently from the one-photon case, the two-photon Dicke model features a four-fold symmetry with the generalized parity operator Π = (−1) N N j=1σ (j) z e iπâ †â /2 [43,48]. The Hamiltonian is invariant under the parity transformation (â → iâ,σ x,y → −σ x,y ). In terms of collective angular momentum operators Ĵ = 1 2 N j=1 σ (i) , we rewrite the Hamiltonian as: whereX =â 2 +â †2 andŶ = i(â 2 −â †2 ). Previous studies have shown that the two-photon Dicke model admits a 'superradiant'-like phase transition, after which the collective pseudospin acquires a macroscopic mean value. At the same time, the bosonic field is driven to a squeezed vacuum state. While the expectation value ofâ remains zero, squeezed quantum fluctuations lead to a phase with non-zero photon number, which is dubbed 'superradiant' [47]. Here we consider this model in an open environment by including a dissipation channel for the bosonic field, with decay rate κ. The system evolution is described by a standard Lindblad master equation: Note that the parity operator Π is still a symmetry of the system, i.e., Π †ρ Π is also a solution of Eq. (3). We will be mainly concerned with the large-N limit when spins can acquire a macroscopic population at coupling strengths comparable to the bosonic field frequency. This justifies a standard mean-field approximation, i.e., we decouple bosonic-qubit correlations as BQ ≃ B Q , obtaining the following set of nonlinear equations: In Eqs. (4-9) we rescaled the expectation values of the spin operators as s = 2 Ĵ /N, and defined the collective spin frequency ω z = Nω q . As only the bosonic field dissipation is taken into consideration, we can restrict the collective spin evolution to the unit sphere, giving a five-dimensional phase space. The above equations contain explicitly N, meaning that system dynamics will be influenced by the number of qubits. The influence will be most obvious on the decoherence time and non-stable dynamics. On the other hand the steady-state is independent of N.
One interesting point about Eqs. (4)(5)(6)(7)(8)(9) is that, at variance with s, the bosonic variables are not scaled with N and can remain small. In other words, the mean-field treatment is still valid for bosonic-field states which do not have a well-defined classical limit, as long as the atomic field has a macroscopic population. In Appendix A we show explicit comparisons to the exact evolution from Eq. (3), confirming the validity of Eqs. (4)(5)(6)(7)(8)(9) in the large N limit.

Stable fixed points
We first discuss the stationary states of the nonlinear dynamics, whose stability can be determined in a standard manner by the Jacobian matrix and Routh-Hurwitz criterion [69,70]. Various stable fixed points appear in our system, whose acronyms are summarized in Table 1. We will also discuss different types of coexistence phases and limit cycles which, for easier reference, are also listed in Table 1. The simplest type of fixed points features zero photon number and trivial spin states: NP ↓ (NP ↑ ) normal phase: all spins in | ↓ (| ↑ ) SP 'superradiant'-like phase, with macroscopic occupation of the atomic field U 0 localized phase: â †â → ∞, originating from a dissipative 'spectral collapse' B ↓ (B ↑ ) bistable phase: the stady-state, SP or NP ↓ (NP ↑ ), depends on initial conditions C ↓ (C ↑ ) U 0 and NP ↓ (NP ↑ ) coexist in this phase, similarly to B ↓ (B ↑ ) LC limit cycle The NP ↑ state is stable for while NP ↓ is only stable for λ < λ t . Therefore, the anisotropic parameter λ can induce a dramatic change of atomic dynamics, where all down-spin states | ↓ transform to up-spin state | ↑ (a pole-flip transition [68]). When λ < λ t there are two additional phase boundaries for the NP ↓ fixed point. The stability conditions are: where we defined: In other words, NP ↓ has an instability window for the intermediate range of couplings g t2 < g < g t3 . The phase boundary with the superradiant state of the closed model is recovered by setting κ = 0 in the above expression of g t2 . Notice also that g t3 → ∞ if λ = 1, i.e., the upper critical line only exists for the anisotropic model. By searching for nontrivial stationary conditions, the following 'superradiant'-like fixed points are found: in which atomic and bosonic fields are spontaneously occupied. By symmetry, a pair of states with opposite value of (± X , ± Ŷ , ±s x,y ) appear in the SP phase. A SP fixed point requires: Beyond Eq. (16) only normal-phase solutions are possible, thus we will always consider values of λ within this range. Furthermore, a physical SP solution only exists for: As it turns out, for λ ≤ 1 the SP fixed point is stable when g t4 < g < g t1 . For λ > 1 the stability condition is more complex, and we can only compute it numerically. Finally, we recall that in the closed system the continuous phase transition to the superradiant phase is followed at larger g by a 'spectral collapse' [43,47], as the coupling strength becomes comparable to the cavity frequency (g c = ω 0 /2 if λ = 1). At this collapse point, the discrete spectrum of H is transformed into a continuous band. This dynamical feature has not been destroyed by dissipation: We see from Eq. (15) that s z → 0 when g → g t1 , leading to a divergent photon number. Therefore, at g t1 the SP states evolve continuously to a pair of localized fixed points: For g > g t1 , initial conditions within the basin of attraction of such localized fixed point lead to rapid growth of photon number, reminding the collapse of numerous energy levels occurring in the closed system. The spectral collapse point of the closed system is recovered by setting κ = 0 in Eq. (14). This critical coupling is shifted by bosonic field dissipation to a larger value (e.g., g t1 = ω 2 0 + κ 2 /2 if λ = 1).

Phase diagram in the (ω z , g) plane
We now discuss explicit phase diagrams of stable fixed points in different parameter regimes. First, we take the qubit collective frequency ω z and the coupling strength g as variables, which for λ < 1 leads to a phase diagram with the structure of Fig. 1(a). Note that Eq. (15) and all the expressions for the phase boundaries (except for λ t ) do not depend on N explicitly if ω z = Nω q is kept constant. Furthermore, for large N also the boundary λ t approaches a well defined limit (λ t → 1). As shown in panels (b) and (c) of Fig. 1, N has an important effect on the timescale to reach the stationary state. However, the phase diagram of Fig. 1(a) and all stationary states are independent of N. In Fig. 1(a), a regular second-order phase transition between the NP ↓ and SP phases occurs at ω z < ω 0 /2, when g t2 and g t4 coincide. At g = g t2 the normal state becomes unstable and 'superradiant'-like solutions appear. This line is also identified as a pitchfork bifurcation point by bifurcation theory [69]. Examples of time evolution in the SP phase are shown in Fig. 1(b) and (c).
The situation is different at ω z > ω 0 /2, when the phase boundaries g t2 and g t4 do not coincide. This determines a bistable region g t4 < g < g t2 . We denote with B ↓ the portion where NP ↓ and SP are both stable, i.e., if the additional constrain g < g t1 is taken into account. As discussed already, the SP fixed point evolves to U 0 at g = g t1 . Therefore, when g t1 < g < g t2 the coexistence is between NP ↓ and U 0 , and is labeled as This phase diagram is applicable for λ ≤ 1. For λ > 1, additional regions with limit cycles appear and a pole-flip transition occurs at λ > λ t (see Sec. 3.2). The two dots mark the intersection of g t2 with g t4 and g t1 , located at ω z /ω 0 = 1/2 and 1, respectively. Right panels: system evolution in SP for g = 0.6, ω z = 0.2, and different values of N . Other parameters are ω 0 = κ = 1 and λ = 0.5. C ↓ . In such coexistence regimes, the final state is determined by the initial condition as each stable point has its own basin of attraction. We refer to Appendix B for a more detailed study of system evolution in the multi-stable regimes.
The last feature of the phase diagram in Fig. 1(a) is the reentrant transition of the normal phase, which occurs for g > g t3 and leads to a second pair of bistable regions B ↓ , C ↓ on the left side of the phase diagram. Since g t3 | λ=1 → ∞, these B ↓ , C ↓ regions will shrink and eventually disappear when approaching the isotropic limit.
From the above discussion we see that, by increasing the coupling strength g from zero to large values, the behavior of the system can be quite different depending on other system parameters. In particular, the critical lines g t1 and g t2 intersect at ω z = ω 0 . Therefore, in the range ω 0 /2 < ω z < ω 0 the system enters the various phases as At variance with the one-photon Rabi model (where a large detuning ω z /ω 0 is favorable to the formation of a superradiant phase), here the system does not support a 'superradiant'-like state when ω z → ∞. In this limit, the bistable region B ↓ shrinks to zero and a direct transition NP ↓ → C ↓ occurs. Such behavior would occur by considering N → ∞ at fixed ω q .

Phase diagram in the (λ, g) plane
To assess the influence of unbalanced rotating and counter-rotating couplings, we show in Fig. 2 phase diagrams of stable fixed points in the (λ, g) plane. We also compare N = 1 phase diagrams (left panels) to N = 10 (right panels), representative of the large-N limit. With respect to Fig. 1(a), an obvious difference is the appearance of the at which stable oscillation phases LC appear (only shown for λ < λ t ). The solid red thin line is g t1 while the other thin boundaries at λ < λ t follows the same line style of Fig. 1(a) for g t2 (dot-dashed), g t3 (dashed), and g t4 (dotted).
fixed point NP ↑ in the regime λ > λ t . The pole-flip transitions λ t are marked by vertical black lines in Fig. 2, dividing each phase diagram into two main parts. For small to moderate N, the position of this pole-flip transition has a sensitive dependence on ω z , which can be seen comparing panels (a1) and (b1). The dependence is non-monotonic, with the maximum λ t occurring at ω z = 2ω 0 . On the other hand, as seen from panels (a2) and (b2), the dependence of ω z is much weaker at large N, when λ t ≃ 1.
The λ < λ t region is quite complex, being determined by the critical lines g t1 , g t2 , g t3 , g t4 discussed already. Therefore, the same phases of Fig. 1(a) appear here. Like in Fig. 1(a), the phase boundaries are independent of N and there is a marked difference between the regimes ω z < ω 0 /2 (upper panels) and ω z > ω 0 /2 (lower panels). Instead, the region λ > λ t has a simpler structure dominated by a NP ↑ fixed point at g < g t1 and a coexistence region C ↑ for g > g t1 . When λ > λ t the fixed point NP ↑ is always stable, thus SP only appears in a nontrivial bistable region B ↑ which shrinks to zero at large N. As seen by a comparison of upper a lower panels of Fig. 2, the B ↑ region also rapidly shrinks with ω z , confirming that a large ω z is detrimental to the SP phase.
In concluding this section, we stress that the knowledge of stable fixed points is not sufficient to characterize the long-time dynamics. In panels (a1) and (b1) of Fig. 2 we also indicate the presence of limit cycles (LC), occurring for 1 ≤ λ < λ t . These limit cycles originate from a Hopf bifurcation of the SP fixed point, which becomes unstable in these regions. Therefore, in panel (a1) the limit cycle does not coexist with any stable fixed point (orange-red region), while in panel (b1) we find the coexistence of the limit cycle and NP ↓ (magenta region). When λ > λ t the system dynamics is more involved: Besides stable fixed points and limit cycles, we find the occurrence of chaos. For the moment, in the regime λ > λ t we have only shown the phases corresponding to stable fixed points. A detailed analysis of limit cycles is presented in the following Sec. 4, while the coexistence of stable fixed points, chaotic dynamics, and limit cycles at λ > λ t is discussed in Sec. 5. A summary of various transitions encountered in parameter space (including chaos and limit cycles) is presented in Table 2.

Stable oscillations
In non-linear dynamics, a Hopf bifurcation is a simple but important type of dynamic bifurcation. It describes the formation of a stable periodic orbit from a fixed point which, by varying system parameters, has lost its stability. In our system we find two Hopf bifurcations, marked by thick curves in Fig.2. The first Hopf bifurcation, g ′ t , is shown as a thick dot-dashed curve and originates from the SP fixed point. As shown in Fig. 3, for g < g ′ t near the transition line the SP fixed point is stable: Panel (a1) shows that no eigenvalue of the associated Jacobian matrix has a positive real part. On the other hand, a pair of conjugate eigenvalues have crossed the imaginary axis in panel (b1), where g > g ′ t . Before the bifurcation, as shown in panels (a2) and (a3), the system photon number and qubit trajectory converge to well-defined values, given by Eq. (15). After the bifurcation, the photon number in panel (b2) shows persistent oscillations and the limit cycle becomes especially obvious from the spin trajectory on the Bloch sphere, shown in panel (b3). Here, for clarity, we only plot the system trajectory for relatively large times, 350 < t < 400 (in units of ω −1 0 ), such that the system has already approached the stable periodic orbit. Due to symmetry, limit cycles bifurcating from SP appear in pairs. Only one of them is shown in Fig. 3.
It is also interesting to consider the case where the g ′ t Hopf bifurcation occurs in a bistable region. An example is the 'LC' region of Fig. 2(b1), where we find a coexistence of limit cycles, linked to the SP phase, and the stable fixed point NP ↓ . In Fig. 4 we give the example of two initial conditions belonging to different basins of attraction in phase space. We see in panel (a) how the stable fixed point NP ↓ and stable periodic orbit  coexist on the Bloch sphere: If the initial value is in the basin of the stable fixed point NP ↓ , the spin trajectory (blue line) converges to the | ↓ state, located on the bottom of the sphere. Instead, for another initial condition the trajectory evolves to a stable periodic orbit encircling the (unstable) SP fixed point. The difference in the asymptotic dynamics is also apparent in the time dependence of â †â , shown in Fig. 4(b). The second Hopf bifurcation occurs at λ t , i.e., coincides with the pole-flip transition and is indicated as thick black lines in Fig. 2. More precisely, for λ > λ t the NP ↓ becomes a limit cycle in phase space, which now coexists with a stable NP ↑ fixed point. Furthermore, we see in Fig. 2 that the g ′ t lines continue to the regime λ > λ t , where a complex coexistence of stable fixed points, different types of limit cycles, and chaotic dynamics occurs. These dynamical features are not shown in the phase diagrams of Fig. 2, which for λ > λ t only consider stable fixed points. While the interplay between different types of dynamics will be discussed in the following section, here we present in Fig. 5 a detailed example of the Hopf bifurcation of NP ↓ . In particular, panels (a1) and (b1) show the variation of eigenvalues near the bifurcation. At λ t , two conjugate eigenvalues cross the imaginary axis and acquire a positive real part. Before the bifurcation, as shown in panels (a2) and (a3), the system evolves to the stable fixed  Figure 6. Evolution from limit cycle to chaotic motion. In panels (a1)-(a3) we use g = 0.9 (taking ω 0 =1). The system evolves on a periodic orbit, forming a single loop in phase space. In panels (b1)-(b3) we use g = 1. After a period-doubling bifurcation, the system evolves on a periodic orbit with two loops in phase space. In panels (c1)-(c3) we use g = 1.1. A cascade of period-doubling bifurcations has lead to the formation of a chaotic attractor. |F | is the power spectral density of â †â . Other parameters: N = 10, ω z = 1.5, λ = 1.4. The initial state has s z = −0.99, s y = 0, and s x = 1 − s 2 z .
point NP ↓ with zero photon number. After the bifurcation, a small persistent oscillation in the photon number arises, see panel (b2). At the same time, a periodic orbit forms on the Bloch sphere near the fixed point NP ↓ , shown in panel (b3).

Chaotic dynamics
Besides the coexistence of multiple fixed points and limit cycles, chaos emerges in the regime λ > λ t . The transition to a chaotic trajectory is exemplified in Fig. 6, where for small coupling strength g the system is in a limit cycle originating from the fixed point SP. Stable oscillations with equal period and constant amplitude are found in the photon number, see panel (a2). At the same time, the qubit trajectory in panel (a1) is a periodic orbit featuring a single saddle loop. In panel (a3), we have computed the power spectral density (PSD) of â †â , which reflects the characteristics of system oscillations and is a common tool in the study of chaos. As seen, the PSD shows several discrete sidebands. At larger coupling strength, a period-doubling bifurcation occurs, seen in Fig. 6(b1). Now the qubit trajectory on the Bloch sphere bifurcates to a new periodic orbit with two saddle loops. The photon number dynamics in (b2) shows periodic oscillations with multiple amplitudes, and the photon number spectrum (b3) has more harmonic peaks. Still, the cyclic motion gives discrete lines in the PSD. By further increasing the coupling strength g, chaos eventually appears. The photon number dynamics of panel (c2) shows irregular oscillations with varying amplitude, giving rise to a continuous spectrum. In phase space, the qubit trajectory becomes a complex orbit formed by an infinite family of loops. This example shows that, in our system, chaos emerges from a cascade of period-doubling bifurcations, which is a typical route to chaos. The appearance of chaos can be explored in detail through a bifurcation diagram, see Fig. 7(a). There, we plot the evolution with g of the stationary points in the photon number oscillations. Periodic motion leads to a finite number of amplitudes, while chaotic motion has an infinite number of amplitudes. In the figure, we see that oscillations at small g have a single amplitude, which successively bifurcates to 2 n amplitudes (n = 1, 2 . . .). This process finally leads to a dense set of points in a certain range of g. The bifurcation diagram is in good agreement with the variation of the Lyapunov exponent λ LE with coupling strength, shown in Fig. 7(b). In the regime with a finite number of amplitudes λ LE is very close to zero (dashed line), consistent with periodic oscillations. In the regime with infinite amplitudes, positive values of λ LE confirm the occurrence of chaos.

Chaotic regions in the phase diagram
Accounting for the presence of chaotic dynamics, we can now give a more complete picture of the phase diagram at λ > λ t . In general, for a given choice of parameters, multiple types of motion coexist in phase space and the system initial values have a significant influence on dynamics. Determining such coexistence for general parameters  Table 1. Furthermore, we indicate with LC NP (LC SP ) limit cycles around the fixed point NP ↓ (SP). Periodic motion in the ultra-strong coupling regime of panel (a) is characterized by large oscillations covering both unstable fixed points, thus is labeled as LC. Chaotic motion occurs in the yellow areas. In panels (c) and (d), the narrow orange strips within the LC regions indicate a SP phase. We used ω 0 = 1, ω z = 1.5, and s x = 1 − s 2 z , s y = 0 for the initial state.
is a quite involved problem. Therefore, we choose to present the phase diagrams in Fig. 8 by fixing two representative initial states. We use different colors to denote various types of motion: The brown area corresponds to dynamics towards the localized point U 0 ; the blue area represents the evolution towards the stable fixed point NP ↑ ; the light blue area is a phase with periodic oscillations (LC), and the yellow area indicates chaotic motion. As seen, we may get a completely different phase diagram by changing the initial state. The coexistence of phases is implied by a comparison of corresponding phase diagrams, e.g., panels (c) and (d).
The left panels of Fig. 8 are for an initial state which is far away from the stable fixed point NP ↑ . Close to λ t we find at smaller values of g a region of periodic motion, formed by the fixed point NP ↓ through Hopf bifurcation. This area is labeled by LC NP . Instead, for λ λ t and larger g, the U 0 phase appears. The NP ↑ phase dominates the right part of the phase diagrams, with large values of λ. Periodic oscillations and chaotic motion appear between the areas of U 0 and NP ↑ . In this case, the periodic motion originates from the fixed point SP, thus we label it as LC SP in Fig. 8(c). In this panel, the two types of oscillatory phases are separated by stable dynamics, converging to the fixed point SP (orange area).  Fig. 8(c).
The number of qubits N has a marked influence on these phase diagrams. In panel (a), where N = 1, the chaotic regime occurs at relatively small values of g and large anisotropy λ. Instead, in panel (c) (where N = 10) we find that chaos is favored for a larger coupling strength. Besides, in panel (c) the phase U 0 appears in a much smaller range of parameters, while the area of fixed points NP ↑ is obviously larger than the area of localized fixed points U 0 and other types of dynamics. We also note that some periodic windows occur in the chaotic phase, which can be seen in the phase diagram (a) and (c) as blue dots within the yellow region. The occurrence of small windows of regular motion can also be seen in Fig. 7.
Finally, the right panels of Fig. 8 are for an initial state close to the stable fixed point NP ↑ . Now the blue area (where the system converges to NP ↑ ) occupies a larger portion of the phase diagram, while other phases have shrunk. Especially, periodic motion and chaos have almost disappeared in panels (b) and (d).

Collision of chaotic attractors
Similar to the one-photon model [68], we also find the interesting phenomenon shown in Fig. 9, where two chaotic attractors merge into a single large attractor. As shown in panel (c1), the large attractor can split again into two separated attractors. Usually, when considering symmetric initial conditions (±s xi , ±s yi , s zi ), the system evolves in separate parts of the phase space without making transitions between the two regions. This behavior is shown in Fig. 9(a1) where, after a long time evolution, the two trajectories form two spatially symmetric chaotic attractors. However, by increasing the coupling strength one may observe a dramatic change in system dynamics, shown in panels (b1) and (b2). Now the system makes frequent transitions between the two previously disconnected attractors, thus symmetric initial conditions give rise to an identical large attractor. In this regime, the evolution of the photon field has frequent jumps between different spaces with positive and negative quadrature components X , see panel (b2). The collision reveals a distinct route of formation for chaotic attractors. By further increasing the coupling strength, the large attractor splits again into two isolated small attractors, shown in Fig. 9(c1). The sequence of collision and fragmentation of attractors will continue at larger values of g, until the boundary of the chaotic region is reached.

Conclusion
In this article, we examined the nonlinear dynamics of the open two-photon quantum Dicke model in the mean-field approximation. As expected, the mean-field treatment becomes accurate in the limit of a large number N of qubits. By allowing unbalanced rotating and counter-rotating coupling strengths, and taking into account decay of the bosonic field, a rich dynamical behavior is found as a function of system parameters. Several phase diagrams are presented, showing a complex interplay of normal and 'superradiant'-like fixed points, as well as localized phases, limit cycles, and chaotic dynamics.
At variance with the one-photon Dicke model, here the 'superradiant'-like phase only displays macroscopic occupation of the atomic degrees of freedom, while the cavity remains in a Gaussian state with finite squeezing and number of photons. A divergence in photon number is found at a critical boundary corresponding to the 'spectral collapse' of the closed system. Furthermore, we find various coexistence phases where multiple dynamical behaviors are allowed. Here, the long-time dynamics reflects the segmentation of phase space into different basins of attraction. Especially interesting is the appearance of chaos from a cascade of period-doubling bifurcations. In this chaotic regime, we have highlighted a distinct mechanism of formation of the chaotic attractor through collision and splitting of symmetric lobes.
The two-photon Dicke model finds a natural realization in chains of trapped ions under bicromatic drive [8,37,43]. Previous proposals focused on λ = 1, but the anisotropic parameter can be controlled in a simple way through the relative strength of the two laser drives. Quantum superconducting circuits leading to a similar two-photon interaction, ∝ (â +â † ) 2σ x , have been discussed in [39,71], and it should be possible to construct alternative schemes realizing Eq. (1). For a given setup, terms neglected in the effective Hamiltonian might have an important effect in certain cases [72]. In particular, the full model is necessary to describe the long-time behavior in the U 0 phase.
While here we restricted ourselves to the semiclassical limit, it would be interesting to address how these nonlinear phenomena are reflected by quantum features beyond the mean-field approach. In particular, the occurrence of limit cycles and chaos might be probed through spectral properties of the stationary state and the Liouvillian [73][74][75][76][77], rather than directly from the time evolution. The exponential growth of the out-of-timeordered correlator (OTOC) is another interesting measure of quantum chaos, recently applied to the closed Dicke model [78]. Furthermore, we have supposed here that the dominant decoherence mechanism is from cavity decay. The role of qubit relaxation and dephasing would considerably enlarge the parameter space and possibly induce even richer dynamics, besides having practical relevance. It would be also interesting to explore possible links between the cyclic dynamics we have described and time-crystals in dissipative settings [79][80][81][82][83][84][85][86][87][88][89][90]. In general, our work provides insight into the rich nonlinear dynamics of the two-photon Dicke model, which may stimulate further investigation of classical and quantum chaos.  Figure A2.
Comparison of mean-field treatment and exact solution in the 'superradiant'-like phase. The curves with symbols in panel (a) and (c) are obtained from the master equation (3) while curves without symbols are the mean-field approximation. More precisely, in panel (a) and (c) we use the master equation to compute 4 Ĵ 2 x /N 2 and 2 â †âĴ z /N , respectively. In the right panels, the circles are stationary values calculated from the master equation. The solid line in panel (b) is a linear fit. The parameters are ω 0 = 1, g = 0.7, ω z = 0.2, and λ = 0.5.

Appendix A. Validity of mean-field treatment
The mean-field approximation is well-known to describe accurately the large-N limit of the one-photon Dicke model. In the two-photon model, however, the cavity does not need to approach a classical state. Here only the atomic variables have a macroscopic population while a non-zero photon number originates from squeezed vacuum fluctuations [47]. In particular, â †â does not scale extensively with system size.
Despite this difference, theĴ x,y,z operators of a large collective spin approach the classical limit, and this is sufficient to justify the mean-field approximation. In particular, we can approximately factorize correlations with the bosonic field (e.g., â †â J z ≃ â †â J z ), which is the basic assumption leading to Eqs. (4)(5)(6)(7)(8)(9). The purpose of this Appendix is to test explicitly this property, by comparison to the exact dynamics obtained from the Lindblad master equation Eq. (3).
We first consider in Fig. A1 the normal phase regime. In the left panels, (a) and (c), the full time evolution is shown, starting from the same eigenstate ofĴ x . We see that the initial agreement between master equation (ME) [91,92] and mean-field treatment (MF) holds for longer times when increasing N. The ME evolution (solid curves) yields a finite saturation value for â †â . However, panel (b) shows that â †â → 0 at large N, in agreement with the MF prediction. Similarly, panel (d) shows that s z → −1 in the stationary state of the ME.
The two initial states have s z = −0.69 and s z = −0.68 (while s x = 0) for the blue and purple lines, respectively. In panels (b1) and (b2) we choose g = 0.45, λ = 1.8, which are in the B ↑ (orange) area of Fig. 2(a1). The two initial states have s z = 0.09 and s z = 0.08 (while s y = 0) for the blue and purple lines, respectively.
. Therefore, the meanfield prediction Ĵ 2 x /N 2 ≃ s 2 x is accurate in the limit of large N. Panel (d) confirms that the decoupling of â †âĴ z become justified at large N, as the difference â †âĴ z M E − â †â Ĵ z M F ∼ O(1) is of subleading order if compared to Ĵ z ∼ O(N).
Even if corrections to mean-field follow the expected scaling with N, Figs. A1 and A2 show that they remain significant for relatively large system size. At finihe N, mean field-theory becomes very accurate when the SP states approach U 0 , due to the macroscopic occupation of the cavity [see Eq. (15), where s z → 0]. On the other hand, the presence of limit cycles and chaos can be only reflected by the transient quantum dynamics, as the master equation always has a well-defined stationary state. Thus, these regimes are challenging to identify from the quantum evolution with moderate N.

Appendix B. Coexistence of stable fixed points
As discussed in the main text, there are regions of parameters where two different types of stable fixed points coexist. The coexistence of (NP ↓ , SP) is indicated as B ↓ and occurs for λ < λ t . Instead, for λ > λ t a bistable phase B ↑ is found. In these phases, as illustrated in Fig. B1, the asymptotic dynamics is sensitive to the initial condition. In panel (a1) we show two trajectories on the Bloch sphere (left panels) which start at nearby points but eventually diverge from each other, to approach either the NP ↓ Figure B2. Basins of attraction on the Bloch sphere, taking the bosonic mode in the vacuum state at the initial time. In the left panel, except for the initial condition of the qubit, the parameters are as in the top panels of Fig. B1 (B ↓ phase). Yellow and orange areas are basins of attraction for SP, while the blue area is the basin of NP ↓ . Right panel: parameters are chosen as in the lower panels of Fig. B1 (B ↑ phase). Here the blue area is the basin of NP ↑ .
fixed point or the 'superradiant'-like one. Accordingly, the photon number in panel (a2) either approaches zero (blue line) or reaches a nonzero stable value (solid purple line). The dashed purple line is the steady-state photon number of SP, given by Eq. (15). Similar behavior can be found in the bistable regime B ↑ , whose dynamics of the qubit and photon number are respectively shown in panels (b1) and (b2).
As each stable fixed point has its basin of attraction, the phase space is divided into distinct regions: Initial states starting from a given basin will all converge to the same stable condition. For definiteness, in Fig. B2 we have chosen initial values â †â = X = Ŷ = 0 for the bosonic mode (i.e., the vacuum state) and N = 1. For this choice, in Fig. B2 we have represented graphically the basins of attraction on the Bloch sphere. While the behavior illustrated by Figs. B1 and B2 is generic, we see that the extension and shape of the basins of attraction is sensitive to system parameters. The area of the SP basin in the B ↓ example (left side of Fig. B2) is much larger than basin of SP in the B ↑ example (right panel). For the bistable regime B ↑ of Fig. B2, the phase space is primarily occupied by the NP ↑ basin of attraction.
For g > g t1 the SP phase turns into U 0 , therefore coexistence phases C ↑/↓ appear. The corresponding behavior for representative trajectories is illustrated in Fig. B3. The discussion is completely analogous to the B ↑/↓ phases, except that now the photon number shows a divergent time dependence for initial conditions in the U 0 basin of attraction. For other initial conditions, see the blue curves of Fig. B3, the system approaches rapidly the normal phase fixed point without displaying any singularity.