Dicke time crystals in driven-dissipative quantum many-body systems

<jats:title>Abstract</jats:title>                <jats:p>The Dicke model—a paradigmatic example of superradiance in quantum optics—describes an ensemble of atoms which are collectively coupled to a leaky cavity mode. As a result of the cooperative nature of these interactions, the system’s dynamics is captured by the behavior of a single mean-field, collective spin. In this mean-field limit, it has recently been shown that the interplay between photon losses and periodic driving of light–matter coupling can lead to time-crystalline-like behavior of the collective spin (Gong <jats:italic>et al</jats:italic> 2018 <jats:italic>Phys. Rev. Lett.</jats:italic>                    <jats:bold>120</jats:bold> 040404). In this work, we investigate whether such a Dicke time crystal (TC) is stable to perturbations that explicitly break the mean-field solvability of the conventional Dicke model. In particular, we consider the addition of short-range interactions between the atoms which breaks the collective coupling and leads to complex many-body dynamics. In this context, the interplay between periodic driving, dissipation and interactions yields a rich set of dynamical responses, including long-lived and metastable Dicke-TCs, where losses can cool down the many-body heating resulting from the continuous pump of energy from the periodic drive. Specifically, when the additional short-range interactions are ferromagnetic, we observe time crystalline behavior at non-perturbative values of the coupling strength, suggesting the possible existence of stable dynamical order in a driven-dissipative quantum many-body system. These findings illustrate the rich nature of novel dynamical responses with many-body character in quantum optics platforms.</jats:p>

section_content: The study of emergent dynamical phenomena in interacting quantum many-body systems constitutes a frontier of research in modern quantum optics and condensed matter physics. In this quest for phases of quantum matter without equilibrium counterpart, time crystals (TC) represent a promising candidate for a novel form of dynamical order out-of-equilibrium. In TCs, observables dynamically entrain at a frequency subharmonic of the one imposed by an external periodic drive [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17], and they have been currently realised with trapped ions [18] and solid state systems [19] [20] [21]. In most previous studies, TCs are realised in closed interacting quantum manybody systems, which are prone to heating towards an infinite temperature state under the action of periodic drive [22, 23], therefore, a slowdown of energy absorption is customarily entailed via a disorder induced many-body localized phase [24] [25] [26] [27], or by fast driving [8, [28] [29] [30] [31] ]. 

An alternative pathway could consist in "cooling" time crystals via coupling to a cold bath, which can absorb the energy pumped by the periodic drive [31]. A natural candidate to explore this avenue is represented by a recent line of inquiry on the exploration of TC-like behaviour in the open Dicke model, which describes an ensemble of atoms collectively coupled to a leaky photon cavity mode. The periodic drive of Dicke light-matter interactions in the superradiant regime can entail sub-harmonic dynamical responses [1], however, the collective nature of interactions renders the dynamics of this class of TCs equivalent to a single body problem consisting of a mean-field p â (t) (â + â † )

FIG. 1. The driven-dissipative many-body Dicke model studied in this work. An ensemble of atoms are collectively coupled to a photon field â with a time-varying strength λ(t), which consistutes a Floquet driving that injects energy into the system. The atoms also interact with each other via a short-range interaction of strength J, which breaks the collective nature of the system. Dissipation of energy into a bath is included through the loss of photons at a rate κ. 

collective spin degree of freedom moving on the Bloch sphere. Our key goal is to understand the stability of the Dicke time crystal when one breaks the mean-field nature of the model. 

To this end, we explore the robustness of Dicke-TCs to local interactions which break the collective coupling of the original model (Fig. 1 ). We observe that this class of Dicke-TCs can remain stable to such mean-field breaking perturbations in certain limits. Crucially, this lifts the phenomenon from an inherently collective, meanfield effect to the steady-state behaviour of a dissipative many-body system. We note however, that unlike the traditional venue for discrete time crystals [9] [10] [11] [12] [13], where short-range interactions are essential for stabilizing time crystalline order, here, the short-range interactions are rather viewed as perturbations to the original mean-field
section_name: I. INTRODUCTION

