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

TL;DR

This paper explores the stability and emergence of Dicke time crystals in driven-dissipative quantum many-body systems, especially when short-range interactions break the mean-field approximation, revealing complex dynamical phases.

Contribution

It demonstrates the existence of long-lived and metastable Dicke time crystals in a driven-dissipative setting with short-range interactions, extending the understanding beyond mean-field models.

Findings

Time crystalline behavior persists with ferromagnetic interactions at strong coupling.

Dissipation can cool many-body heating, stabilizing time crystals.

Rich dynamical responses emerge from the interplay of drive, dissipation, and interactions.

Abstract

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 are 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. In this work, we investigate whether such a Dicke time crystal 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 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 time crystals, 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.

Figures

Equations

$${\hat{H}}_{int}=J\sum _{i=1}^N{\hat{s}}_i^x{\hat{s}}_{i+1}^x,$$

$${\hat{H}}_{ac}=\omega {\hat{a}}^{†}\hat{a}+{\omega }_0{\hat{S}}_z+\frac{2\lambda (t)}{\sqrt{N}}(\hat{a}+{\hat{a}}^{†}){\hat{S}}_x,$$

$$\dot{\hat{\rho }}=-i[\hat{H},\hat{\rho }]+\frac{\kappa }{2}(2\hat{a}\hat{\rho }{\hat{a}}^{†}-{\hat{a}}^{†}\hat{a}\hat{\rho }-\hat{\rho }{\hat{a}}^{†}\hat{a}),$$

$$\dot{\theta }$$

$$+J{\delta }_{pp}\sin \theta \sin \varphi \cos \varphi -J{\delta }_{pq}\sin \theta \cos \theta {\cos }^2\varphi ,$$

$$\dot{\varphi }$$

$$+J{\delta }_{qq}\cos \theta {\cos }^2\varphi -4J{\delta }_{pq}\sin \varphi \cos \varphi ,$$

$$\dot{a}$$

$${\delta }_{\alpha \beta }$$

$$ϵ$$

$${\Delta }_k^{\alpha \alpha }$$

$${\Delta }_k^{pq}$$

$$\begin{array}{rl}{\dot{\Delta }}_k^{qq}& =-\frac{4\lambda (t)}{\sqrt{N}}(a+a^{\ast })\frac{\cos \varphi }{\sin \theta }{\Delta }_k^{pq}\\ & -2J({\cos }^2\varphi -\cos k{\sin }^2\varphi ){\Delta }_k^{pq}\\ & -2J\cos k\cos \theta \sin \varphi \cos \varphi {\Delta }_k^{qq},\\ {\dot{\Delta }}_k^{pp}& =\frac{4\lambda (t)}{\sqrt{N}}(a+a^{\ast })\frac{\cos \varphi }{\sin \theta }{\Delta }_k^{pq}\\ & +2J({\cos }^2\varphi -\cos k{\cos }^2\theta {\cos }^2\varphi ){\Delta }_k^{pq}\\ & +2J\cos k\cos \theta \sin \varphi \cos \varphi {\Delta }_k^{pp},\\ {\dot{\Delta }}_k^{pq}& =\frac{2\lambda (t)}{\sqrt{N}}(a+a^{\ast })\frac{\cos \varphi }{\sin \theta }({\Delta }_k^{qq}-{\Delta }_k^{pp})\\ & +J({\cos }^2\varphi -\cos k{\cos }^2\theta {\cos }^2\varphi ){\Delta }_k^{qq}\\ & -J({\cos }^2\varphi -\cos k{\sin }^2\varphi ){\Delta }_k^{pp}.\end{array}$$

$$\begin{array}{r}{H}^{\prime }=H_{(LMG)}+{\hat{H}}_J,\text{with}{H}_{(LMG)}=-\frac{\lambda }{N}{\hat{S}}_x^2-g{\hat{S}}_z.\end{array}$$

$$\hat{U}={\hat{U}}_{kick}\exp \left[-i{\hat{H}}_{(LMG)}T\right],$$

$$\hat{H}=\omega {\hat{a}}^{†}\hat{a}+\frac{{\omega }_0}{2}\sum _i{\hat{\sigma }}_i^z+\frac{\lambda }{N}({\hat{a}}^{†}+\hat{a})\sum _i[e^{j\mathbf{k}\cdot {\mathbf{x}}_i}{\hat{\sigma }}_i^{+}+h.c.],$$

$$\hat{H}$$

$$n_{cav}$$

$${\sigma }_A^x$$

$$n_{cav}$$

$${\sigma }_A^x$$

