Entangled time-crystal phase in an open quantum light-matter system

Time-crystals are nonequilibrium many-body phases in which the state of the system dynamically approaches a limit cycle. While these phases are recently in the focus of intensive research, it is still far from clear whether they can host quantum correlations. In fact, mostly classical correlations have been observed so far and time-crystals appear to be effectively classical high-entropy phases. Here, we consider the nonequilibrium behavior of an open quantum light-matter system, realizable in current experiments, which maps onto a paradigmatic time-crystal model after an adiabatic elimination of the light field. The system displays a bistable regime, with coexistent time-crystal and stationary phases, terminating at a tricritical point from which a second-order phase transition line departs. While light and matter are uncorrelated in the stationary phase, the time-crystal phase features bipartite correlations, both of quantum and classical nature. Our work unveils that time-crystal phases in collective open quantum systems can sustain quantum correlations, including entanglement, and are thus more than effectively classical many-body phases.

In this paper, we consider a paradigmatic model describing atoms coupled to a light field via an excitationexchange interaction [60][61][62][63], see sketch in Fig. 1(a).This system can be realized in current cavity-atom ex- (a) An ensemble of two-level atoms, with ground state |g⟩ and excited state |e⟩, is driven by a laser with Rabi-frequency Ω and interacts via exchange of excitations with the light field of a cavity (coupling constant λ).The cavity is subject to photon loss at rate κ.(b) Phase diagram in terms of the time-averaged magnetization mz as a function of Ω and λ.It features a bistable regime terminating at a tricritical point, (Ω/κ, λ/κ) ≈ (0.17, 0.41).Here, the transition becomes of second order (see inset).Below the dashed line, the system does not possess a well-defined stationary state.(c) In the stationary phase the magnetization mz(t) approaches a constant value.(d) In the time-crystal phase, mz(t) undergoes persistent oscillations and the atoms and the light field are collectively entangled.
periments [64][65][66][67][68][69][70][71][72] and realizes a boundary time-crystal [62,[73][74][75] when adiabatically eliminating the light field.It features both a stationary superradiant and a timecrystal phase [63].Within the parameter regime in which one may expect the adiabatic elimination to hold, the nonequilibrium transition between the two phases is a second-order one.However, in contrast to the boundary time-crystal, the system features a bistable regime, characterized by the coexistence of a limit cycle and a stationary phase.This signals that the phase transition eventually becomes of first order [cf.Fig. 1(b)].Explicitly taking into account the light field further allows us -to best of our knowledge for the first time -to observe the emergence of quantum correlations, including entanglement, in a time-crystal regime.The existence of these correlations may motivate the development of alternative strategies for exploiting these phases for enhanced metrological applications [57,58].
For concreteness, we focus on a realization of the model in a cavity setup, as depicted in Fig. 1(a).The spins describe N two-level atoms with ground state |g⟩, excited state |e⟩, and energy splitting ω at .The bosonic operators a, a † are associated with the light field inside the cavity (frequency ω cav = ω at ).For later convenience, we define the quadrature operators q = i(a − a † )/ √ 2 and p = (a + a † )/ √ 2, such that [q, p] = i.The atoms are resonantly driven by a laser and, in the rotating frame, the system Hamiltonian is given by with Ω being the laser Rabi-frequency and λ the coupling constant providing the "rate" of the coherent exchange of excitations.The collective atom operators S ± are defined as − , with σ − = |g⟩⟨e|, S + = S † − .The factor 1/ √ N in front of the coupling term ensures a welldefined thermodynamic limit [13,76].Photon losses, at rate κ, are described by the dissipator [77,78] The full quantum state of the system, ρ, thus evolves according to the quantum master equation ρt = −i[H, ρ t ]+ L[ρ t ] and allows for the calculation of the expectation value of any operator O, as ⟨O⟩ t := Tr(ρ t O).

III. TIME-CRYSTAL PHASE TRANSITION
To demonstrate the emergence of a phase characterized by non-stationary asymptotic dynamics, we analyze the long-time behavior of our system, in the thermodynamic limit (N → ∞).To this end, we introduce the average "magnetization" operators m r /N for the atoms, with σ r being the Pauli matrices constructed from the states |g⟩ and |e⟩.For the light field, we consider the rescaled quadratures m N q = q/ √ N and m N p = p/ √ N .In the thermodynamic limit, both the atom and the lightfield operators m r = lim N →∞ m N r describe average properties of the system [53,76] and provide suitable order parameters.

A. Mean-field equations and fixed points
Since we are interested in the long-time regime, we derive the evolution equations for the average operators.We focus on physically-relevant initial states of the system [112], i.e., states with sufficiently short-range correlations.For such initial states, following the derivation put forward in Ref. [76], it is possible to show that, in the thermodynamic limit, the order-parameter dynamics is exactly captured by a set of nonlinear differential equations [76,113].These equations are the so-called mean-field equations and for the model considered they are given by The latter show that m 2 x + m 2 y + m 2 z is a constant of motion, which we set to one, and that assuming an initial state for which m x (0) = m q (0) = 0, results in having m x (t) = m q (t) = 0, for all times t > 0. The remaining operators evolve via the equations ṁy We note that, adiabatically eliminating m p (t), by setting the last of the equations above to zero and substituting the result in the other two equations, leads to the equations of motion for the boundary time-crystal model [19,53].A similar mapping holds also at an operatorial level [62].By setting the derivatives in the above equations to zero and using the constant of motion (m 2 x + m 2 y + m 2 z = 1), we find the stationary solutions to the mean-field equations, given by The stability of the stationary solutions can be analyzed by looking at the Jacobian matrix J, obtained by linearizing the mean-field equations around the stationary values.This matrix can be obtained by writing m c (t) ≈ m * c + δm c , with δm c being small, and considering perturbations only up to first-order.The linearized Jacobian matrix takes the form A stationary solution is stable if the real part of all the eigenvalues of the matrix J is smaller or at most equal to zero.The eigenvalues µ i of the matrix J are given by which immediately shows that the stationary state with positive m * z is unstable.The only stable stationary mean-field solution is the one with negative m * z [see also Fig. 1(c)].Such a stationary solution is physical only when Ω ≤ λ 2 /κ.Here, the light field becomes macroscopically occupied, ⟨a † a⟩ ∝ N (m * p ) 2 , denoting the superradiant character of the phase [13].For Ω > λ 2 /κ, no stationary solution exists (within the sector identified by the choice of the conserved quantities) and the system belongs to a time-crystal phase, as shown in Fig. 1(b-d).

B. Proof of existence of the limit cycle
The non-stationary behavior of the system in the timecrystal phase is, as we will show analytically by closely following the derivation in Section 8.5 of Ref. [114], the result of an emergent limit-cycle dynamics.To show the existence of limit cycles for the mean-field equations [cf.Eq. ( 3)], we first bring them into a more convenient form.We make use of the fact that m 2 y +m 2 z = 1 is a conserved quantity and thus Eq. (3) describes an evolution taking place on the surface of a cylinder.The dynamics of the system is then captured by which can be obtained by exploiting the ansatz obeying ṁy(z The above equation is closely related to the differential equations for the dynamics of the phase difference across a Josephson junction (see Section 8.5 of Ref. [114]).
With the restriction |m y | ≤ 1, we find again that the stationary solutions of Eq. ( 6) only exist for Ωκ < λ 2 .Above the critical value Ω = λ 2 /κ, we find persistent oscillations of the system witnessing a stable limit cycle to which all trajectories are attracted.To analyze the longtime behavior of the system in this parameter regime, we consider the nullcline Y = 2λ 2 κ m y − 2Ω, with |m y | ≤ 1, which defines a regime with vanishing derivative Ẏ = 0.For smaller (larger) values of Y , the derivative Ẏ is positive (negative) so that for long times all trajectories end up in a regime restricted to the strip [114].Given the periodicity of m y [cf.Eq. ( 5)], it is sufficient to consider values 0 ≤ θ ≤ 2π.For Ω > λ 2 /κ, we can fix y 2 < 0, such that the derivative of the angle θ = Y < 0 does not change its sign inside the strip.Thus, in the long-time limit a periodic solution can only exist within this strip.A limit cycle is a trajectory that starts at a point Y * and ends after one period at the same point P (Y * ) = Y * , where P is called Poincaré map [114].In order to show the existence of such a point inside the strip, we make use of the fact that P (y 1 ) > y 1 , ∀ y 1 < − 2λ 2 κ − 2Ω, which is due to the fact that the derivative Ẏ is strictly positive for y 1 < − 2λ 2 κ − 2Ω and thus Y cannot go back to the value y 1 [114].Similarly, we have Since the Poincaré map P is continuous and monotonic, there must thus exist a value Y * such that P (Y * ) = Y * , implying the existence of the limit cycle [114].It is also possible to show that the closed orbit is unique (for details we refer to Section 8.5 of Ref. [114]).
In Appendix A we further demonstrate that the emergent limit-cycle dynamics is associated with the spontaneous breaking of continuous time-translation symmetry.This shows that the considered system features a proper time-crystal phase [36].

C. Phase diagram and bistability
Having established the existence of a non-stationary regime, we now analyze in detail the nonequilibrium phase diagram of the system.We observe that, also within the parameter regime in which the stable stationary state of Eq. ( 4) is well-defined, the system can approach a limit cycle.This implies the existence of a region where the stationary phase [cf.Eq. ( 4)] and the time-crystal one coexist.Such bistable regimes usually occur for stationary phases (see, e.g., Refs.[120,121]) and are characterized through a stability analysis.However, in our case one of the two asymptotic solutions is a limit cycle.As such, to fully explore the bistable region we take an approach which exploits the coexistence between the two phases.To treat the latter on an equal foot- ing, we will focus on the time-averaged order-parameter mz = 1 t t 0 du m z (u), which converges to the stable value in Eq. ( 4) within the stationary regime while it gives an average over the oscillations in the time-crystal phase.
When Ω > λ 2 /κ, the system can only be found in the time-crystal phase.The curve Ω = λ 2 /κ thus provides one of the boundaries of the bistability region.To find the other boundary, i.e., the line separating the bistable regime from the stationary phase [cf.Fig. 1(b)], we probe coexistence behavior.The idea is as follows.We start at a point (Ω, λ) in parameter space, with Ω ≫ λ where only the time-crystal phase is stable [cf.Fig. 1(b)].We initialize the system in the state |ψ⟩, with all atoms in the excited state |e⟩ and the light field in the vacuum, and let it relax towards the asymptotic limit cycle.We then increase λ, in small discrete steps, in an adiabatically slow fashion, i.e., always giving the system sufficient time to accommodate into the new limit cycle.In this way, we can enter the bistable regime lying within the basin of attraction of the time-crystal phase.For sufficiently large λ, only the stationary phase is eventually stable.As shown in Fig. 1(b), this makes the second (upper) spinodal line emerge as the line where mz jumps from positive values, attained in the time-crystal phase, to the negative ones given by Eq. (4).A similar sweep through the phase diagram can be done by fixing λ.This procedure also shows coexistence of the two phases as apparent from Fig. 2(a).The two spinodal lines meet at a tricritical point, highlighted in Fig. 1(b).Beyond this point, the phase transition does not switch to a crossover as it usually happens, but it rather changes nature and becomes a second-order one, see Fig. 2(b).Note that the curve in Fig. 2 point with an infinite derivative [cf.Eq. ( 4)] and ii) as we calculate below and anticipate in Fig. 2(c), approaching the critical point from the stationary regime the system features a diverging susceptibility.

D. Characterization of the phase transitions in terms of bifurcations
The different phase-transition behavior can be related to the different types of bifurcations [114] occurring at the transition lines, see also animations provided as Supplemental Material [122].
In the regime where the adiabatic elimination is valid [lower left corner of the phase diagram in Fig. 1(b)] the system undergoes a phase transition at the critical line Ω = λ 2 /κ.Crossing the latter from the time-crystal phase the periodic solution is disrupted by the emergence of a pair of fixed points, a saddle and a node [122].Here, an infinite-period bifurcation [see also Fig. 3] occurs and the behavior of the system is analogous to that of the boundary time-crystal [115].
Above the tricritical point, the system can be found in a bistable regime, in which the stable stationary solution and the stable limit cycle coexist.When starting from the time-crystal phase and moving adiabatically slow inside and within the bistable regime, it is possible to remain within the basin of attraction of the time-crystal phase.In this case, when approaching the upper spinodal line [the line separating the bistable regime from the stationary one in Fig. 3] the limit cycle eventually hits the unstable (saddle) stationary solution.Here, a homoclinic bifurcation takes place and the system "jumps" into the stable stationary solution [122].On the other hand, coming from the stationary phase and increasing the parameters adiabatically slow, the system stays in the basin of attraction of the stable stationary solution, even within the coexistence regime.In this case, approaching the lower spinodal line [the line separating the bistable regime and the time-crystal phase in Fig. 3], stable and unstable stationary solutions coalesce (saddle-node bifur-cation).Beyond this line, the only attractor is the limit cycle.
The presence of different types of bifurcations (infinite-period bifurcation below the tricritical point [115], saddle-node and homoclinic ones above) explains the appearance of different phase-transition behavior [cf.Fig. 2(a-b)].Approaching the critical line below the tricritical point, the limit cycle acquires an infinite period and spends an infinite amount of time close to where the stable solution emerges.In this way, when passing from the limit cycle to the stationary solution the timeaveraged magnetizations change continuously, i.e., the system undergoes a second-order phase transition.On the other hand, above the tricritical point, when passing from one phase to the other, the system experiences sudden jumps between two already existing solutions, which live in different regions of the "phase space".This fact gives rise to a first-order phase transition with the associated jump of the order parameters.

IV. DYNAMICS OF QUANTUM FLUCTUATIONS
Average operators converge, in the thermodynamic limit, to multiples of the identity [123] and thus cannot carry information about correlations.The natural next step is thus to consider suitable susceptibility parameters.In analogy with classical central limit theorems, for the atoms we introduce the quantum fluctuation operators [51,[124][125][126][127][128][129]] whose variance, χ rr = ⟨F 2 r ⟩, provides the fluctuations of the order parameter m N r , that is, its susceptibility.The operators in Eq. ( 7) retain a quantum character in the thermodynamic limit.To understand this, let us consider the state with all atoms in |e⟩.The commutator [F N x , F N y ] = im N z is proportional to an average operator and, thus, converges in the thermodynamic limit to the multiple of the identity im z , with m z = 1, due to our choice of the state.This commutation relation identifies the limiting fluctuation operators, q A = lim N →∞ F N x and p A = lim N →∞ F N y , as two (bosonic) quadrature operators such that [q A , p A ] = i.Together with these atom fluctuations, we consider the light-field fluctuation operators q L = q − ⟨q⟩ and p L = p − ⟨p⟩ [54].The emergent two-mode bosonic description formed by the fluctuation operators R = (q A , p A , q L , p L ) T can be used to analyze correlations between the atoms and the light field [54].
To this end, we introduce the covariance matrix Σ uv = ⟨{R u , R v }⟩/2 and investigate its time evolution.Due to the dynamics of average operators, the commutation relation between the fluctuation operators associated with the atoms generically depends on time [51].To "remove" this dependence, we move to the frame rotating with the time-evolving average operators.Here, we can derive the Lindblad generator for the dynamics of the two-mode bosonic system related to quantum fluctuations (see Appendix B).The time-dependent Lindblad generator is of the form with the Hamiltonian where In the generator above, the map L * L is the dual of the map which is analogous to the one in Eq. ( 2) but with jump operator a L = (p L − iq L )/ √ 2. Under this dynamics the covariance matrix evolves according to the differential equation with s being the symplectic matrix of a two mode bosonic system [51] and the matrices encoding the dissipative dynamics of the system.
The emergent two-mode Hamiltonian, which can also be written as is time-dependent as a consequence of the timedependence of the instantaneous magnetization m z (t) and encodes both an excitation-exchange and a two-mode squeezing process.To show this, we represent the fluctuation operators q A , p A , q L , p L in terms of bosonic creation and annihilation operators.Due to the definition of the original quadrature operators of the light field, we write For the atoms, we instead recall that q A is the limit of F N x and p A the one of F N y and, in order to associate the annihilation operator with S − , we write Substituting these definitions into the Hamiltonian in Eq. ( 9), we find which makes apparent that the Hamiltonian can be decomposed into an excitation-exchange term, proportional to 1 − m z (t), and a two-mode squeezing term, proportional to 1 + m z (t).These contributions provide the only coupling between the atoms and the light field and, as we show below, can generate quantum correlations between the two subsystems.In contrast to the boundary time-crystal model [53], the dynamics of fluctuations is not fully dissipative due to the collective Hamiltonian in Eq. ( 1).Since the emergent dynamical generator is quadratic, quantum fluctuations remain Gaussian [130].

V. QUANTUM CORRELATIONS AND ENTANGLEMENT
From the time evolution of the covariance matrix Σ, we can calculate classical correlation, quantum discord [116,117,[131][132][133], as well as bipartite (collective) entanglement between the atoms and the light field [118,119], in the thermodynamic limit.Within the stationary phase, the asymptotic covariance matrix can be computed exactly as with the stable m * z [cf.Eq. ( 4)].This expression shows that the light field (described by the operators q L , p L ) is in the vacuum state, while the collective atom operators q A , p A are in a squeezed state.Eq. (11) shows no correlations between the atoms and the light field in the stationary phase.Yet, the atoms display spin-squeezing, with a squeezing parameter, ξ = |m * z |, which diverges (to zero) on the spinodal line separating the bistable regime from the pure time-crystal phase.The divergence (to infinity) of Σ 22 ∝ |1/m * z | is related to the divergence of the susceptibility close to the second-order phase transition [cf.Fig. 2(c)].Since fluctuations are in the frame aligned with the direction of the stable state in Eq. ( 4), Σ 22 , in the stationary regime and close to the phase transition line, is essentially the susceptibility of the orderparameter m z .
We now turn to the time-crystal phase.Here, there is no significant spin-squeezing in the atom ensemble.Moreover, it can be shown that the determinant of the covariance matrix increases indefinitely with time, which indicates that the state of the system becomes more and more mixed.Nonetheless, in this regime the atoms and the light field are correlated.This can be seen, for instance, through the one-way classical correlation.This quantity is a measure of the maximal information about one of the two subsystems, let us say the atoms, that can be gained by performing measurements on the other subsystem, in our case the light field.This one-way classical correlation, denoted as J A←L , is shown in Fig. 4(a) and demonstrates the existence of correlations in the time-crystal phase.Even more interestingly, also correlations of genuine quantum nature can be observed in this regime, as measured by the (one-way) quantum dis-cord D A←L = I − J A←L , with I being the mutual information between the atoms and light field.The quantum discord quantifies the amount of correlations which are not of classical nature.In Fig. 4(b), we show that in the time-crystal phase the quantum discord is nonzero throughout.(We report results for J A→L , D A→L in Appendix C.) Remarkably, a fraction of these quantum correlations is related to bipartite entanglement between the atom ensemble and the light field, which can be quantified through the logarithmic negativity E shown in Fig. 4(c).Both classical and quantum correlations display coexistence behavior, as for instance shown in Fig. 4(d), due to the coexistence of the uncorrelated stationary phase and the correlated time-crystal.
To conclude we note that Fig. 4(a-c) clearly shows that increasing the coupling strength λ between the atoms and the light field does not always lead to increased correlations.Indeed, for fixed Ω and κ, a too large coupling strength λ brings the system into the stationary uncorrelated phase.

VI. DISCUSSION
The system we have investigated is related to the wellknown boundary time-crystal model [19] through an adiabatic elimination of the light field [62,63,[73][74][75].For what concerns the atoms, it shows features which are similar to those of the boundary time-crystal.That is, we observe spin-squeezing in the stationary regime, and absence of quantum correlations among the atoms in the oscillatory phase [38,52,53,[55][56][57][58].However, explicitly considering the light field allowed us to uncover the existence of genuine quantum correlations in the timecrystal regime, even though the latter is characterized by a mixed state, established between atoms and light field.From a fundamental perspective our results demonstrate that time-crystal phases can display quantum correlations and are thus certainly not classical.Given the Gaussian character of the quantum state of the atoms and the cavity mode, the correlations we have investigated here may be accessed experimentally via measurements of two-point correlation functions.Our findings are valid in the thermodynamic limit.For a finite system, they are accurate up to a timescale t * (diverging for N → ∞).Beyond this timescale, the oscillations in single realizations of the dynamics dephase [115].The average state thus consists of the sum over all possible dephased limit-cycles and becomes asymptotically time invariant.This phenomenology is related, in the thermodynamic limit, to mode-softening and phase diffusion in time-crystals [37,134,135].
Finally, we note that the time-crystal phase appears to be related to lasing -since there is an inversion of population signalled by a positive magnetization mz > 0 [13] -even though the model does not possess a U (1) symmetry.This is due to the fact that the "pumping" is not incoherent but rather is implemented through external laser driving.However, the oscillations established are not harmonic [see Fig. 1(c)] and, deep in the time-crystal phase, |Ω| ≫ |λ|, the time-averaged magnetization mz tends to zero, i.e., there is no inversion of population [13].
Sec. III B through the existence of the limit cycle we will now focus on the latter.With our ansatz in Eq. ( 5) and an initial angle α, where m y (0) = sin α and m z (0) = cos α, we find m y (t) = sin Θ(t) and m z (t) = cos Θ(t), with the phase angle Θ(t) = θ(t) + α.For randomly sampled initial conditions the system can assume all phases in the limit cycle [see Fig. 5].Similarly to the discussion in Ref. [36], this witnesses the breaking of continuous symmetry in the time domain.Additionally, this also shows that the system always approaches the time-crystal phase demonstrating its robustness against varying initial conditions.

Appendix B: Dynamics of quantum fluctuations
In this Appendix, we give details on the derivation of the evolution of the covariance matrix for fluctuation operators, as well as on the transformation to the frame which rotates solidly with the main direction of the atom average operators.We then explicitly derive the dynamical generator for quantum fluctuations in such rotating frame.

Time evolution of the covariance matrix of quantum fluctuations
The derivation of the time evolution of the covariance matrix for fluctuation operators follows closely the one presented in Ref. [54].We start by introducing the vector of fluctuation operators RN = (F N x , F N y , F N z , q L , p L ) T in the time-independent frame.The covariance matrix of these fluctuation operators can be written as Σ = lim N →∞ K N + (K N ) T /2, where we have defined Here, the expectation ⟨•⟩ denotes expectation with respect to the state at time t.
We now consider the time evolution of this correlation function K N uv .First, we note that as well as q L = − ⟨q⟩ and ṗ L = − ⟨p⟩, and that ⟨ RN u ⟩ = 0. Taking the time derivative of K N uv then leads to where L * is the dissipator in the Heisenberg picture, i.e., the map dual to L. To proceed, we compute the commutators in the above equations which can then be rewritten in terms of fluctuation operators exploiting again that ⟨ RN u ⟩ = 0.A similar calculation also applies to the dissipative part in the above equation.As in Ref. [54], this gives rise to products of fluctuation operators and average operators.Making use of the fact that, in the thermodynamic limit, lim The evolution for the case considered in the main text is obtained by setting m x = m q = 0.

Covariance matrix in the rotating frame
We now focus on the case in which the system is initialized in the state with all atoms in |e⟩ and the light field in the vacuum.This gives m x (t) = m q (t) = 0 and m 2 y (t) + m 2 z (t) = 1.This is the initial state considered for producing the plots in the main text.Our task is now to find the time evolution of the covariance matrix in the frame which rotates solidly with the direction identified by the average operators.To rotate the reference frame of the atom operators back to the initial one, we need to find the rotation matrix which maps the instantaneous vector of the average operators m = [0, m y (t), m z (t), m q (t), 0] T into the one m = [0, 0, 1, m q (t), 0] T .Exploiting the conservation law m 2 y (t)+m 2 z (t) = 1, this matrix can be found to be the matrix Under this transformation, the symplectic matrix becomes .
The time evolution of the covariance matrix in the rotating frame can be calculated by taking the derivative of Σ = U (t) Σ(t)U T (t), which gives where For the considered initial state, the covariance matrix is given by the diagonal matrix Σ(0) = 1/2 diag(1, 1, 0, 1, 1).Starting from this covariance matrix, it is possible to see that the third row and the third column of the covariance matrix are not coupled with the remainder of the matrix.We thus define Σ to be the covariance matrix of the fluctuation operator q A (which is the limiting operator of the fluctuation F N x in the rotating frame), p A (which is the limiting operator of the fluctuation F N y in the rotating frame) coupled to the fluctuations q L ,p L (see also main text).
For such a matrix, the time evolution is given by the equation where X, s, c are the 4 × 4 matrices obtained by removing the third row and third column in Q, S, C respectively.

Dynamical generator for the quantum fluctuation dynamics in the rotating frame
We now want to find the generator for the dynamics of the two-mode bosonic system described by the vector of bosonic operators R = (q A , p A , q L , p L ) T .As done in Ref. [53], to this end we consider a time-dependent Lindblad generator on bosonic operators of the form with an ansatz for the Hamiltonian given by The dissipative part of the generator is essentially equivalent to the one of the original system, except that it now features the "rescaled" fluctuation operators of the light [cf.Eq. ( 8) in the main text].Using the generator W * A−L (t) to calculate the time evolution of the covariance matrix one finds Here, we have that the 4 × 4 matrix b is obtained by removing the third row and the third column from the matrix B introduced above.Comparing the above equation with Eq. (B2) gives that the generator correctly captures the dynamics of the covariance matrix if the relation In this Section, we describe how to calculate the correlation measures that we analyze in our work and we further present additional results on these.For details on the derivation of these measures, we refer to Refs.[116][117][118][119].
Given the two-mode covariance matrix Σ(t), we now show how to compute measures for the classical correlation, for the quantum discord, and for the logarithmic negativity.To start, we identify the relevant 2 × 2 minors of Σ as the matrices α, β and γ, such that 2Σ(t) = α γ γ β .
Here, the matrix α contains the variances of the atom fluctuations, β those of the light-field fluctuations, while as well as the function For a two-mode Gaussian state an expression for the one-way classical correlation, quantifying the information on the first mode obtained by measurements performed on the second mode, is given by while the quantum discord is where E min is defined as The quantities ν − and ν + are the symplectic eigenvalues of the matrix 2Σ, with ν − < ν + .These are found as the positive eigenvalues of the matrix 2isΣ.To compute the correlations J A→L and D A→L , quantifying the information about the light field that can be obtained from a measurement on the atoms, one can exploit the same definitions as above but exchanging the roles of α and β in all of the above relations.
In order to quantify the amount of bipartite entangle-ment between the atoms and the light field, we compute the logarithmic negativity.This is defined as where ν− is the smallest symplectic eigenvalue of the partially transposed covariance Σ P T = ΛΣΛ, where Λ = diag(1, 1, 1, −1).The latter is computed as the smallest positive eigenvalues of the matrix 2isΣ P T .

FIG. 1 .
FIG. 1. System and nonequilibrium phase diagram.(a) An ensemble of two-level atoms, with ground state |g⟩ and excited state |e⟩, is driven by a laser with Rabi-frequency Ω and interacts via exchange of excitations with the light field of a cavity (coupling constant λ).The cavity is subject to photon loss at rate κ.(b) Phase diagram in terms of the time-averaged magnetization mz as a function of Ω and λ.It features a bistable regime terminating at a tricritical point, (Ω/κ, λ/κ) ≈ (0.17, 0.41).Here, the transition becomes of second order (see inset).Below the dashed line, the system does not possess a well-defined stationary state.(c) In the stationary phase the magnetization mz(t) approaches a constant value.(d) In the time-crystal phase, mz(t) undergoes persistent oscillations and the atoms and the light field are collectively entangled.

FIG. 2 .
FIG. 2. Coexistence and critical behavior.(a) Timeaveraged mz as a function of Ω, for λ/κ = 0.8 [upper dotted line in Fig. 1(b)].The solid (dashed) curve is obtained by starting from the time-crystal (stationary) phase and moving "adiabatically" Ω towards the stationary (time-crystal) one.(b) Time-averaged mz as a function of Ω, for λ/κ = 0.3 [lower dotted line in Fig. 1(b)].Here, the phase transition is of second order [see also inset of Fig. 1(b)].(c) Critical behavior of the spin-squeezing and of the susceptibility in the stationary regime, as a function of (Ωκ/λ 2 ) 2 .The coordinates x, y refers to the frame in which the z-axis is aligned with the direction of the vector identified by the stable stationary values in Eq. (4).
FIG. 3. Phase diagram and bifurcations.Sketch of the phase diagram of the model specifying the types of bifurcation occurring at the critical and at the spinodal lines.

FIG. 4 .
FIG. 4. Quantum and classical correlations.(a) Time-averaged classical correlation J A←L , as a function of λ and Ω.For each value of Ω, the data are obtained by initializing the system in state |ψ⟩, evolving it with the smallest value of λ, and then adiabatically increasing the interaction parameter λ in discrete steps up to the largest values.The evolution time for each value of λ is κt = 5000 and coincides with the averaging window for the correlation measure.(b) Same as in panel (a) for the quantum discord DA←L .The latter shows that the time-crystal phase features quantum correlations.(c) Same as in panels (a) and (b) but for the logarithmic negativity quantifying the amount of entanglement between the atoms and the light field.(d) Coexistence of different bipartite entanglement behavior, as measured by the time-averaged logarithmic negativity Ē, for λ/κ = 0.8 [cf.upper dotted line in Fig. 1(b)], and different values of Ω slowly varied in discrete steps both starting from the stationary phase and the time-crystal phase.

FIG. 5 .
FIG. 5. Continuous symmetry-breaking in the timedomain.Distribution of the phase angle Θ in the timecrystal phase for fixed parameters (Ω/κ, λ/κ) = (0.5, 0.5) and 200 random initial conditions encoded in the angle α.The time at which Θ is evaluated is fixed for all initial values.As shown in the plot, Θ can assume all values between 0 and 2π, witnessing a continuous-time symmetry breaking.

,
through which we can define C = (D + D T )/2 and B = (D − D T )/(2i), and satisfied.Exploiting that s 2 = −I, we can invert the relation to find the Hamiltonian reported in the main text.Appendix C: Quantum and classical correlations

FIG. 6 .
FIG. 6.Additional results on quantum and classical correlations.(a) Time-averaged classical correlation J A→L , as a function of λ and Ω.For each value of Ω, the data are obtained by initializing the system in state |ψ⟩, evolving it with the smallest value of λ, and then adiabatically increasing the interaction parameter λ in discrete steps up to the largest values.The evolution time for each value of λ is κt = 5000 and coincides with the averaging window for the correlation measure.(b) Same as in panel (a) for the quantum discord DA→L .The latter shows that the time-crystal phase features quantum correlations.