section_content: The interplay between Floquet driving, dissipation, and interactions results in a rich set of dynamical responses. In particular, we find regimes where TCs are stabilised by the bath, which counteracts the energy pumped into the system by the drive. We also observe the emergence of metastable dissipative TCs, characterised by a slowly decaying envelope evolving eventually into a trivial steady state dominated by dissipation. In addition, we find a family of ferromagnetic driven-dissipative TCs with strong resilience to many-body heating.
section_name: Dicke-TC.

section_content: We consider an chain of N two-level atoms with shortrange interactions among each other

where ŝx,y,z = σx,y,z /2, and σx,y,z are Pauli matrices. 

The atoms are collectively coupled to a photon field, e. g., by placing them inside an optical cavity (Fig. 1 ), which can be described by the hamiltonian [32, 33]

where Ŝx,y,z = N i ŝx,y,z i. We allow the light-matter coupling to be varied in time, λ(t). Dissipation occurs when photons leak out of the cavity, as encoded by the quantum master equation

for the total density matrix of the system, ρ, where Ĥ = Ĥac + Ĥint, and κ characterises the rate of photon loss. 

When J = 0, the above reduces to the well-known open Dicke model [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47]. As the coupling is only between the single photon mode and the collective spin operator, Ŝx, the Dicke model is exactly solvable in the thermodynamic limit N → ∞: its dynamics can be described by the mean-field motion of the photonic amplitude, a, coupled to three classical degrees of freedom, Ŝx,y,z (t), evolving on the Bloch sphere. When J = 0, short-range atom-atom interactions break the exact solvability of H ac, spoiling the collective character of the Dicke hamiltonian. In addition to the collective mode S, which corresponds to the k = 0 Fourier mode s k ≡ N j=1 e -ikj s j, all other k = 0 modes could also be excited. Hence H int introduces quantum fluctuations in the spin (or atomic) degrees of freedom, which require treating the dynamics in Eq. (3) as a quantum manybody problem. 

We simultaneously account for dissipation and quantum fluctuations using a time-dependent spin-wave approach, which has been demonstrated effective in capturing dynamical quantum many-body effects [48] [49] [50]. Specifically, we first perform a time-dependent rotating frame transformation to align the time-dependent ẑ-axis with the collective Bloch vector S(t), and then in this comoving 'frame' we perform a Holstein-Primakoff transformation in order to expand the spin operators in H int to the lowest-order in the density of spin-wave excitations, (t) (see Appendix A for details). The many-body effects introduced by H int are encoded in the dynamical coupling between spin-waves, and the collective spin as well as photon field. Excitation of spin waves leads to a depletion of the k = 0 mode, 2| S(t)|/N = 1 -(t), similarly in spirit to approaches which incorporate selfconsistently the effect of quantum fluctuations in the dynamics of a "condensate" [51]. The spin-wave density, (t), representing the total population of all k = 0 spinwave excitations, is required to remain small at all times in order to have a self-consistent lowest-order Holstein-Primakoff expansion. By monitoring the growth of spinwave density, we can identify regions of "heating", where the collective spin order shrinks under effect of strong many-body interaction, accompanied by a large value of (t) in dynamics (see Appendix A).
section_name: II. THE MODEL

section_content: We impose a periodic modulation on the light-matter coupling λ(t): during a first "bright-time", nT ≤ t < (n + 1/2)T, we set λ(t) = λ = Ω, and during the "dark-time", (n + 1/2)T ≤ t < nT, we switch off λ(t). 

Here, Ω = 2π/T = (ω 0 + ω)/2 denotes the driving frequency. This periodic driving counteracts with the energy loss through photon leaking, and as previously found in Ref. [1] for the collective Dicke model ( Ĥint absent), it entails a period-doubling response in spin observables. The presence of such TC-like behaviour can be understood from the Z 2 symmetry of the Dicke model with constant λ, under the parity operator P = e iπ(â † â+ Ŝz+ N 2 ). For λ > λ c = 1 2 (ω 2 + κ 2 /4)ω 0 /ω, a quantum phase transition that breaks the Z 2 symmetry occurs, and the system enters a superradiant phase, featuring two steady states, with spin projection S x /N = ±X, S y = 0 and non-vanishing photon amplitude a = ∓ √ N λX/(ω-iκ/2)

c /λ 4 ). When δ ≡ (ω 0 -ω)/Ω = 0, the free evolution during the 'dark-time' accumulates a phase of π for both S x = S x and a, and the system switches from one steady state to another, i. e., S x → -S x and a → -a. As a result, the dynamics repeats after two cycles of the driving and a sub-harmonic response with period 2T appears. In such a picture, atoms are simply described as a single classical spin, which becomes invalid in the presence of Ĥint, and thus it deserves a careful investigation whether the Dicke-TC exists in a many-body system with nonzero J. 

We explore the spin dynamics for various interaction strength J and dissipation rate κ, using the timedependent spin-wave approach. Our analysis shows that the Dicke-TC order can exist beyond the collective case (J = 0) and survive many-body interactions. For a range of finite J/Ω 1 and 0 < κ/Ω 1, we observe a stable subharmonic response in S x (t), as plotted in Fig. 2 (a) for an instance, where the spin-wave density remains small and therefore it is robust to heating (see Appendix B). Upon increasing the values of J, the inelastic scattering induced by many-body interactions becomes efficient, provoking a sizeable growth of spin-wave density, which invalidates the lowest-order spin-wave expansion and makes the collective spin, S, crumble. In this regime of strong interactions, the system is prone to heating under the action of the Floquet driving (see also [52] ). However, as shown in Fig. 2 (b), for moderate strengths of J/Ω 1, such effect of heating can still be remediated with sufficient dissipation κ/Ω 1 (see also Appendix B and D). In this case, strong dissipation acts as a 'contractor' for the dynamics, guiding swiftly the system towards the desired non-equilibrium steady state and enables stable oscillations of S x. Meanwhile, the combination of strong dissipation and driving leads to remarkable disturbance of the spin state within each period, resulting in a larger value of S y compared to the Dicke-TC at small κ (cf. Figs. 2(a) and 2(b)). Here, the enhanced fluctuation of photon field associated with strong dissipation can induce large phase noise and destroy the Dicke-TC as well; however, for sufficiently large atom number N, the phase noise can still be suppressed, since the amplitude of photon field scales as ∼ √ N in the superradiant regime (see Appendix A), and thus allows the observation of the Dicke-TC. 

While period doubling is a fragile dynamical response in one-particle periodically driven systems, i. e. it disappears as a tiny δ = 0 is switched on, the collective (J = 0) Dicke-TC is robust in a range of small δ = 0, thanks to the macroscopic S x -order built during the superradiant "bright-time" and thanks to the "contractive" role of dissipation which guides the system towards the desired non-equilibrium steady state. Such robustness to deviations from the δ = 0 limit, persists upon inclusion of many-body interactions, J = 0 (as an example, δ = -0. 12 is used for Fig. 2 ). The persistence of many-body Dicke-TCs on time scales much longer than t κ ∼ 1/κ (see the exemplary dynamics in Fig. 2(a) ), indicates that they represent a long-lived phenomenon, since in the presence of dissipation, relaxation is typically expected to occur on time scales inversely proportional to the system-bath coupling, κ. Indeed, we never observe decay of the Dicke-TC order on the longest timescales accessible to our numerical study (e. g, see Fig. 2 (a) and footnote [53] ). The existence of Dicke-TC is also insenstive to initial conditions. As we checked in our numerical simulations, similar responses are observed over a wide range of initial conditions (see Appendix B). 

Here, dissipation is capable of stabilising rather than destroying the Dicke-TC order in S x (t) since it acts on the photon mode and therefore collectively on spin order. On the contrary, when local losses are introduced, dissipation is intrusive and detrimental, and it destroys TCs (see for instance Ref. [54] ). Note that in conventional discrete TCs, instead, the interactions serve to lock single-spin dynamics into a stable subharmonic response that is robust to perturbations or imperfections of the drive [13, 18]. 

When the rate of dissipation is intermediate, κ/Ω ∼ 1, we recognise a region of irregular (IR) dynamics, where the trajectory of S(t) is scattered on the Bloch sphere (Fig. 2(c) ). In this case, the photon amplitude is sizeably reduced, and since it contributes to building the S x -order via the light-matter coupling term ∝ λ, the system does not develop a S x component sufficiently strong in order to counteract the dephasing induced by the "transverse field" ∝ ω 0 Ŝz during the 'dark-time', and this results into a featureless dynamical response lacking of period-doubling. We also notice that in this regime a relatively small value of J can lead to a proliferation of spin-wave excitations, suggesting a tendency to heating (see Appendix D). 

