Time crystallinity in open quantum systems

Time crystals are genuinely non-equilibrium quantum phases of matter that break time-translational symmetry. While in non-equilibrium closed systems time crystals have been experimentally realized, it remains an open question whether or not such a phase survives when systems are coupled to an environment. Although dissipation caused by the coupling to a bath may stabilize time crystals in some regimes, the introduction of incoherent noise may also destroy the time crystalline order. Therefore, the mechanisms that stabilize a time crystal in open and closed systems are not necessarily the same. Here, we propose a way to identify an open system time crystal based on a single object: the Floquet propagator. Armed with such a description we show time-crystalline behavior in an explicitly short-range interacting open system and demonstrate the crucial role of the nature of the decay processes.

Time crystals are genuinely non-equilibrium quantum phases of matter that break time-translational symmetry. While in non-equilibrium closed systems time crystals have been experimentally realized, it remains an open question whether or not such a phase survives when systems are coupled to an environment. Although dissipation caused by the coupling to a bath may stabilize time crystals in some regimes, the introduction of incoherent noise may also destroy the time crystalline order. Therefore, the mechanisms that stabilize a time crystal in open and closed systems are not necessarily the same. Here, we propose a way to identify an open system time crystal based on a single object: the Floquet propagator. Armed with such a description we show timecrystalline behavior in an explicitly short-range interacting open system and demonstrate the crucial role of the nature of the decay processes.

Introduction and outline
Statistical mechanics has been extremely successful in describing the behavior of systems at equilibrium and, occasionally, even the relaxation towards it. During the last decades, countless efforts have been devoted to genuinely non-equilibrium systems. In particular, a lot of attention has been drawn to non-equilibrium Floquet systems, i.e. systems undergoing coherent time-periodic dynamics. Those systems have found numerous applications that go from thermal machines and transport [1][2][3][4][5], to Floquet engineering [6][7][8], as well as the discovery of non-equilibrium phases of matter [9][10][11][12][13][14][15][16][17][18]. This work focuses on the latter, the so-called time crystals.
A system has discrete time-translational symmetry if the generator of the evolution, at any time t, is invariant under the transformation t → t + T , where T is the period. The pioneering ideas of time-crystals proposed in [19] and polished by subsequent discussions in [20], led to the concept of discrete time crystals (DTCs), first put forward in [11,15] and then experimentally realized in [13,14]. A discrete time crystal is a system that breaks discrete time-translational symmetry showing robust subharmonic response (to be precised below). Closed quantum systems might display subharmonic oscillations in a wide variety of scenarios; e.g., from Rabi oscillations of quantum optical systems to Bloch oscillations in lattices. Hence, criteria on how to identify a time crystal are essential to understand this non-equilibrium phase of nature. In closed systems, a discrete time crystal phase is characterized by an observable, O, acting as an order parameter whose expectation value f (t) = tr [Oρ(t)] must fullfill the three following conditions: [10,21,22]: (I) Time-translation symmetry breaking: the order parameter is less symmetric than the Hamiltonian, i.e. f (t + T ) = f (t) when H(t + T ) = H(t). For the DTC we recover f (t) after an integer number of periods N > 1, f (t + N T ) = f (t).
(II) Rigidity of the oscillations: f (t) shows a fixed oscillation period N T without fine-tuned Hamiltonian parameters. Equivalently, the oscillations should lock-off at frequency 2π/(N T ).
(III) Persistence to infinite time: the non-trivial oscillation with period T must persist for infinitely long time in the thermodynamic limit.
Nevertheless, it is not known whether this phase of matter survives the action of an environment or, even more interesting, if the environment can help to stabilize it [12]. For instance, the discrete time crystal appearing in a disordered one dimensional Ising spin chain (the so-called πSG) cannot survive the coupling to an environment [23]. As a result, efforts have been redirected towards the study of open mean-field models with long-range interactions [17,18,24] or short-ranged perturbations from them [25] coupled to an environment. In such models, some signatures of time crystallinity have been theoretically predicted, but there is a still a controversy on whether the mean-field description used on those models could sweep under the carpet part of the relevant physics which would destroy the time crystalline order. Part of the controversy is related to the fact that there is no well-posed description of what an open system time-crystal should be.
In this work we propose plausible criteria to define and characterize discrete time crystals in open systems governed by a Lindblad master equation. Our study is based on the analysis of the so-called Floquet propagator and the properties the associated Liouvillian must have in order to support time crystallinity. Further, we investigate the stability mechanisms and analyze how time crystals can be implemented in open systems. Also, there has been some discussion around the possibility that only long-range models can exhibit time crystalline order in open quantum systems. Our analysis also shows that long-range interactions (or mean-field models) are not crucial features to observe timecrystalline behavior. Nonetheless, our findings show that long-range jump operators are relevant in order to have subharmonic oscillations that are more robust to rotation errors.
The outline of the article is as follows: In Sec. 2, we present in detail the most relevant tools and concepts used throughout. In Sec. 3, we introduce the definition of time crystals in open systems and derive some important properties to further characterize them. In section Sec. 4, a collection of low dimensional open system examples is introduced in order to build up some intuition on how to reach stability of time crystals. Sec. 5 and Sec. 6 form the main body of this work, there we present a short range many-body model (XY model) and evaluate its properties and validity as an open system time crystal using the definitions proposed in Sec. 3. To this aim we derive its corresponding master equation and solve it numerically. Finally we present our conclusions.