Full text

I Introduction

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 Shapere and Wilczek (2012); Wilczek ; Bruno (2013); Nozieres ; Sacha ; Watanabe and Oshikawa (2015); Else et al. ; Khemani et al. (2016); von Keyserlingk et al. ; von Keyserlingk and Sondhi ; Khemani et al. (2017); Yao et al. (2017, ); Gong et al. (2018); Iemini et al. (2018); Tucker et al. (2018); Gambetta et al. (2019), and they have been currently realised with trapped ions Zhang et al. and solid state systems Choi et al. ; Rovny et al. (2018a, b). In most previous studies, TCs are realised in closed interacting quantum many-body systems, which are prone to heating towards an infinite temperature state under the action of periodic drive Lazarides et al. (2014); D’Alessio and Rigol (2014), therefore, a slowdown of energy absorption is customarily entailed via a disorder induced many-body localized phase Ponte et al. (2015, ); Zhang et al. (2016); Bordia et al. , or by fast driving Abanin et al. (2015); Mori et al. (2016); Else et al. ; Ho et al. (2017); Else et al. (2017).

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 Else et al. (2017). 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 Gong et al. (2018), 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 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, mean-field effect to the steady-state behaviour of a dissipative many-body system. We note however, that unlike the traditional venue for discrete time crystals Khemani et al. (2016); von Keyserlingk et al. ; von Keyserlingk and Sondhi ; Khemani et al. (2017); Yao et al. (2017), 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 Dicke-TC.

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.

II The model

We consider an chain of $N$ two-level atoms with short-range interactions among each other

[TABLE]