Complementarily, with excessively large κ/Ω 1, while dissipation is sufficient to cool the system and prevents many-body heating, it also destroys Dicke-TC order: in this case, dissipation overdamps dynamics, and the collective spin of the system relaxes to a trivial steady state where all spins point down towards the south pole of the Bloch sphere [55] (see Fig. 2(d) ). This overdamped regime (OD) is thus characterised by a vanishing magnitude of the spin projection, |S x + iS y |. 

The different dynamical responses discussed above are summarised in a qualitative cartoon, Fig. 9, in Appendix D. It is interesting to note that, with increasing J, we find the system tends to be overdamped with a smaller κ, which is indicated by the downwards bending of the boundary between the OD and Dicke-TC regions shown in Fig. 9. Hence fast "cooling" does not always protect, but can instead destroy Dicke-TC behaviour. Nevertheless, dissipation plays a crucial role in establishing the rich phenomena in Fig. 1 (see Appendix E, where the impact of H int on TCs of the Lipkin-Meshkov-Glick model is addressed). We note that, in this work, the frequency of the drive, Ω, only plays the role of an overall energy scale.
section_name: III. DYNAMICAL RESPONSES

section_content: For intermediate values of both J and κ (see for instance Fig. 9 ), our system hosts another type of nontrivial behaviour: a dissipative metastable time crystal (MTC) characterised by a slowly decaying envelope, which deteriorates, in the long time, into a trivial asymptotic state dominated by dissipation with vanishing S x (see upper panel of Fig. 3 ). This behaviour is distinct from the Dicke-TC: in the lower panel of Fig. 3 we plot the associated spin-wave density (blue line), which exhibits a discontinuous jump at a non-vanishing value of J, when the lifetime of the TC starts also to decrease (purple line); this suggests that the metastability is not expected to manifest for small values of J, and therefore the conventional Dicke-TC represents a longlived phenomenon (see footnote [56] ). The lifetime, τ, of this MTC gradually decreases with J, following the empirical law, τ ∝ exp [-A(J/Ω) 1. 6 ] (with A a positive prefactor), and vanishes when the system enters the OD regime, where also the spin-wave density becomes small since the system reaches the fully polarised state in the negative ẑ-direction of the Bloch sphere (corresponding to spin-waves vacuum). The MTC found here appears FIG. 3. Top: The stroboscopic dynamics of the MTC for κ/Ω = 2. 5, J/Ω = 0. 08, δ = -0. 12. The dynamics of Sx(t) (blue, normalised by N ) appears indistinguishable from a conventional Dicke-TC response on time windows of the order of a few decades of cycles (inset); on longer time scales it displays instead a slowly decaying envelope. Bottom: The lifetime of the MTC (defined as the time when the stroboscopic amplitude of Sx decays to 0. 1), and spin-wave density of the MTC after 5 × 10 3 cycles. We find that the lifetime falls to zero following an empirical law ∝ exp [-A(J/Ω) 1. 6 ], with A a positive factor. 

as a genuine interplay of period driving, dissipation and interaction. We remark here that it doesn't result from a high-frequency expansion [28, 29], and thus is distinct from the "Prethermal" Floquet TCs found in previous studies [19, 31], since Ω is an overall energy scale in our system (see Fig. 1 ). A possible explanation for the phenomenon is suggested by the dynamics of (t): during the metastable evolution, the density of spin waves strongly fluctuates and is out of phase with the dynamics of the photon amplitude, accumulating at every cycle a tiny dephasing, which eventually leads S x (t) to collapse.
section_name: IV. METASTABLE DISSIPATIVE TIME CRYSTAL