Open quantum system dynamics: concepts and tools
In this section, we introduce the tools of -Markovian-open quantum systems and fix the notation used throughout the article. In what it follows, H denotes a Hilbert space of dimension d H , Op(H) denotes the set of operators acting on H and S(H) ⊂ Op(H) the set of states. Also, SOp(H) denotes the set of all superoperators i.e., all linear maps between elements of Op(H).

Open systems dynamics
Consider a closed time-independent physical system whose dynamics is generated by a Hamiltonian H ∈ Op(H). Then, the Schrödinger equation: and the state |Φ undergoes unitary dynamics. Often, one is interested only on the dynamics of a reduced set of the degrees of freedom commonly referred to as system (S). Complementary to those, there are bath (B) or environment degrees of freedom. A partition H = H S + H SB + H B is always possible, where the subscript indicates the degrees of freedom of the system (S), the bath (B) or both at the same time (SB) are present. Under fairly general conditions, the evolution of the system is generated by the so-called Lindblad master equation: where L denotes the Liouvillian superoperator and L α are the jump operators. Note that this is a hermiticity preserving equation, i.e. Lρ † = (Lρ) † . In particular, Eq. (2) can be derived under three key approximations [26,27]: the state of the system is initially uncorrelated from that of the environment, and the coupling between the system and the environment is weak (Born/weak-coupling approximation), and the environment equilibrates fast (Markov approximation). In addition, the fast rotating terms are usually disregarded (secular approximation). Such evolution generates then a family of completelypositive and trace-preserving (CPTP) maps satisfying E(t + s) = E(t)E(s) and E(0) = I, where I is the identity map. Remarkably, different Lindblad equations can give rise to the same CPTP map. For instance, Eq. (2) is invariant under the gauge transformation, for a ∈ C, Since linear operators form a vector space, a useful tool to characterize operators is the so-called vectorization procedure that maps any linear operator to a vector, i.e., O = The same transformation holds for linear superoperators, being now regular operators on H ⊗ H. In an abuse of notation, we denote with the same symbol both, O ∈ SOp(H) and O ∈ Op(H ⊗ H). The inner product in the extended Hilbert space is given by the Hilbert-Schmidt product defined by ⟪X|Y⟫ = tr X † Y . It automatically introduces the notion of adjoint superoperator fulfilling with the property O = X·Y ⇒ O ‡ = X † ·Y † . In particular, the adjoint Lindblad equation yields Using the vectorization procedure, L admits a matrix representation of dimension d 2 H × d 2 H . Hence, there is a set of d 2 H eigenvalues {λ µ } that can be found as the roots of the characteristic polynomial P L (λ) = det(L − λI). The ordinary left and right eigenvectors are defined as non-trivial solutions of the equations In the following, we outline some important mathematical properties of L that will be needed later. For simplicity, we assume that L is diagonalizable and refer the reader to App. A and references [26,[28][29][30] for an extended discussion.
(i) The eigenvalues of L are either real or come by conjugate pairs. Also, positivity requires the eigenvalues to have negative real part Reλ µ ≤ 0.
(ii) Ordinary eigenvectors of different eigenvalue are linearly independent.
(iii) The ordinary eigenvectors of L and L ‡ can be chosen bi-orthogonal, i.e. ⟪l µ |r ν ⟫ = δ µν . More compactly, for diagonalizable L we have W ‡ l W r = I.
(iv) For t ≥ 0, the solution of Eq. (2) for an initial state ρ(0) is given by where we allow for time-dependent generator L(t), and the symbol T stands for the time-ordering operator.
(v) For any L there is always one eigenvalue λ 0 = 0, with left eigenvector ⟪1|. The corresponding right eigenvector ρ ∞ = r 0 fulfills E(t)ρ ∞ = ρ ∞ and is often referred to as the steady-state. The steady-state might not be unique.
As stated in property (v), the steady-state may not be unique. Hence, it is useful to introduce the asymptotic subspace: In some contexts, the set of eigenvalues associated to the asymptotic subspace are also called the peripheral spectrum. Note that the elements of the asymptotic subspace are not always steady, but rather non-decaying. Also, they are general elements of Op(H) and not always proper quantum states. The asymptotic space can be always diagonalized, otherwise the dynamics will explode as t → ∞. From now on, we denote the elements of the asymptotic space by Ψ µ ∈ As(H) and we just refer to the orthogonal projection as the decay space (see Fig.1). We also define the dissipative gap ∆ as: that fixes the time-scale of convergence towards the steady state of the system.

Conserved quantities of the evolution
In closed systems, any symmetry of the Hamiltonian is a conserved quantity of the evolution. In dissipative systems, this is not always the case, and the relation between symmetries and conserved quantities is, in general, more complex. Here we formulate a simplification of the correspondence between steady-states and conserved quantities presented in [31]: Basically, given |Ψ µ ⟫ the eigenvectors of L spanning As(H) with eigenvalues λ µ = iω µ there is a corresponding set of conserved quantities ⟪j µ | such that, for any initial state |ρ(0)⟫ , asymptotically we have where ⟪j µ |ρ(0)⟫ is the imprint of the initial state on the conserved quantities. For instance, when there is only one steady-state, ⟪j 0 | = ⟪1| is the only conserved quantity and the trace tr [ρ(0)] is conserved throughout the evolution. An analytic expression for the conserved quantities that we will use throughout Sec. 4 was given in [31].

Periodic time-dependent Hamiltonians in open systems
When the Hamiltonian depends parametrically on time, i.e. H → H(t), the dynamics is governed by the unitary operator Solving the dynamics of arbitrarily driven time-dependent systems is, in general, quite demanding. The situation can be simplified if we allow the Hamiltonian to have a symmetry in time H(t + T ) = H(t). For this type of systems, the generator of the evolution over one period, namely U(t 0 + T, t 0 ) can be used to study the stroboscopic evolution of the system, i.e. at times t = t 0 + nT . The initial time t 0 can be regarded as a gauge freedom and we will set it to t 0 = 0 by convenience. We define U F = U(T, 0) as the Floquet generator, which will be of crucial importance in this work. Studying the reduced dynamics of the system, requires further tracing out the bath degrees of freedom. As a first approach to periodic driven-dissipative dynamics, we focus on a class of driving protocols for which this operation is analytically doable: the kicked protocols. Consider a constant system-bath Hamiltonian H, which we kick periodically with with a Hamiltonian H K that only contains system degrees of freedom. More precisely, The stroboscopic evolution of the full system-bath is governed by the Floquet propagator U F = exp(−igH K ) exp(−iH 1 T ). Hence, the state of the system after one period is given: where E(T ) is the CPTP map giving rise to the static one period evolution, and U K = exp(−igH K ) is the unitary kick operator. Therefore, for a kicked open system, we use the dissipation model for a time-independent problem. Note that, the Born approximation guarantees the state of the bath to remain unchanged after one period of the evolution. Therefore, we can iterate this process to obtain the state of the system at stroboscopic times t = mT with m ∈ N. The kicked protocol is used, for instance, in the models studied in [17,24,25].

Discrete time crystals beyond closed systems
As discussed in the introduction, time crystals in closed systems are identified by exhibiting: (I) discrete time-translational symmetry breaking, (II) rigidity on the subharmonic response of the order parameter and (III) the infinite persistence of the subharmonic response in the thermodynamic limit. Based in these conditions, we propose to characterize an open system time crystal using a single object, the Floquet propagator [17]: The spectrum {ε µ } of E F fully characterizes the stroboscopic properties of a periodically driven systems. Note that all ε µ lie within the unit circle. Focusing on: (i) The asymptotic space that we can now identify As(H) = span {Ψ µ ∈ Op(H) : |ε µ | = 1}, i.e. those eigenstates lying on the unit circle; and (ii) the dissipative gap ∆ F := max µ {1 − log |ε µ | : |ε µ | = 1}, the identification criteria for an open discrete time crystal can be formulated as: (I) Time-translation symmetry breaking: there exists at least one µ for which ε N µ = 1, but ε µ = 1. Note that if ε µ = −1 it comes with its conjugate pair.
(III) Persistence to infinite time: Strictly speaking, this property is ensured by (i). In general, for finite systems, 0 < ε N µ ≈ 1 which will fix the decay rate of the oscillations.
A system fulfilling (I)-(III) displays long-lived and robust subharmonic oscillations and therefore it is in a time-crystal phase. In the opposite situation, the system will typically arrive to equilibrium reaching the thermal phase.
Let us illustrate the typical scenario of a system in the time-crystal phase in the case of period doubling, i.e. N = 2. There should be, at least, two robust eigenvalues of E F with modulus one such that ε 2 µ = 1. We denote them by , ε 0 = 1 and ε π = −1. After a (large) number of periods m, the state of the system will be well approximated by The expectation value of the order parameter after m periods will be where we have used ⟪l 0 |ρ(0)⟫ = tr [ρ(0)] = 1. Therefore, we will be able to observe the long-lived subharmonic oscillations if: (a) the choice of ρ(0) is such that ⟪l π |ρ(0)⟫ = 0, i.e. the initial state has some overlap with the π eigenspace, and (b) the choice of O is such that ⟪O|r π ⟫ = 0, i.e. the order parameter is sensible to the subharmonic oscillations.

Useful observations for open system time crystals
Before proceeding further, let us remark some general properties of the Floquet propagator for open system time crystals. We focus on the dynamics consisting on a dissipative evolution under L 1 during a time t 1 followed from a unitary transformation U K (kick).
Observation 1 (kicked system propagator) The Floquet propagator E F of a system under a kicked protocol takes the simple form: and has always an eigenvalue ε 0 = 1. Proof It follows from Eq. (14) the form of the Floquet propagator. Notice that and therefore ⟪1| E F = ⟪1| has eigenvalue 1. Alternatively, it is a consequence of the fact that a concatenation of trace preserving maps is also trace preserving.
Observation 2 (unitarily transformed spectrum) The spectrum of a CPTP map E is unchanged under the unitary transformation

Observation 3 (kicked spectrum)
The spectrum of a CPTP map E is generically changed under a kick Observation 4 (asymptotic space of a kicked evolution) Consider {|ψ k } and {|φ k } two basis of H and a CPTP map E. If E is such that E |ψ k ψ k | = |φ k φ k | for a subset of tuples (k, k ) ∈ S, then it exists a unitary kick U K such that E F = U K ⊗ U * K E has at least |S| elements in its asymptotic space. Proof. Consider the unitary kick U K := k |ψ k φ k |. Then for all (k, k ) ∈ S, E F |ψ k ψ k | = |ψ k ψ k | and are therefore part of As(H) of the map E F . Note that the combinations of k and k that form S are not completely arbitrary. Since the map E is positive, the tuples (k, k ) with k = k can only be part of S if also (k, k) ∈ S and (k , k ) ∈ S.
Remark. Note that, in general, the dimension of the asymptotic space of a CPTP map E can increase, decrease or stay equal after a unitary kick.

Observation 5 (multistability and time-crystals)
A CPTP Floquet map E F supports time crystalline behavior only if dim As(H) ≥ 1. The converse is not true.
Observation 6 (general protocol for sub-harmonic response) A general kicked protocol on a CPTP map E with a unitary U K that gives rise to subharmonic response can be obtained by demanding: (i) the map E exhibits static multistability, i.e. the peripheral spectrum contains only eigenvalues ε µ = 1.
(ii) the unitary kick U K = U K ⊗ U * K acts independently on the asymptotic and decay spaces, i.e. U K = U As ⊕ U D .
(iii) the unitary kick U As has eigenvalues u α = exp(in α 2π/N ) with n α ∈ Z, and at least one eigenvalue is different from one.

Exemplary: Few-body systems
As previously shown, the structure of the asymptotic subspace As(H) is crucial to identify when a physical system can support subharmonic response. Here we present a collection of one and two qubit models (specified by H S and L α ) that give rise to different multistable structures of As(H). However, we do not refer to them as time-crystals since this multistability will be, in general, fine tuned.

A single qubit: Dephasing
A well known example of dissipative evolution supporting more than one steady state is pure dephasing. In this scenario, As(H) is able to support only one bit of classical information. We consider: which has dim As(H) = 2, and As(H) = span(1, Z) both with eigenvalue λ µ = 0. Note that there are two conserved quantities in this evolution, dual to the two steady states, that can be written as j 0 = 1 and j 1 = Z. Equivalently, one may think about alternative conserved quantities and steady states Ψ µ = j µ = |µ µ| with µ = 0, 1. Consider now U K = X, which can be implemented, for instance, as: such that H S (t + T ) = H S (t). Then E F = X ⊗ X exp(LT ) with the spectrum: The basis of As(H) for E F is again given by and bistability is lost. Since χ (1) = 1/2 1, we see that condition (II) is not fulfilled. Instead, consider the deformation of committing an error in the rotation angle: . The spectrum of E F (η), to lowest order in η, is given by which yields χ (1) = 0 and makes the subharmonic response linearly robust to rotation errors.

Two qubits: Local jumps
We have studied a one-qubit scenario where As(H) could support one classical bit. In this subsection, we consider the larger Hilbert space of two two level systems H = C 2 ⊗ C 2 . A new structure of dim As(H) that supports quantum information can now arise. These instances of asymptotic space are known as decoherence free subspaces, and have been recently studied in the literature as dissipation protected memories (see for instance [32,33]). We will consider two different subcases, and focus on the differences between those [34,35]. For convenience, we also introduce the Bell basis: We start considering no Hamiltonian and local noise operators: It is easy to see that As(H) = span(Ψ αβ ) with Ψ αβ = |ψ α ψ β | and, therefore, dim As(H) = 4. Note that we can encode the state of a qubit inside As(H) and it will be protected from dissipation. The conserved quantities read [31] for α, β = 0, 1. Two of those four conserved quantities correspond to the populations of α = β = 0 and α = β = 1, while the other two α = β correspond to the coherences only within the ψ-block. The coherences in the φ-space are destroyed by dissipation.

Two qubits: Collective jumps
If now we consider a collective noise operator L = α L α = α |ψ α φ α |, the same asymptotic space As(H) = span(|ψ α ψ β |) is found. The key difference can be spotted by looking at the conserved quantities. We encounter and therefore, also the coherences within the φ-space are now preserved. It turns out that preserving the coherences will become crucial in order to construct robust time-crystals.

Two qubits: Coherence suppression by a jump
We now look at the effect of the decay operators within the φ-space. Consider a jump L(η) = α |ψ α φ α | + η n · σ φ , where n is the Bloch vector and σ φ is the Pauli vector in the φ subspace (e.g.
This gives rise to a suppression factor for the coherences that depends on the particular choice of n. In general, For instance, for the particular choice of: indicates that dephasing within the decay space is translated into a coherence suppression of order ∼ |η| −2 within the φ -block.

Two qubits: Coherence suppression by a Hamiltonian
Finally, we will look at the effect of having a Hamiltonian that acts locally in the ψ and φ subspaces. This translates into having residual Hamiltonian evolution in the steady-state. Consider, and L = α |ψ α φ α |. The conserved quantities now read Hence, the effect of a coherent evolution within the asymptotic and decay space is to suppress the coherences between the φ-block. However, the strength depends on the effective detuning δ between the two spaces. Now, we have several possibilities of jumps operators, all of them giving rise to the same As(H) but with different conserved quantities and different eigenvalues λ µ . Regardless of the one we chose, a TC protocol can be implemented by kicking the system periodically with the rotation U K = X ψ + X φ . Note that U 2 K = 1 and U K = U As ⊕ U ⊥ , giving rise to sub-harmonic response for N = 2. It can be implemented by: For E F , the asymptotic space is As(H) = span(Ψ αβ ), with with eigenvalues ε αβ = (−1) α+β . The same linear combination gives rise to the conserved quantitiesj αβ . The order parameter O can be any observable that does not transform trivially under the not-gate X ψ , e.g. O = Z ψ fulfilling X ψ Z ψ X ψ = −Z ψ . In the local basis, this looks like O = |01 10| + |10 01| which is a coherence measure between the excitation of the two qubits. Hence, condition (I) is fulfilled. Regarding condition (II), exact diagonalization of the perturbed protocol is, in general, analytically demanding. However, if an error is made in the rotation angle, that is , the susceptibility can be computed. Using perturbation theory, one can verify In fact, this particular case can be computed at all orders to see that the spectrum for any η is {ε αβ } = {1, 1, −e iη , −e −iη }. A similar calculation, but with the generator of the kick operator given by M x = X 0 + X 1 , shows instead a robust rotation (i.e. χ (1) = 0). This last case, can be mapped to an open Ising chain of 2 particles that are kicked with the magnetization M x .

The open XY model
Consider a system described by the Hamiltonian: representing a 1D lattice of L spins that interact anisotropically. We impose periodic boundary conditions ( i.e. σ r = σ L+r ), restrict ourselves to L even, and gather the Hamiltonian parameters as ξ = (J, γ, h). This model has several symmetries [36]: (i) a rotation by π/2 along the z-axis interchanges the x and y spin interactions and it is equivalent to γ ↔ −γ, (ii) a reflection of the spins with respect to the x-y plane is equivalent to h ↔ −h. Hence, it is sufficient to study the phase diagram for γ, h ≥ 0. It is well known that this system undergoes a quantum phase transition at h = 1, that goes from an ordered phase for h < 1 to a disordered phase for h > 1. The full Hamiltonian of system and bath is given by: where k labels the potentially multiple reservoirs and η 1 is a small parameter ensuring the weak coupling limit. A similar system-bath Hamiltonian was used in [37] to study transport properties. The first step towards the derivation of the master equation is diagonalizing H ξ . This can be achieved using the Jordan-Wigner, Fourier and Bogoliubov transformation [36,38,39]. In the following, we sketch the diagonalization procedure both for completeness and to fix the notation used subsequently. First we note that H ξ can be broken in parity sectors since [P, H ξ ] = 0 and P = r Z r with possible eigenvalues p = ±1. With the Jordan-Wigner transformation, we map spins into fermions using: Note that this transformation does not depend on ξ. After some manipulation and imposing the appropriate boundary conditions in each parity sector we find: We now take advantage of the fact that H ± ξ are translationally invariant and perform the Fourier transform where BZ ± stands for the Brillouin zone specified by: This brings the Hamiltonian to the explicit block quadratic form We can now perform the ξ-dependent Bogoliubov transformation to diagonalize H ± ξ . Since the matrix H ξ,q is a combination of the Pauli matrices in the x and z directions, it can be diagonalized via a rotation such thatH ξ,q = R ξ,q H ξ,q R † ξ,q is diagonal. Setting the off-diagonal terms to zero requires: Finally, the Hamiltonian takes the expression: where ± stands for the even and odd parity sectors, d ξ,q are the Bogoulibov fermions and the dispersion is given by: Now we are ready to derive the master equation. We sketch here the crucial parts of the derivation, while details are given in App. E. For simplicity we assume only one reservoir, since the extension to multiple baths is analogous. The starting point is the Redfield equation in the rotating frame of H ξ + H B : The jump operators will arise from the decomposition of the M z into the eigenmodes of the system Hamiltonian. For every quasimomentum q, the magnetization can be divided into three different rotating frequencies with energies labeled by α: E α ξ,q = {0, ±2ω ξ,q } for α = 0, ↑, ↓ respectively. Then, the magnetization can be decomposed as: where L 0 ξ,q = 2 cos θ ξ,q (d † ξ,q d ξ,q −d ξ,−q d † ξ,−q ) and L ↑ ξ,q = 2 sin θ ξ,q d † ξ,q d † ξ,−q = L ↓ † ξ,q . Inserting this expression into Eq. (48), it followṡ The next step consists in using the secular approximation, that selects only those terms that fulfill the resonant condition E α ξ,q − E α ξ,q = 0. We now see that two different situations arise: If ξ is such that the dispersion of the energy as a function of q is approximatly flat (e.g. γ = 1, h = 0) the resonant condition gives α = α . If, on the other hand, the dispersion is large enough, only those terms with the same q will be resonant. These two conditions lead to collective and local decay processes in the sense of Subsec.4.3. Then, the final master equation for the collective decay processes is given by: where κ ξ are the decay rates and L ξ = q>0 L ↓ ξ,q . If, instead, the decay processes are local we haveρ Then, the jump operators destroy/create a pair of fermions of momentum ±q (see App. D), i.e. L ↓ ξ,q |m q ∝ m q |m q − 1 . At zero temperature, we have κ ↑ ξ = κ ↑ ξ,q = 0 and As(H) = span{|p, GS p , GS|} for p, p = ±1 which allows for a multistable region. In analogy to the examples of Sec. 4, it is the nature of the jump operators that discriminates the dynamics in Eq. (51) from that of Eq. (52).

The battle against decoherence
Before analyzing the stability of the open XY model as a time crystal, our aim is to identify when such a system will thermalize. We start noting that parity symmetry in the XY model gives rise to the four-dimensional asymptotic space As(H) = span{|p, GS p , GS|} with p, p = ±. In the ordered phase, the dynamics described by Eq. (51) and Eq. (52), periodically concatenated with the unitary kick U K = P, can display subharmonic oscillations. In particular, when the two ground states are perfectly degenerate, for any ρ ∈ As(H) we haveρ = 0, i.e. it is a steady asymptotic space. Consider the local order parameter O = m x = 1/L r X r and its expectation value m x = m x . Consider ρ ∈ As(H) with entries ρ pp , then This has an important implication: the oscillations are detected only when the asymptotic state retains some coherence between parity blocks. Then, when and why will coherence be lost? Consider the vectorized space H ⊗ H, and a partition between parity eigenblocks and coherences between them: where L ± are bona fide Liouvillians acting on the positive and negative parity blocks. It remains to check the action on the coherence part ρ c . We consider the case at zero temperature, the local master equation in Eq. (51) applied to a general coherence ρ c = |+, m −, m| (see App. D) yields giving rise to an exponential decay of the coherence of the excited modes. Note that the coherences in the ground subspace do not decay since m q = 0 ∀q. Therefore, subharmonic oscillations can be seen as long as the initial state has some coherence in As(H). Of course, this is true assuming U K = P. If we commit an error in the parity, the coherence in the asymptotic subspace can be mapped to coherence between some higher excited states. Then, the subsequent dissipation will decohere the state leading to a thermalized state in the (very) long time. Here is where we take advantage of the collective dissipation process. As we identified in Sec. 4, coherences can be preserved for collective decay processes. Hence, we expect the collective equation to be more robust against errors. In the following section, we analyze this possibility and compare the local and collective decay processes.

Characterization of the open-XY time crystal
In this section, we aim at characterizing the open-XY time crystal. We will do it combining analytic and numerical tools and will be primarily interested in evaluating properties (I)-(III) stated in Sec. 3. Of course, this analysis becomes costly very fast, as the number of spins L increases. Recall that the dimension of H . Therefore, we analyze the dynamics for relative small system size and expect that the results will not depend strongly on the system size.

Time translation symmetry breaking
We start by analyzing for which parameters ξ the response of the system is subharmonic. As it is shown in Fig. 2a, the subharmonic response survives to long times when the system is kicked properly, that is, with H K = π 2 n δ(t−nT )M z at any point of the phase diagram with h < 1. However, if the two ground-states are not exactly degenerate the oscillations may show beats. Hence, either along the factorization line (see App. C) for finite L or at any point of the ordered phase in the thermodynamic limit L → ∞ subharmonic oscillations are predicted to be observed (see Fig. 2).

Coupling to higher states via an error in the parity
Time crystals are expected to have some robustness against an error in the protocol. This behavior is corresponded with condition (II) when the -small-error in the protocol is generated by a too-long (π + η)-pulse. Instead of implementing the unitary rotation U K (0) = P, the implemented unitary rotation for small η reads: Therefore the quantum channel E F (η) will depend on the rotation error η. Expanding perturbatively the channel we obtain where we have defined V = −i/2[M z , E F (0)]. Therefore, we can get an idea of the robustness of this channel by evaluating the distance between the expected state after one evolution period for a perfect protocol and with an error η. As a figure of merit, we consider 1 − F(ρ(η, T ), ρ(0, T )), where F(ρ, σ) is the fidelity between quantum states, which upperbounds the trace distance squared. Taking as initial state the ground state |0, GS at every point of the diagram, and computing 1 − F(ρ(η, T ), ρ(0, T )) for η = 0.05π, we obtain the result shown in Fig. 2e. Because the ground states are invariant under the dissipative dynamics, the only effect comes from the error in the unitary part. We see that the error is smaller closer to the isotropic line γ = 0. This can be understood easily since, for γ = 0, the system has a continuous symmetry generated by the magnetization M z . Therefore, the action of V on eigenstates of the system will be negligible. On the other hand, as h grows, the ground states of the system become closer and closer to the magnetization eigenstate |GS ≈ |↑ ⊗L , which is invariant under any rotation generated by M z . In the rest of the parameter space, the error is approximately linear in η. Even though this seems to indicate that the time crystal is more robust around γ = 0, the most relevant feature for robustness is the collective nature of the jump operators. Hence, we show most of the results at the Ising point where the collective decay processes are encountered.

The phase diagram
The time crystal phase is often associated with a region of non-zero measure of the parameter space. In Fig. 2, we show the collective phase diagram for the open XY chain which exhibits two distinct phases of matter: a thermal phase and a time-crystal phase. The brightness of the heat map indicates the distance δ F = min µ |ε µ + 1|, that is, the spectral distance to the ideal subharmonic state.
In the thermal phase, the error produced in the kick operation induces decoherence in the system leading to a thermalized final state. When the baths are at zero temperature, thermalization leads to the statistical mixture ρ th = 1/2(|+, GS +, GS|+|−, GS −, GS|), for which oscillations in m x are no longer visible. Contrarily, in the time-crystal phase the collective dissipation protects the system from thermalization, leading to the observation of robust subharmonic oscillations on the order parameter m x . Ideally, the state on the time-crystal phase is given, after n periods, by ρ(mT ) = |k ⊕ m, GS k ⊕ m, GS| with k = 0, 1.

Preserving coherence for a Long-time TC
As we have briefly discussed in Sec. 5, the main process for which the oscillations in the time-crystal may die out is decoherence. We also have discussed in Sec. 4 that collective jump operators preserve coherence during the evolution. Hence, it is interesting to compare the robustness of collective and local master equations. In panel (c) of Fig. 2, we compare the dynamics evolution of the collective (blue) and local (orange) master equations for a system of L = 6 spins when an error η = 0.1 is committed during the rotation. As discussed previously, the oscillations are robust when the dissipation is generated by collective jump operators while it clearly decays when the jump operators are local.
In panel (d) we show the scaling with η of the distance δ F = min µ |ε µ + 1| for collective (blue) and local (orange) decay processes. The derivative of this plot at η = 0 corresponds to the susceptibility χ (1) .

Conclusions
In this article we have presented some self-contained results concerning the existence and properties of time crystals in open systems whose evolution is described with a Lindblad master equation. After introducing the tools of Markovian quantum open system dynamics, we have provided a compact definition of an open system time crystal derived from the properties the spectrum of the Floquet propagator. We have, as well, identified which are the most relevant properties of this object with special emphasis on the relevant aspects of the asymptotic subspace structure and the conserved quantities. We have analytically solved the kicked dynamics of an exemplary set of one and two-qubit open system models and exploit such analysis to provide key features on the properties and stability of time crystals in open systems. Finally, we have derived and analyzed the open XY model as a time-crystal. There has been some discussion around the possibility that only longrange models can exhibit time-crystalline order in open quantum systems, our analysis shows that this is not the case and we conclude that long-range interactions are not crucial features to observe time-crystalline behavior. Nonetheless, our findings show that longrange (or collective) jump operators are crucial in order to have subharmonic oscillations that are more robust to rotation errors. Intuitively, the collective jump operators help to preserve coherence in the dissipation process, and the time-crystalline oscillations are usually coherent in the Hamiltonian eigenbasis.
To conclude, we believe that a promising direction of investigation is that of non-Markovian environments, where the backflow of information to the system can be controlled to achieve sub-harmonic response.
A Mathematical properties of L For completeness, we include here some discussion about the mathematical properties outlined in Sec. 2 of the main text.
If a particular eigenvalue λ µ of L has algebraic multiplicity m µ ≥ 1, the number of non-trivial solutions of Eq. (6) lies between one and m µ . If there is strictly one solution |r µ ⟫ associated to λ µ but m µ > 1, higher rank generalized eigenvalues can be found as solutions of the recursive equation (L − λ µ I) |r µ (s)⟫ = |r µ (s − 1)⟫, where |r µ (1)⟫ = |r µ ⟫ and s denotes the rank. Note that (L − λ µ I) k |r µ (s)⟫ = 0 only if k ≥ s.
(i) Spectrum of the Liouvillian: Consider r µ ∈ Op(H) such that Lr µ = λ µ r µ . Then, hermiticity preservation Lr † µ = (Lr µ ) † guarantees that is, either the eigenvalues are real or come by conjugate pairs. Note that if r µ = r † µ λ µ ∈ R. The converse is true, at least, when λ µ is non-degenerate.
and, therefore, eigenvectors of different eigenvalues can be chosen biorthonormal. For a diagonalizable matrix, the biorthogonal relation can be compactly written as W ‡ l W r = I where the columns of W r are the eigenvectors |r µ ⟫.
(iii.2) Normal Jordan form: Given L, it exists a similarity transformation W r such that W −1 r LW r = J where J is in Jordan canonical form such that the columns W r are the generalized eigenvectors |r µ (s)⟫. This corresponds to solving the generalized eigenvalue equation LW r = W r J .
(iv) Time-ordered propagator: This is the well-known solution of a linear equation with time-dependent coefficients when the generator may not commute with itself at different times.
(v.1) Existence of the steady-state: The trace preserving condition for an arbitrary state ρ, together with Eq. (5) lead to: guarantees at least one eigenvalue λ 0 = 0. The corresponding eigenvector ρ ∞ = r 0 fulfills that ∂ t ρ ∞ = 0, and it is often referred to the steady-state. In general, however, the steady-state may not be unique.
(v.2) Contractivity of the evolution. Convergence to As(H): We include it here the proof given in [40]. Given an Hermitian operator where O ± are positive matrices, and a CPTP map E, it follows that In particular, given ρ, σ ∈ S(H) we see that the trace distance D(Eρ, Eσ) ≤ D(ρ, σ), which provides a convergence towards the asymptotic subspace.

B Generalized susceptibilities and higher order robustness
In the main text we considered the linear susceptibility χ (1) . However, higher order measures of robustness can be obtained as we show this in this section. Consider a general quantum map E = E 0 + ηE 1 with η 1. We aim at finding its spectrum defined via the equations: for a particular eigenvalue ε. We assume that r, l and ε can be expanded in powers of η It follows that Eq. (62) can be written which leads to the recurrence relation: The correction to the eigenvalues can be computed by projecting onto ⟪l 0 |, with the relation Perturbed eigenvectors can be also computed from Eq. (65), however its expression is quite involved and dependent on the choice of the inverse of the operator E 0 − λ 0 [41] and we do not include them here.

C XY chain: A sub-manyfold of product ground states
There exists a particular sub-manyfold ξ p of the parameter space ξ (within the ordered phase h < 1) such that the ground space of the system can be analytically found as product of rotated spin states. This manyfold is often refered as the factorization line. This sub-manyfold is described by the radius-one circle h 2 + γ 2 = 1, or equivalently, [36]. The exactly degenerated ground states are found: where cos 2 (2ζ) = (1 − γ)/(1 + γ). Note that this includes the Ising Hamiltonian at zero transverse field. Also, both states are connected via P |k, GS = |k ⊕ 1, GS where ⊕ means sum modulo two. This allows to define a ground state within each parity sector as such that P |p, GS = p |p, GS . The associated dispersion relation is for each of these ground states is:

D XY chain: Pseudo-spin representation
We start be rewriting the Hamiltonian of the system summed over only the positive quasimomentum part. From Eq. (46), we find where q > 0 = {1/2, · · · , (L − 1)/2}, a total of L/2 values; and q > 0 = {1, · · · , L/2 − 1}, a total of L/2 − 1 values, for the even and odd parity sectors respectively. Folding the BZ into the q > 0 part, the full Fock space for a given quasimomentum q is spanned by the four states |0 q , c † q |0 q , c † −q |0 q , and c † q c † −q |0 q . Within the vanishing total quasimomentum subspace, defined by Q = q qd † q d q , the states can be labeled with binary numbers collected in the vector m, such that with the subtlety that the q = −π mode should be unoccupied for the odd parity sector when L is even. This comes from the fact that the mode −q = q = −π goes into itself at the borders of the Brillouin zone. If the system is prepared in a state of the subspace of vanishing total quasi-momentum q, for instance |±, GS , a pseudo-spin representation is possible for each block q. This is because H ± ξ only connects the Fock vacuum (of the physical fermions) |vac with the state c † q c † −q |vac for each q. We introduce the notation: Then for any operator O = q O q , that acts independently on the different subspaces of quasi-momentum q, we can decompose it in this basis as: s|O q |s |s s .
In the second-quantization, the expression of O is given by: