A dissipative time crystal with or without $\mathbb Z_2$ symmetry breaking

We study an emergent semiclassical time crystal composed of two interacting driven-dissipative bosonic modes. The systems has a discrete $\mathbb Z_2$ spatial symmetry which, depending on the strength of the drive, can be broken in the time-crystalline phase or it cannot. An exact semiclassical mean-field analysis, numerical simulations in the quantum regime, and the spectral analysis of the Liouvillian are combined to show the emergence of the time crystal and to prove its robustness against quantum fluctuations.


Introduction
The advances in preparing and manipulating quantum matter in the laboratory during the past decades has led to a growing interest in out-of-equilibrium quantum phases [1][2][3][4][5][6][7][8][9][10]. Simultaneously, a great degree of attention has been devoted to the search for the spontaneous breaking of time-translational invariance [11][12][13][14][15][16][17][18]. Both efforts has converged on the realization of a time-crystalline phase of matter, the so called discrete or Floquet time crystals [19,20]. There, a many-body quantum system self-organizes and responds with a period different from the one imposed by the time-periodic external drive, breaking the discrete time-translational symmetry. An important characteristic of Floquet time crystals is that strong disorder is needed to induce many-body localization preventing the system from absorbing energy from the drive and heating up towards a featureless thermalized state.
Another type of time crystal dissipates energy to the environment instead of relying on disorder [21][22][23][24][25][26][27][28][29][30]. Interestingly, a subgroup of these are not driven at all [22] or the driving is such that the time-dependence can be completely eliminated by moving to a rotating frame [24,26,29,30]. In theses cases, the symmetry breaking is assessed with respect to the time-independent dynamical generator, usually the Lindblad superoperator L of a Markovian master equation, ∂ tρ = L(ρ). Such dissipative quantum systems have one or more attractive steady states that respect the time-translational invariance of L. In the thermodynamic limit, however, the steady state might never be reached, signalling that the continuous time-translational symmetry is broken. If one finds a time-periodic response on a function f (τ ) = lim t,N →∞ Tr[Ôρ(t + τ )] -N being the system size-for a suitable system operatorÔ one can say that a time crystal has been formed. Crucially, the period is not constrained to integer multiples of the drive period and can vary continuously with the system's parameters.
Most of the dissipative time crystals with continuous time-translation symmetry studied so far rely on long-range interactions [22,24,26], which occur naturally in systems with dipolar interactions or can be engineered, for instance, by coupling matter to a common resonant mode of a lossy cavity. This means those systems can be well described by mean-field equations in the infinite-volume thermodynamic limit, and so be regarded as emergent semiclassical time crystals. A few exceptions [29,30] rely on a different notion of 'thermodynamic limit' in effectively zero dimensions [31], where the number of bosonic excitations diverge in a system with one or a few bosonic modes and quantum fluctuations become negligible. These therefore belong to the same kind of time crystals.
The majority of the time crystals which have been studied so far have an underlying symmetry in addition to the time-translational symmetry, and the two are broken together. In Floquet time-crytals, for example, it is normally a global parity symmetry (Z 2 ), and this leads to long-range correlations in both time and space. For this reason the time-crystalline order is sometimes dubbed spatio-temporal order [16,17,[32][33][34].
In this work we study an emergent semiclassical time crystal in a dissipative system with two interacting bosonic modes which are driven and posses a minimal spatial Z 2 symmetry. We show that, depending on the strength of the drive, the time-crystalline phase is either accompanied by the Z 2 -symmetry breaking or it is not. This only occurs in a well-defined 'thermodynamic limit' in which the numbers of bosonic excitations diverge. By analysing how our quantum system scales towards this limit, we show the emergence of these time-ordered phases (i) proving that they are robust against quantum fluctuations as well as (ii) providing insight on the feasibility of observing long-lived oscillations in an experimental realization in which the system might be far from the 'thermodynamic limit'. A crucial feature in our approach is that we consider a system with nonlocal dissipation (also called dissipative coupling or collective dissipation), in which the two bosonic modes dissipate together. Nonlocal dissipation is also encountered in other models of dissipative time crystals [22,26,35]. An intuitive explanation might be that such collective process enhances synchronization between different parts of the system, which is in a sense what happens in continuous time crystals. In previous work [29] we have shown that collective dissipation causes a collective mode to be underdamped; this mechanism is the origin of the time crystalline behaviour. Even though it is often neglected in theoretical modelling, nonlocal dissipation occurs naturally in many systems in which there is one environment weakly-interacting with the whole system, being a crucial requirement for obeying quantum detailed balance [36]. Here we highlight the importance of nonlocal dissipation in the context of time crystals.

The model
We consider an open two-mode bosonic system with nonlocal dissipation which evolves according to the Lindblad equation ( = 1) where In a frame rotating with the pump frequency ω p the Hamiltonian readŝ Here,â i is the annihilation operator of the i-mode, ∆ = ω p −ω c is the frequency detuning between the pump frequency and the resonant frequencies ω c of the two modes, U is the interaction strength, J is the coupling between the two modes, F is the pump amplitude, and γ is the decay rate.
To better appreciate the effects the drive and the dissipation have over the two modes, we can rewrite the Lindblad equation in terms of bonding and antibonding modesâ B = ( In this new basis we can note that only the bonding mode dissipates and only the antibonding mode is driven, while the interaction (U ) couples both modes. This bosonic dimer (although with local dissipation) has been engineered using exciton-polaritons in microcavity pillars [37]. We illustrate our system in the context of micropillars in figure 1(a) showing how the bonding and antibonding modes look in the two coupled pillars. The bonding (red colored) and antibonding (blue colored) modes resemble the symmetric and antisymmetric wavefunctions, respectively, of the hydrogen molecule. Semiconductor microcavities are not the only platform available to study dissipative and interacting bosonic modes; circuit QED or optomechanical devices are also suitable. The drive we consider in this work can be achieved in any of these three platforms by using two coherent drives, one for each one of the 1 and 2 bosonic modes, with a π-phase difference, as illustrated in figure 1(a). Our nonlocal dissipation could be engineered in microcavity pillars by deliberately introducing defects in the overlap region of the two pillars in order to increase nonradiative losses in the bonding mode. Alternatively, it could be engineered in circuit QED devices by coupling two resonant modes to a microwave resonator at the same position, such that the amplitude and phase of the linear coupling between the resonator and each mode are the same. Moreover, reservoirengineered nonlocal dissipation has already been used in optomechanical circuits to achieve non-reciprocity [38,39]. There is an exact semiclassical limit for the system governed by equation (1) in which the occupation numbers of modes 1 and 2 diverge and the system's criticality becomes manifest [40,41]. It is defined by taking F → +∞ while keeping F √ U fixed. Is thus convenient to introduce a scaling parameter N to define F =F √ N and U =Ũ /N . From (1) we obtain the semiclassical mean-field equations where we have definedα A,B = α A,B / √ N . We can see that these equations are scale invariant and, in particular, they are exact ( â 1,2 = α 1,2 ) in the weak-coupling N → +∞ limit, which we study in this work. Note that in the laboratory, the interaction strength U is not easily controlled, yet is usually weaker than all other energy scales. Therefore, a large N limit is achieved solely by increasing the pump amplitude up to the level where the modes population become large enough such that the effective interaction energy U â † iâ i is relevant. For this reason F is the only parameter we vary throughout this work. In a previous work [29] we considered the same model but with a pump acting over the bonding instead of the antibondig mode. There we found that the timetranslational symmetry of the steady state was broken and it was accompanied by a first-order dissipative phase transition in the form of bistability. This means there was a region of parameter space where one could see long-lived oscillations with one of two different frequencies, depending on the initial conditions. In the current work the symmetries are different, so bistability is now replaced by a second-order phase transition due to the breaking of a Z 2 symmetry, as we explain in the following.

Symmetries
Our system has a discrete (also dubbed weak, see [42]) Z 2 symmetryâ 1 ↔ −â 2 (orα 1 ↔ −α 2 ) described by the bonding parity operatorP = (−1)â † Bâ B = ∞ n 1 ,n 2 =0 (−1) n 1 +n 2 |n 1 , n 2 n 2 , n 1 | (written in the Fock basis of the 1,2 modes). Note that due to the coherent drive, there is no U (1) phase invariance. For any finite value of N , the dynamical equation (1) has a unique steady state (i.e., time-independent) which is symmetricρ ss =Pρ ssP † . However, both the Z 2 symmetry and the time-translational symmetry of the steady state can be broken in the N → +∞ limit where the number of bosons in the system diverge. Naturally, then, our order parameters to witness spatial and time symmetry breakings in the semiclassical limit should be lim t→∞αB (t) = lim t,N →∞ N −1/2 Tr[â Bρ (t)] and any time-dependent function respectively. If f (t) is periodic, then the system would be in a time-crystalline phase. In the next section we will show that we find periodic oscillations for all values of the pump ampltiudeF , but the Z 2 symmetry is broken only for a particular region.

Mean-field semiclassical dynamics and symmetry breakings
In this section we analyse the dynamical behaviour of our semiclassical model (4). We look for fixed points and their stability, as well as the formation of limit cycles. In the first part we present the numerical results we obtain by solving (4) and then we give an analytical explanation. But first, since the nomenclature for nonlinear dynamics varies from one physics community to another, we first give a brief summary. The fixed points are the stationary solutions, i.e., ∂ tα fp A,B = 0, and their local stability is deduced from the eigenvalues of the Jacobian matrix obtained linearising the equations around them. If we express (4), and their complex conjugates, in vector notation as ∂ t α = G(α) whereα = (α B ,α * B ,α A ,α * A ) T , the expansion α → α fp + δ for δ the fluctuations vector leads us to the linear equation ∂ t δ = M δ where M = (∂G/∂α)| fp is the Jacobian matrix. Depending on the real part of the eigenvalues, the fixed point can be (locally) attractive, stable but non-attractive, or repulsive, corresponding to having all eigenvalues negative, at least one equal to zero, or at least one positive, respectively. A limit cycle is a periodic orbit in phase space; when a trajectory enters a limit cycle it remains there forever.

Numerical results
We find five regions in parameter space which are shown in figure 1(b), highlighted in different colors and annotated using Roman numerals. We plot the rescaled order parameter |α B |/ √ N of the Z 2 symmetry as a function of the rescaled driving amplitude F/ √ N γ, showing that in regions II, III, and IV the Z 2 symmetry is broken. The timetranslational symmetry is broken in all regions. In the following we explain the dynamics in each one of them.
In regions I and V there is a single non-attractive fixed point (dotted black line in figure 1(b)) which preserves the symmetry. This fixed point is never reached; all the trajectories go into a family of limit cycles on the antibonding mode alone, as the bonding mode amplitude goes to zero. These can be seen in panel (c) of the same figure, where we show a phase space portrait forF = 0.5 in region I. For panel (b) [as well as for (c), (d), and(e)] we allow 100 random initial conditions to evolve until they have become stationary and then we sampleα A,B (t) with a rate (γ∆t) −1 = 10 −2 .
In regions II and IV there are three repulsive fixed points (dashed black line); one of them preserves the symmetry and two of them belong to a symmetry-related pair (α B ↔ −α B ) breaking the symmetry. Here there is also a family of symmetrypreserving limit cycles revolving around the fixed point with α B = 0. Additionally, there are symmetry-breaking attractive limit cycles: in the low end of region II, a pair of them emerge from the fixed point as the pump amplitude rises, converging eventually to the fixed points with symmetry-breaking when region III is reached (see below). The same happens in regions IV in reverse. The different limit cycles can be seen in figures 1(d.i) and (d.ii), where we plot a phase space portrait forF = 0.95. Note we now plot the phase space of both modes as the bonding mode does not vanish.
Finally, in region III, the two symmetry-breaking repulsive fixed points of II and IV become attractive. There are only symmetry-preserving limit cycles in this region and they revolve around the repulsive fixed point with α B = 0. The fixed point (indicated by an arrow) and the limit cycles can be seen in panel (e), where we plot a phase portrait forF = 1.2. We do not show the bonding mode in this case, as it goes to zero (for the limit cycles) or to the finite fixed point value indicated by the order parameter in panel (b).
We can connect the image in figure 1(a) with the dynamics just described to have a pictorial understanding. Succinctly we can say: In regions I and V the red (symmetric) mode is not occupied and the blue (antisymmetric) mode always oscillates; in regions II and IV both red and blue modes always oscillate; and in region III either the blue mode oscillates and the red mode is empty or both modes are occupied without any oscillations.
In figure 2 we show examples of limit cycles for regions I and II by plotting the modes' amplitudes as a function of time. For the moment we just focus on the curves for N = ∞, which correspond to the semiclassical case. In panel (a), for low pump (F /γ = 0.5), we can see that the limit cycle is reached when the bonding mode goes to zero on a short time scale. In panels (b) and (c) we show the same quantities for a pump amplitude F/ √ N γ = 0.95, inside region II. In panel (b) the initial condition is such that a symmetry-preserving limit cycle is reached, while in (c), it is such that a symmetry-breaking limit cycle is obtained. On the right hand side of the figure, we show the Fourier transforms of the antibonding mode α A / √ N and we observe there are equidistant frequency peaks akin to frequency combs. This implies the corresponding periods are commensurate with each other, meaning a single common period governs the oscillations. Nevertheless, the elementary separation between frequency peaks has a nontrivial dependence on all parameters.

Analysis
Some intuition about the dynamical behaviour can be gained by noting that, if J−∆ > 0, which is our case, a fixed point of (4) is given bỹ whereα fp A is real and negative and the solution of the depressed cubic equatioñ U (α fp A ) 3 + (J − ∆)α fp A + √ 2F = 0. This is the symmetry-preserving fixed point found for allF in the numerical analysis and shown in figure 1(b). The eigenvalues of the Jacobian matrix which determines its linear stability are given by The eigenvalues λ ± A are purely imaginary meaning (5) cannot be attractive. Instead, a family of limit cycles can be reached: wheneverα B = 0 the two modes decouple, and the dynamical equation ofα A (t) corresponds to the semiclassical equation of a coherently driven, dissipationless nonlinear harmonic oscillator The limit cycle attained depends on all parameters and the value ofα A (t) at the time whenα B vanishes. A similar effective decoupling mechanism has been studied in [29].
There is an instability signalled by Re[λ B + ] > 0 where two symmetry-breaking fixed points emerge. This inequality is satisfied in a region where the following two conditions are fulfilled: For the parameters we use in figure 1 (same throughout the manuscript) these two inequalities are simultaneously satisfied in the regionF /γ ∈ [0.927, 1.596], which corresponds to regions II, III and IV. Interestingly, the bonding mode amplitude need not be zero for the effective decoupling of the two modes. When the amplitude is small (|α B | 1), equation 4(b) decouples from 4(a), resulting in equation (8) which drives the antibonding mode into periodic oscillations. At the same time, equation (4a) becomes linear in the bonding mode while the antibonding mode acts as a drive giving Note that only the time dependence ofα A is written explicitly. This is to highlight that, since α A (t) is periodic, it will cause the bonding mode to oscillate with the same period -this is, it is acting as a Floquet driving. This explains why in regions II and IV a family of symmetry-breaking limit cycles organise around the unstable fixed points emerging from the instability outlined in (9). We remark on the importance of a pure nonlocal dissipation. If local dissipators in the form of κ 2 i=1 D[â i ](ρ) (with small κ γ) were to be added to (1), the oscillations in the bonding and antibonding modes would be damped for all valuesF . Even though the system would preserve the Z 2 symmetry, one would be able to observe long-lived oscillations only up to a time ∼ κ −1 .
Having given an analytical explanation for the phase space of our system in the N → +∞ limit, we remark that an analysis that does not go beyond mean-field is incomplete: quantum fluctuations may well destroy the semiclassical limit cycles. A transparent example can be found in [43], where a linearisation method is developed to show that quantum fluctuations smear out the semiclassical limit cycle of the Van der Pol oscillator. In order to prove the robustness against quantum fluctuations, we proceed in the next section by solving our system exactly deep in the quantum regime (N = 1) and carrying out a 'finite-size' expansion in terms of the scaling parameter N .

Quantum dynamics
For our numerical calculations in this section, we have appropriately truncated the Fock space in the 1-2 basis ensuring convergence in the results for each value of the scaling parameter N . For the time evolution of expectation values, we solve (1) using a photoncounting unravelling [44] of the master equation and averaging over 5 × 10 3 quantum jump trajectories, recovering to an excellent accuracy the full dynamics of the density matrix. We also use the Wigner phase space representation and the Truncated Wigner Approximation (TWA) [44] -we invite the reader to see Appendix A for details.
By studying the time evolution of the rescaled expectation values â A,B (t)/ √ N for increasing values of N , we gain insight on both symmetry breakings we are interested in. Firstly, the emergence of periodic oscillations in | â A |(t)/ √ N for all the pump amplitudes we consider would highlight the time-translation symmetry breaking. Secondly, if the Z 2 symmetry is to be broken for anF window (i.e., regions II-IV), the response of | â B |(t)/ √ N should provide evidence for critical slowing down in this region due to the necessary closing of the Liouvillian eigenvalue gap [45] (more to this below).
We compare the previously shown semiclassical limit cycles with the dynamics of the quantum regime in figure 2. It clearly depicts the emergence of periodic oscillations in the antibonding mode as N is increased. Recall that panel (a) corresponds to a low pump amplitude in region I, while panels (b) and (c) correspond to region II with different initial conditions. Looking at the oscillations in the bonding mode in (b) and (c), we can see that also in the quantum regime a symmetry-preserving or symmetrybreaking limit cycle is approached with increasing N . Furthermore, we can observe in the Fourier spectrum of figures 2 (a), (b), and (c) that the long-lived oscillations have frequency peaks matching those of the semiclassical limit. For N up to 50 we have used quantum jump trajectories while for the very high N = 200 we have used the TWA which is a very good approximation for weak interactions (U =Ũ /N = 0.005).
Thanks to the linearity of the Lindblad equation we can expand its formal solution in eigenmodes of the superoperator L aŝ ρ(t) = e Lt (ρ(0)) = j≥0 c j e λ j tρ j =ρ ss + j≥1 c j e λ j tρ j , whereρ j are the eigenmodes, λ j the complex eigenvalues, and c j coefficients that depend on the initial condition. As our notation suggests, λ 0 = 0 is associated with the steady state eigenmodeρ 0 =ρ ss . The other eigenvalues have negative real part, hence |Re λ j | correspond to decay rates of the transient modes, while Im λ j are frequencies. Since the steady state preserves the symmetry of L, this picture suggests that to have a phase transition at least one non-zero eigenvalue must vanish in the N → +∞ limit (or in a thermodynamic limit, in general). This is analogous to the gap closing in a second-order phase transition of a closed quantum system [46]. In order to obtain the eigenvalues of a general Liouvillian superoperator L, one normally proceeds by first writing it as a matrix,Ľ, and then diagonalizing it. In our case, thanks to the discrete Z 2 symmetry, theĽ obtained from (1) decomposes into two block matrices asĽ where (±) refers to the two eigenspaces corresponding to eigenvalues ±1 of the symmetry operatorP ≡P ⊗P * that commutes withĽ. Details can be found in Appendix B.
The usefulness of this is twofold. Firstly, it reduces the computational cost of finding the eigenvalues. Secondly, eigenmodes in the (+)-subspace are Z 2 -symmetric while eigenmodes in the (−)-subspace are antisymmetric. This means we should be able to find evidence of the time-translational symmetry breaking in the spectrum ofĽ + alone, and of the Z 2 symmetry breaking in the spectrum ofĽ − . In figure 3 we show the non-zero eigenvalues of L with smallest absolute real part as a function of the pump amplitude for different values of N . In panel (a), we show an eigenvalue (λ + 1 ) ofĽ + which is purely real and goes to zero as N increases. In panel (b), for the same symmetry subspace, we show a second eigenvalue (λ + 2 ) which is complex. When N is increased, its real part tends to zero while its imaginary part converges to a finite value. This eigenvalue is responsible for the emergent oscillations for all the pump values we consider. In order to appreciate this, we have included in panel (b) the frequency predicted from the semiclassical linearised equations, i.e., from (5) and (7). Clearly |Im λ + 2 | approaches the semiclassical frequencies. The smallest eigenvalue (λ − 1 ) ofĽ − tends to zero with increasing N in the region where the semiclassical analysis predicts the broken Z 2 phase. This is shown in 3(c), where we have delimited between vertical blue dashed lines regions II, III and IV (in the plot for the real part). The real part of the eigenvalue has an inverted parabolic shape approaching zero with increasing N , while the imaginary part is zero throughout this region. The asymptotic vanishing of this eigenvalue explains the critical slowing down in the convergence to the steady-steady expectation values of non-symmetric operators, as shown for the bonding mode in figure 2(c). Figure 2(b) does not show the same critical slowing down in the bonding mode in spite of sharing the same pump amplitude. This can be understood by recalling that the coefficients c j in (12) depend on the initial conditions; in the latter case the initial state does not couple to the eigenmode associated to λ − 1 . Figure 4 shows the scaling of Re λ − 1 and Re λ + 2 with N in a log-log scale. For Re λ − 1 , we take the largest value of each curve shown in figure 3(c), and we extract Re λ + 2 for the same values of the pump amplitude. We see that both eigenvalues are well fitted The two smaller eigenvalues ofĽ + . λ + 1 is real and λ + 2 is complex, with its imaginary part approaching the semiclassical frequency (7) (black dashed line) as N is increased. (c) The smaller eigenvalue ofĽ − . Its real part approaches zero by increasing N in a region where its imaginary part vanishes. The region between vertical dashed blue lines corresponds to regions II, III, and IV where the semiclassical analysis predicts breaking of the Z 2 symmetry. The markers are just a guide; each curve contains 76 points.  by an algebraic scaling λ ∝ N β with β = −0.52 for λ − 1 and β = −0.64 for λ + 2 , showing that their real part vanish in the N → +∞ limit.
In region III the semiclassical prediction is not completely accurate. It predicts there is either a steady state with Z 2 broken or a limit cycle only in the antibonding mode. Yet from figures 3 and 4 one can deduce limit cycles will always be found in both modes throughout region III. Even though one could fine tune the initial condition such that a symmetric or non-symmetric steady state is reached -in the N → ∞ limit when the gaps are closed-any perturbation or time-delayed two-point measurement on the system would be enough to make the bonding and antibonding occupations oscillate forever. To illustrate this point, we show in figure 5 the second-order coherence of the bosonic modesâ 1 (left panel) andâ B (right panel) forF = 1.2 in region III, and two values of N . In the long time limit, the coherence is given by where j can be 1, 2, A, or B, and we have definedρ j =â jρssâ † j / Tr[â † jâ jρss ]. The last equality in the above equation emphasises the following: the first detection transformŝ ρ ss intoρ j , which is then evolved up to a time τ where the second detection occurs. Noting thatρ j can be expanded as a linear combination of L's eigenmodes, as in (12), we can conclude that the long-lived eigenmodes are probed by the measurement [47].

Discussion and outlook
The time crystal discussed in this work pertains to the class of emerging semiclassical time crystals also discussed in [22][23][24][25][26][27][28][29][30]. The mechanism behind it can be traced to the effective dynamical decoupling between the bonding and antibonding modes. The key ingredients are that these modes are nonlinearly coupled and that only one of them is explicitly damped (c. f. (3)). The other (non-damped) mode can evolve autonomously when the population in the damped one is small. We note that this behaviour is closely related to the one found in frequency combs in the weak-lasing regime of dissipatively coupled condensates [48]. The mean-field model for the two condensates considered in that work, can also be obtained by the semiclassical limit of the Lindblad master equation which considers incoherent (W ) instead of coherent pump. The dissipative coupling is crucial, and the limit cycles are found when the population in the two modes are small and the nonlinearity α becomes inefficient, like in our case. This phenomenon can be generalised to spatially extended configurations of bosonic modes. In a ring, for instance, one would need that the different linear modes (with well defined angular momentum) have different dissipation rates. Then, coherently pumping one of the linear modes would result in periodic long-lived oscillations in the mode(s) with smallest decay rate(s). If one of them has no decay channel -like in the case discussed in this work-then a time crystal would be found.

Acknowledgments
We are thankful to the developers of the Quantum Toolbox in Python (QuTiP) [49,50], as we have used it for most of our calculations in the quantum regime.
Our system has the discrete Z 2 symmetry which can be expressed as L(PρP † ) = P L(ρ)P † for anyρ, or equivalently, (P ⊗P * )Ľ =Ľ(P ⊗P * ). This is, they commute. SinceP has eigenvalues ±1, we can split the Liouvillian into two blocks, whereĽ ++,++ corresponds to projecting onto the eigenspaces ofP andP * with eigenvalue +1 on the right and left hand sides ofĽ, and so on.