section_content: As discussed above, the steady states underlying the Dicke-TC order possesses a ferromagnetic nature. When the many-body interactions H int is also ferromagnetic (J < 0), inter-spin interactions can reinforce the ordering along the ±x-direction, giving rise to robustness against heating and the overdamping caused by dissipation. Indeed, for non-perturbative values of |J| ∼ O(1), we find that a ferromagnetic (J < 0) dissipative Dicke-TC can Comparison between the stroboscopic dynamics of spin-wave densities of a ferromagnetic (cyan dashed line, J/Ω = -2. 4) and of an anti-ferromagnetic (blue line, J/Ω = 0. 2) Dicke-TC for δ = -0. 12 and κ/Ω = 0. 4. The orange area denotes the region of values of spin-wave density where the lowest order Holstein-Primakoff breaks and the system becomes prone to many-body heating (see footnote [53] ). 

be stabilised at intermediate dissipation rate without significantly heating up the system (see footnote [57] ). This is shown in Fig. 4, where we plot the spin-wave density for anti-(J > 0) and ferro-(J < 0) magnetic TCs: a sizeable spin-wave density denotes fragility to many-body interactions and a dynamics prone to heating, and such effects are expected to be pronounced at large J. Remarkably, this does not occur for ferromagnetic inter-spin interactions which develop tiny values of (t) even for |J| ∼ O(1). The emergence of such a Dicke-TC response within the coexistence of significant many-body interactions and dissipation rate, appears to us a strong incarnation of TC-like behaviour in drivendissipative platforms: it is a novel form of dynamical order out-of-equilibrium, which significantly departs both from the mean-field Dicke TC response (where J 0), and from conventional (many-body) discrete TCs where dissipation is not a constitutive ingredient (κ = 0). In this perspective, such ferromagnetic TC represents a nonequilibrium state of strongly coupled, driven-dissipative, quantum matter exhibiting rich dynamics that can trigger motivation towards the search of other non-trivial dynamical phases in many-body quantum optics.
section_name: V. FERROMAGNETIC DRIVEN-DISSIPATIVE TIME CRYSTAL

section_content: The Dicke model is currently engineered in several experimental platforms [58] [59] [60] [61] [62] [63]. We expect our results to be qualitatively insensitive to the details of the microscopic structure of the interaction term H int, and to hold in a broader set of models, and thus would be relevant for experiments where collectivity of the system is inevitably broken by inhomogeneous fields, or spatially varying light-matter couplings, or genuine interparticle interactions such as using Rydberg atoms [64] [65] [66] [67]. The stabilisation of Dicke-TC order seen for strong ferromagnetic spin-spin interactions can also be generalised to systems where the inter-spin interactions have a anti-ferromagnetic character (J > 0), given the capability to control atom-light coupling in cavity experiments [68]. For coupling ∝ λ(â + â † ) j (-1) j σx j, a similar Dicke-TC response exists in this case but with anti-ferromagnetic ordering. Hence a Ĥint with certain J > 0 would be expected to extend the Dicke-TC to a many-body regime. Another interesting possibility offered by the control of light-matter coupling consists in realising Dicke-TC responses with higher integer periods (nT with n > 2) without employing high spin atoms (see for instance Ref. [19] ). In the Appendix F we show that a Dicke model with coupling of the form ∝ λ(â + â † ) j [(-1) j/2 σ+ j + h. c. ] realises a dynamical response with quadruple period. 

The setup analysed in our work is closely related to quantum optics platform, but can also be relevant for condensed matter and solid state platforms, where the system can be modeled as a quantum spin chain coupled to a phonon bath. We believe that the outreach of our results has the potential to motivate a new generation of experiments on TCs in many-body systems, where the presence of a bath plays a crucial role at variance with current realisations [18, 19]. 

We use the shorthand O ≡ Ô for expectation values of operators Ô, and define the coordinates for the collective spin vector, S = N 2 (sin θ cos φ, sin θ sin φ, cos θ), in terms of the polar, θ, and azimuthal angle, φ, on the Bloch sphere, which allows writing compactly the equations of motion:

The effect of H int enters the dynamics of the collective spin vector via the spin-wave correlations δ αβ (α, β = q, p) and the spin-wave density, which are in turn dynamically coupled to the collective spin and photon field. 

Here,

where q k and p k are the canonically conjugated bosonic variables associated with spin-waves with wave-vector k = 0. Following the procedure in [50], we derive the equations of motion for the spin-wave correlations