where $\hat{s}^{x,y,z}=\hat{\sigma}^{x,y,z}/2$, and $\hat{\sigma}^{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 Carmichael et al. (1973); Ritsch et al. (2013)

[TABLE]

where $\hat{S}_{x,y,z}=\sum_{i}^{N}\hat{s}_{i}^{x,y,z}$. We allow the light-matter coupling to be varied in time, $\lambda(t)$. Dissipation occurs when photons leak out of the cavity, as encoded by the quantum master equation

[TABLE]

for the total density matrix of the system, $\hat{\rho}$, where $\hat{H}=\hat{H}_{ac}+\hat{H}_{int}$, and $\kappa$ characterises the rate of photon loss.

When $J=0$, the above reduces to the well-known open Dicke model Hepp and Lieb (1973); Wang and Hioe (1973); Narducci et al. (1973); Duncan (1974); Garraway (2011); Nagy et al. (2010, 2011); Keeling et al. (2010); Bhaseen et al. (2012); Dalla Torre et al. (2013); Lang and Piazza (2016); Dalla Torre et al. (2016); Gelhausen et al. (2017); Shchadilova et al. (2018). As the coupling is only between the single photon mode and the collective spin operator, $\hat{S}^{x}$, the Dicke model is exactly solvable in the thermodynamic limit $N\to\infty$: its dynamics can be described by the mean-field motion of the photonic amplitude, $\langle a\rangle$, coupled to three classical degrees of freedom, $\langle\hat{S}_{x,y,z}(t)\rangle$, evolving on the Bloch sphere. When $J\neq 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 $\vec{S}$, which corresponds to the $k=0$ Fourier mode $\vec{s}_{k}\equiv\sum^{N}_{j=1}e^{-ikj}\vec{s_{j}}$, all other $k\neq 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 many-body 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 Lerose et al. (2018a, b, 2019). Specifically, we first perform a time-dependent rotating frame transformation to align the time-dependent $\hat{z}$-axis with the collective Bloch vector $\vec{S}(t)$, and then in this co-moving ‘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, $\epsilon(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|\vec{S}(t)|/N=1-\epsilon(t)$, similarly in spirit to approaches which incorporate self-consistently the effect of quantum fluctuations in the dynamics of a “condensate” Citro et al. (2015). The spin-wave density, $\epsilon(t)$, representing the total population of all $k\neq 0$ spin-wave 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 spin-wave 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 $\epsilon(t)$ in dynamics (see Appendix A).

III Dynamical responses

We impose a periodic modulation on the light-matter coupling $\lambda(t)$: during a first “bright-time”, $nT\leq t<(n+1/2)T$, we set $\lambda(t)=\lambda=\Omega$, and during the “dark-time”, $(n+1/2)T\leq t<nT$, we switch off $\lambda(t)$. Here, $\Omega=2\pi/T=(\omega_{0}+\omega)/2$ denotes the driving frequency. This periodic driving counteracts with the energy loss through photon leaking, and as previously found in Ref. Gong et al. (2018) for the collective Dicke model ($\hat{H}_{int}$ absent), it entails a period-doubling response in spin observables. The presence of such TC-like behaviour can be understood from the $\mathbb{Z}_{2}$ symmetry of the Dicke model with constant $\lambda$, under the parity operator $P=e^{i\pi(\hat{a}^{\dagger}\hat{a}+\hat{S}_{z}+\frac{N}{2})}$. For $\lambda>\lambda_{c}=\frac{1}{2}\sqrt{{(\omega^{2}+{\kappa^{2}}/{4})\omega_{0}}/{\omega}}$, a quantum phase transition that breaks the $\mathbb{Z}_{2}$ symmetry occurs, and the system enters a superradiant phase, featuring two steady states, with spin projection $\langle S_{x}\rangle/N=\pm X$, $\langle S_{y}\rangle=0$ and non-vanishing photon amplitude $a=\mp\sqrt{N}\lambda X/(\omega-i\kappa/2)$ ($X=\frac{1}{2}\sqrt{1-\lambda_{c}^{4}/\lambda^{4}}$). When $\delta\equiv(\omega_{0}-\omega)/\Omega=0$, the free evolution during the ‘dark-time’ accumulates a phase of $\pi$ for both $S_{x}=\langle S_{x}\rangle$ and $a$, and the system switches from one steady state to another, i.e., $S_{x}\rightarrow-S_{x}$ and $a\rightarrow-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 $\hat{H}_{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 $\kappa$, using the time-dependent 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/\Omega\ll 1$ and $0<\kappa/\Omega\lesssim 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, $\vec{S}$, crumble. In this regime of strong interactions, the system is prone to heating under the action of the Floquet driving (see also  111In this work, we choose as upper threshold for spin wave density the value, $\epsilon\sim 0.2$, in order to delimit the heating region from the other dynamical responses). However, as shown in Fig. 2(b), for moderate strengths of $J/\Omega\lesssim 1$, such effect of heating can still be remediated with sufficient dissipation $\kappa/\Omega\gtrsim 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 $\kappa$ (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 $\sim\sqrt{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 $\delta\neq 0$ is switched on, the collective ($J=0$) Dicke-TC is robust in a range of small $\delta\neq 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 $\delta=0$ limit, persists upon inclusion of many-body interactions, $J\neq 0$ (as an example, $\delta=-0.12$ is used for Fig. 2). The persistence of many-body Dicke-TCs on time scales much longer than $t_{\kappa}\sim 1/\kappa$ (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, $\kappa$. 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 222In principle, we cannot exclude that collapse of the order parameter might occur at much later times, and establishing the stability to infinite times may require another approach, which is beyond the scope of this work. However, this effect, even if present, would be of no practical relevance on experimentally accessible time scales.). 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. Lazarides and Moessner (2017)). 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 Yao et al. (2017); Zhang et al. .

When the rate of dissipation is intermediate, $\kappa/\Omega\sim 1$, we recognise a region of irregular (IR) dynamics, where the trajectory of $\vec{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 $\propto\lambda$, the system does not develop a $\langle S_{x}\rangle$ component sufficiently strong in order to counteract the dephasing induced by the “transverse field” $\propto\omega_{0}\hat{S}^{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 $\kappa/\Omega\gg 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 333$\hat{S}_{z}\gtrsim-N/2$, since part of the collective spin is dissipated in the bath of spin waves. (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 $\kappa$, 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, $\Omega$, only plays the role of an overall energy scale.

IV Metastable dissipative time crystal

For intermediate values of both $J$ and $\kappa$ (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 long-lived phenomenon (see footnote 444We have checked with a finite-size analysis that these conclusions do not depend on $N$.). The lifetime, $\tau$, of this MTC gradually decreases with $J$, following the empirical law, $\tau\propto\exp{[-A(J/\Omega)^{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 $\hat{z}$-direction of the Bloch sphere (corresponding to spin-waves vacuum). The MTC found here appears as a genuine interplay of period driving, dissipation and interaction. We remark here that it doesn’t result from a high-frequency expansion Mori et al. (2016); Abanin et al. (2015), and thus is distinct from the “Prethermal” Floquet TCs found in previous studies Else et al. (2017); Choi et al. , since $\Omega$ is an overall energy scale in our system (see Fig. 1). A possible explanation for the phenomenon is suggested by the dynamics of $\epsilon(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.

V Ferromagnetic driven-dissipative time crystal

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 $\pm\hat{x}$-direction, giving rise to robustness against heating and the overdamping caused by dissipation. Indeed, for non-perturbative values of $|J|\sim O(1)$, we find that a ferromagnetic ($J<0$) dissipative Dicke-TC can be stabilised at intermediate dissipation rate without significantly heating up the system (see footnote 555In this case, we switch off the $H_{int}$ term during the “dark-time”, in order to avoid inhomogeneous dephasing resulting from inter-spin interactions. ). 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 $\epsilon(t)$ even for $|J|\sim 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 driven-dissipative 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\simeq 0$), and from conventional (many-body) discrete TCs where dissipation is not a constitutive ingredient ($\kappa=0$). In this perspective, such ferromagnetic TC represents a non-equilibrium 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.

VI Summary and outlook

The Dicke model is currently engineered in several experimental platforms Black et al. (2003); Baumann et al. (2010, 2011); Brennecke et al. (2013); Klinder et al. (2015a, b). 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 inter-particle interactions such as using Rydberg atoms Norcia et al. (2018); Vaidya et al. (2018); Goban et al. (2014); Epple et al. (2014). 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 Welte et al. (2018). For coupling $\propto\lambda(\hat{a}+\hat{a}^{\dagger})\sum_{j}(-1)^{j}\hat{\sigma}_{j}^{x}$, a similar Dicke-TC response exists in this case but with anti-ferromagnetic ordering. Hence a $\hat{H}_{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. Choi et al. ). In the Appendix F we show that a Dicke model with coupling of the form $\propto\lambda(\hat{a}+\hat{a}^{\dagger})\sum_{j}[(-1)^{j/2}\hat{\sigma}^{+}_{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 Choi et al. ; Zhang et al. .

VII Acknowledgements

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’).

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 Lerose et al. (2018a, 2019) for a comprehensive discussion on the method.

We use the shorthand $O\equiv\langle\hat{O}\rangle$ for expectation values of operators $\hat{O}$, and define the coordinates for the collective spin vector, $\vec{S}=\frac{N}{2}(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, in terms of the polar, $\theta$, and azimuthal angle, $\phi$, on the Bloch sphere, which allows writing compactly the equations of motion:

[TABLE]

The effect of $H_{int}$ enters the dynamics of the collective spin vector via the spin-wave correlations $\delta_{\alpha\beta}(\alpha,\beta=q,p)$ and the spin-wave density $\epsilon$, which are in turn dynamically coupled to the collective spin and photon field. Here,

[TABLE]

with

[TABLE]

where $q_{k}$ and $p_{k}$ are the canonically conjugated bosonic variables associated with spin-waves with wave-vector $k\neq 0$.

Following the procedure in Lerose et al. (2019), we derive the equations of motion for the spin-wave correlations

[TABLE]

These quantities intertwine with the equations of motion (4), and represent the feedback of the non-equilibrium Gaussian fluctuations of spin waves on the motion of the collective spin (the $k=0$ mode). The self-consistent solution of Eqs. (4) and (9) yields the dynamics of the model. We notice that these equations of motion are derived in the thermodynamic limit Lerose et al. (2019), 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. (4), 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 $\propto 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. 4(c).

The spin wave density $\epsilon$ 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. (4) and (9). 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 $\epsilon(t)=0$ at any time, while for $J\neq 0$, the length of the collective spin is shrinked via $2|\vec{S}(t)|/N=1-\epsilon(t)$. A pictorial description of the method is provided in Fig. 5.

Appendix B Persistence of Dicke-TC

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 time-dependent spin-wave approach at stroboscopic times, when increasing the interaction strength $J$. As discussed in the main text, when the system is prone to many-body 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, $\epsilon\sim 0.2$. From the spin dynamics, we can also calculate the variance $\chi$ of the Fourier spectrum of $S_{x}+iS_{y}$, defined as $\chi=\int_{\nu>0}d\nu p(\nu)(\nu-0.5\Omega)/\int_{\nu>0}d\nu p(\nu)+\int_{\nu<0}d\nu p(\nu)(\nu+0.5\Omega)/\int_{\nu<0}d\nu p(\nu)$, where $p(\nu)$ is the height of the Fourier spectrum at frequency $\nu$. TCs feature stable subharmonic responses, with the spectrum composed of two sharp peaks located at $\nu=\pm 0.5\Omega$, and thus a small value of $\chi$. The variance for increasing $J$ is plotted in the bottom panel of Fig. 6. With a small $\kappa/\Omega<1$, the spin-wave density quickly grows with $J$, which is reduced with a large $\kappa/\Omega>1$. The variance $\chi$ 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 $\theta=0.5\pi$ and $\phi=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.

Appendix C Mean-field analysis of dynamical responses

To understand the dynamical responses in our driven-dissipative 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 dynamical 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.

Appendix D Summary of dynamical responses from time-dependent spin-wave analysis

To account for quantum many-body effects, we solve Eqs. (4) to (9) in the time-dependent spin-wave approach to obtain the spin dynamics. We explore over a range of parameters $\kappa$ 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 $\epsilon(t)\geq 0.2$ as heating (H), where the spin-wave treatment breaks down. IR regime is characterised with a large variance $\chi$ of the Fourier spectrum together with $\epsilon(t)<0.2$. As discussed in Ref. Gong et al. (2018) 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 $\kappa/\Omega\sim 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 $\chi$ and $\epsilon(t)$ remaining below $0.2$. The boundary between the OD and the Dicke-TC regime resembles the one in a mean-field 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 $\kappa$, 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 $\epsilon(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}$,

[TABLE]

The dynamics entailed by $H^{\prime}$ is periodically perturbed by a collective rotation along the $\hat{x}$-axis; the evolution operator reads in a period

[TABLE]

with $\hat{U}_{kick}\equiv\exp\left[-i\phi S_{x}\right]$. This protocol has been shown in Ref. Russomanno et al. (2017) to display TC behaviour when $\phi=\pi$ (for instance, in the dynamics of the stroboscopic transverse magnetization). The TC in this case is robust to displacements around the $\phi=\pi$ point, i.e. for angles $\phi=\pi\pm\delta$ (with $\delta>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 (11) 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

[TABLE]

with the coupling strength between atoms and cavity controlled by tuning the phase factor $\propto{\bf k}\cdot{\bf x}_{i}$; such coupling can be realized by choosing proper laser dressing in a quantum optics setup. When $e^{j{\bf k}\cdot{\bf x}_{i}}=1$, the hamiltonian (12) becomes the conventional Dicke model. In the case of $e^{j{\bf k}\cdot{\bf x}_{i}}\neq 1$, the coupling strength varies from site to site. We choose $e^{j{\bf k}\cdot{\bf x}_{i}}=(-1)^{i/2}$, therefore, the hamiltonian (12) becomes

[TABLE]

We divide the system into four sublattices, A, B, C and D, with ${\rm mod}[i,4]=\{0,1,2,3\}$, respectively. Within each sublattice, atoms are collectively coupled to cavity photons. In the limit of $N\rightarrow\infty$, the stationary solution of the corresponding master equation hosts multiple possible steady states, depending on the value $Q=\frac{\lambda^{2}\omega/\omega_{0}}{\kappa^{2}/4+\omega^{2}}$. The case $Q<1/4$ corresponds to the normal state without superradiance. When $1/4\leq Q\leq 1/2$, the system is superradiant, with photon occupation

[TABLE]

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

[TABLE]

with $X=\sqrt{Q^{2}-1/16}/Q$, and both with $\sigma_{A}^{z}=\sigma_{B}^{z}=\sigma_{C}^{z}=\sigma_{D}^{z}=-\sqrt{1-X^{2}}$. 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 (15), while the photon occupation in the second group is given by

[TABLE]

This group includes 8 different spin configurations

[TABLE]

with $X^{\prime}=\sqrt{Q^{2}-1/4}/Q$, and $Z^{\prime}=\sqrt{1-X^{\prime 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\leq Q\leq 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. (17). A possible procedure is as follows: We initialise spins close to $\sigma_{A}^{x}=\sigma_{B}^{y}=-\sigma_{C}^{x}=\sigma_{D}^{y}=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 $\tilde{\theta}$ at the end of each period $T$. $\tilde{\theta}$ consists of single-site rotations. Specifically, we apply $\hat{R}_{z}(\pi)$ rotation to all odd sites (A and C); even sites (B and D) are rotated depending on the measurement outcome of $\sum_{i\in A}S_{i}^{x}$ prior to the pulse: $\hat{R}_{z}(\pi)$ for negative outcome, and $\hat{R}_{y}(\pi)$ for positive outcome. The pulse strength is kept equal to cavity detuning.

If we measure $\sum_{i}S_{i}^{y}$ 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 $\tilde{\theta}$, we still see oscillations at a stable $4T$ period. An alternative order parameter is the inversion $\sum_{i}S_{i}^{z}$ 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. Gong et al. (2018) to persist for times longer than the decay time, therefore we can expect similar conclusions to hold for the quadruple period dynamics discussed here.