These quantities intertwine with the equations of motion (A1), and represent the feedback of the non-equilibrium Gaussian fluctuations of spin waves on the motion of the FIG. 5. In the time-dependent spin wave theory, quantum fluctuations introduced by the short-range interaction term Hint are included as a self-generated bath which couples to the order parameter, S(t). The latter couples also to the photonic mode which is cooled by a zero temperature bath. Therefore, the dynamics of the order parameter results from the competing interactions with an internal spin wave bath (represented by spin waves' density, (t), in the sketch above), and with an external cold bath, mediated by the cavity photon, a(t). 

collective spin (the k = 0 mode). The self-consistent solution of Eqs. (A1) and (A6) yields the dynamics of the model. We notice that these equations of motion are derived in the thermodynamic limit [50], therefore we expect our results to be insensitive to the choice of N, as we checked for the main results of our work. 

In the last equation of Eqs. (A1), the phase fluctuation accompanied with photon loss has been neglected. This is valid in the large N limit, since the time crystals under scrutiny here exist in the superradiant regime, where the photon number is ∝ N, and thus photon noises become subleading. For small N, the effect of phase fluctuation can be incorporated in the spin-wave approach by adding a Langevin noise term in Eqs. A1(c). 

The spin wave density has to remain small during the course of the evolution, in order to render consistent the lowest order Holstein-Primakoff expansion, employed to derive Eqs. (A1) and (A6). This restriction of the spin-wave approach limits resolving possible structures in the dynamical responses when the heating from J is significant, which might be interesting to explore in future works. For J = 0, one can readily see from the above equations that (t) = 0 at any time, while for J = 0, the length of the collective spin is shrinked via 2| S(t)|/N = 1 -(t). A pictorial description of the method is provided in Fig. 5. 

FIG. 6. The maximum spin-wave density (t) within 5000 periods of dynamics (top) and the variance of Fourier spectrum (bottom) for various interaction strengths J at fixed dissipation rate κ corresponding to the Dicke-TC shown in Fig. 2(a) (blue dotes), and (b) (green squares) and δ = -0. 12. The variance is not calculated for those values of J corresponding to spin-wave density above the threshold (orange region). As a reference, the dotted line in the bottom panel shows the variance obtained for the IR dynamics in Fig. 2(c ). 

defined as χ = ν>0 dνp(ν) (ν -0. 5Ω)/ ν>0 dνp(ν) + ν<0 dνp(ν) (ν + 0. 5Ω)/ ν<0 dνp(ν), where p(ν) is the height of the Fourier spectrum at frequency ν. TCs feature stable subharmonic responses, with the spectrum composed of two sharp peaks located at ν = ±0. 5Ω, and thus a small value of χ. The variance for increasing J is plotted in the bottom panel of Fig. 6. With a small κ/Ω < 1, the spin-wave density quickly grows with J, which is reduced with a large κ/Ω > 1. The variance χ remains negligible for values of J with low spin-wave density, indicating the existence of the Dicke-TC with finite J. 

In Fig. 2 an initial state with θ = 0. 5π and φ = 0 has been used. We have also checked in our numerical simulations that the Dicke-TC is insensitive to initial conditions, thanks to the dissipative nature of dynamics: we observe similar sub-harmonic dynamics over a wide range of different initial states, as demonstrated in Fig. 7. analysis, except that at a small value of J it is interrupted by the emergence of other dynamical behaviours and heating: for increasing J, the boundary is set by a smaller, rather than a larger κ, suggesting that faster dissipation does not always protects the system against the heating from many-body interactions. As noted in Appendix A, limited by the choice of truncation in (t), a quantitative identification of the parameter regime for all dynamical responses and the nature of their boundaries is beyond the scope of this work. 

Appendix E: Comparison with integrability breaking of the periodically kicked 'LMG' model

In order to exemplify the non-trivial interplay that dissipation can have with many-body interactions, we have considered a similar analysis in a case governed by purely unitary dynamics. We have studied the Lipkin-Meshkov-Glick (LMG) model perturbed by H J,

The dynamics entailed by H is periodically perturbed by a collective rotation along the x-axis; the evolution operator reads in a period

with Ûkick ≡ exp [-iφS x ]. This protocol has been shown in Ref. [69] to display TC behaviour when φ = π (for instance, in the dynamics of the stroboscopic transverse magnetization). The TC in this case is robust to displacements around the φ = π point, i. e. for angles φ = π ± δ (with δ > 0). We have chosen this system as a comparison for the Dicke dynamics studied in the main text, since the latter effectively reduces to the LMG model via adiabatic elimination of the photon mode for large cavity detunings. 

The periodically driven unitary dynamics in (E2) does not entail a rich set of dynamical responses as in Fig. 9 : the TC persists (most likely in a pre-thermal fashion) for values of J smaller than a certain critical threshold, J c, above which the system develops sizeable spin-wave density and therefore crosses over into a regime of 'heating'. This suggests that the presence of dissipation enriches the dynamical responses of a periodically driven interacting quantum many-body system. 

Appendix F: Quadruple-period dynamical response

We consider the Hamiltonian

(F1) with the coupling strength between atoms and cavity controlled by tuning the phase factor ∝ k • x i ; such coupling can be realized by choosing proper laser dressing in a quantum optics setup. When e jk•xi = 1, the hamiltonian (F1) becomes the conventional Dicke model. In the case of e jk•xi = 1, the coupling strength varies from site to site. We choose e jk•xi = (-1) i/2, therefore, the hamiltonian (F1) becomes We divide the system into four sublattices, A, B, C and D, with mod[i, 4] = {0, 1, 2, 3}, respectively. Within each sublattice, atoms are collectively coupled to cavity photons. In the limit of N → ∞, the stationary solution of the corresponding master equation hosts multiple possible steady states, depending on the value Q = λ 2 ω/ω0 κ 2 /4+ω 2. The case Q < 1/4 corresponds to the normal state without superradiance. When 1/4 ≤ Q ≤ 1/2, the system is superradiant, with photon occupation

and similar to conventional Dicke superradiance, there are two possible spin states:

In this case the spin state at each site is in phase with the corresponding photon state. When Q > 1/2, the steady state consists of two groups. The first group is the same as in (F4), while the photon occupation in the second group is given by

This group includes 8 different spin configurations

with X = Q 2 -1/4/Q, and Z = √ 1 -X 2. These states correspond to having a "defect" in the spin configurations, and thus results in a lower photon number. A linear stability analysis suggests that the above states can all be stable. The existence of multiple steady states provides the possibility of producing subharmonic responses to external driving. 

When 1/4 ≤ Q ≤ 1/2, we can have a period-doubled dynamical response if we apply a Floquet driving scheme similar to the one discussed in the main text. 

When Q > 1/2, in addition to period doubling, we can have a dynamics with period 4T. However, this would require driving that can convert one steady state to another in Eq. (F6). A possible procedure is as follows: We initialise spins close to σ x A = σ y B = -σ x C = σ y D = 1, and let the atom-cavity system interact for some time to reach steady state, which can be diagnosed via monitoring the photon emission from cavity. Then, we start to apply a driving pulse θ at the end of each period T. θ consists of single-site rotations. Specifically, we apply Rz (π) rotation to all odd sites (A and C); even sites (B and D) are rotated depending on the measurement outcome of i∈A S x i prior to the pulse: Rz (π) for negative outcome, and Ry (π) for positive outcome. The pulse strength is kept equal to cavity detuning. 

If we measure i S y i for i in all even sites, we will observe a quadruple period in dynamics, as plotted in Fig. 10, demonstrating that the system indeed undergoes several steady states during the Floquet dynamics. Including small photon loss during θ, we still see oscillations at a stable 4T period. An alternative order parameter is the inversion i S z i in sublattices A and B, the dynamics of which also exhibits a period 4T (see again Fig. 10 ). 

The above analysis assumed that atoms within each sub-lattice remain in the collective manifold and that therefore there are no quantum correlations. When atom number is sufficiently large, this represents a reliable approximation; notice, however, that even for a small number of atoms, period doubling has been observed in Ref. [1] to persist for times longer than the decay time, therefore we can expect similar conclusions to hold for
section_name: VI. SUMMARY AND OUTLOOK

section_name: VII. ACKNOWLEDGEMENTS

section_content: The authors acknowledge discussions with D. Abanin, S. Gopalakrishnan, W. W. Ho, R. Moessner, C. Nayak, V. Oganesyan, M. Schleier-Smith, D. Sels and S. Yelin. The authors would particularly like to thank F. Machado for a careful reading of the draft and many helpful comments. JM acknowledges A. Lerose, A. Gambassi, A. Silva and B. Zunkovic for previous collaborations on the subject of time-dependent spin-wave expansions. This work is supported by the DARPA DRINQS program (award D18AC00014 ), the NSF ( PHY-1748958 ), the AFOSR-MURI Photonic Quantum Matter (award FA95501610323 ), the Harvard-MIT CUA, the National Science Foundation (NSF), the David and Lucille Packard Foundation, and the Vannevar Bush Faculty Fellowship. BZ is supported by the NSF through a grant for the Institute for Theoretical Atomic, Molecular, and Optical Physics at Harvard University and the Smithsonian Astrophysical Observatory. 

JM is supported by the European Union 's Framework Programme for Research and Innovation Horizon 2020 under the Marie Sklodowska-Curie Grant Agreement No. 745608 (' QUAKE4PRELIMAT '). 

the quadruple period dynamics discussed here.
section_name: VII. ACKNOWLEDGEMENTS

section_name: Appendix B: Persistence of Dicke-TC

section_content: Appendix A: Time-dependent spin-wave theory

In this Section we provide further information on the time-dependent spin-wave expansion, referring the reader to [48, 50] for a comprehensive discussion on the method.

section_content: In this section we provide more details of the Dicke-TC behaviour discussed in Sec. III. In the top panel of Fig. 6, we plot the maximum spin-wave density over the 5000 periods of dynamics obtained from the above timedependent spin-wave approach at stroboscopic times, when increasing the interaction strength J. As discussed in the main text, when the system is prone to manybody heating, sizable amount of spin-waves develops, and the Holstein-Primakoff expansion breaks down. Here, we choose as upper threshold for spin-wave density at the value, ∼ 0. 2. From the spin dynamics, we can also calculate the variance χ of the Fourier spectrum of S x +iS y,
section_name: Appendix B: Persistence of Dicke-TC

section_content: To understand the dynamical responses in our drivendissipative many-body system, here, we first apply a mean-field treatment to solve dynamics from the master equation Eq. ( 3 ). In Fig. 8 we display the various dynam- ical responses of the system. The boundaries between the IR and the Dicke-TC regions are dictated by a transition in the variance of the Fourier spectrum of S x (t), which is negligible in the case of the Dicke-TC. The OD regime is identified with a vanishing |S x + iS y | at 5000 periods, which exhibits a second-order phase transition when crossing into the Dicke-TC regime. In a mean-field analysis, quantum fluctuations are neglected and thus it does not predict the MTC, the heating region, and the enhanced robustness for ferromagnetic Dicke-TC, which are the genuine many-body results of our study. Instead, it shows as an artefact the persistence of Dicke-TC despite strong interactions.
section_name: Appendix C: Mean-field analysis of dynamical responses

section_content: To account for quantum many-body effects, we solve Eqs. (A1) to (A6) in the time-dependent spin-wave approach to obtain the spin dynamics. We explore over a range of parameters κ and J and find rich dynamical responses. Fig. 9 shows a qualitative sketch of regions where different dynamical behaviours are observed. We label the region with large spin-wave density (t) ≥ 0. 2 as heating (H), where the spin-wave treatment breaks down. IR regime is characterised with a large variance χ of the Fourier spectrum together with (t) < 0. 2. As discussed in Ref. [1] for the case of J = 0, chaotic dynamics may arise in the collective Dicke model, which can result in large numerical errors. Here, we also find that when κ/Ω ∼ 1 the numerical integration tends to be unstable, and we associate these instances with IR as well in Fig. 9. The MTC region is identified with a slow decay of S x to a nonzero value at 5000 periods, with a small χ and (t) remaining below 0. 2. The boundary between the OD and the Dicke-TC regime resembles the one in a mean-field
section_name: Appendix D: Summary of dynamical responses from time-dependent spin-wave analysis