Time crystallinity in dissipative Floquet systems

We investigate the conditions under which periodically driven quantum systems subject to dissipation exhibit a stable subharmonic response. Noting that coupling to a bath introduces not only cooling but also noise, we point out that a system subject to the latter for the entire cycle tends to lose coherence of the subharmonic oscillations, and thereby the long-time temporal symmetry breaking. We provide an example of a short-ranged two-dimensional system which does not suffer from this and therefore displays persistent subharmonic oscillations stabilized by the dissipation. We also show that this is fundamentally different from the disordered discrete time crystal previously found in closed systems, both conceptually and in its phenomenology. The framework we develop here clariﬁes how fully connected models constitute a special case where subharmonic oscillations are stable in the thermodynamic limit.

Introduction.Understanding how statistical mechanics emerges in closed quantum many-body systems undergoing coherent dynamics with time-independent Hamiltonians has been one of the major themes of physics research over the last few decades.More recently, attention has been focussed on closed systems with timeperiodic ("Floquet") Hamiltonians, where fundamentally novel out-of-equilibrium phases describable in macroscopic terms have been discovered; none more prominent than the π-spin glass also termed the discrete time crystal (DTC) [1][2][3][4][5][6][7][8][9][10].
Generically, a major obstacle to working with Floquet systems is that they suffer heat death due to unbounded increase of entropy, approaching an infinite-temperature state at long times [11][12][13].This heating can be avoided by introducing disorder-induced localisation [14,15], or by coupling the system to an external environment which drains energy from the system.The former approach, used in Ref. [1], additionally endows the Floquet eigenstates with a discrete-symmetry broken spatial glassy order.Crucially, the eigenstates connected by the symmetry (and having the same spatial ordering pattern) are separated in quasienergy by π/T with T the driving period, leading to the subharmonic oscillation of an appropriate local observable.It was shown in Ref. [16] that an external Markovian environment, unless explicitly finetuned, destroys such a delicate coherence required for the subharmonic oscillations; the system is driven towards a mixture of various Floquet eigenstates all with uncorrelated patterns of the spatial glassy order.
Here we analyse general dissipative Floquet systems in a different setting where the ordered phases, timecrystalline or otherwise, are stabilised by the dissipation, and which would be entirely absent without it.The mechanism, manifestly different from that of the π-SG, FIG. 1. Schematic of how a dissipative DTC can (a) survive and (b) die: Evolution of the distribution and the mean (black line) magnetisation of a quantum spin system with time.Within a period (of T ) a unitary process effecting global spin rotation and non-unitary cooling processes act together for tR whereas for tC , the cooling processes act alone.(a) and (b) show an example where the magnetisation distribution does not broaden in time.Consequently in (a) the weight is always on the correct magnetisation sector leading to the persistent subharmonic oscillations whereas in (b), the cooling processes split the distribution as such there is weight left behind in the wrong sector which eventually kills the DTC.
involves a periodic rotation between two "sectors" of a Hilbert space followed by dissipative "cooling" to states distinguishable by a measure such as magnetisation.
This intuitively appealing picture ignores the possibility that the dissipation in addition to the desired cooling also generates noise and eventual loss of phase coherence in the oscillations.In terms of the density matrix of the system, the noise potentially translates to the probability distribution of the observable over individual realisations being broad, which we argue destroys the DTC.We show that while this is indeed the case for 1-d spin chains with local interactions and dissipation, in 2-d the dissipative processes can naturally lead to a narrow distribution and hence persistent time-crystalline order unlike in 1-d.We note that this broadening mechanism can be absent altogether in mean-field dynamics such as for fully-connected models [17][18][19][20].
How a dissipative DTC can exist. . .First, let us describe the general arguments for the stability, or lack thereof, of DTCs to dissipation, before demonstrating them concretely using specific quantum spin-models.For a visual schematic, see Fig. 1.
• The Hilbert space is divided into different (for concreteness: two) sectors, H ± , which could be symmetry sectors or simply based on an empirical criterion based on the expectation value of an observable, and • each of these sectors has a manifold of states {|G ± } which posses quantum order characterised by an observable, M , which may for example be due to symmetry-broken order.Crucially, • the two ordered manifolds, {|G ± }, cannot be connected to each other via local operators, and • the expectation value m of the order parameter, M , is sufficiently narrowly distributed over the states within each of the manifolds that the two distributions for the two manifolds do not overlap (up to exponentially small corrections).
In Fig. 1, the two sectors are the positive-and negative-m halves of the vertical axis.Within this setting consider a two-step Floquet protocol in the presence of dissipative processes such that, • In the first step ("rotation phase"), the system evolves under the simultaneous action of a Hermitian rotation operator H R and the dissipative terms dissipation.In the absence of dissipation the evolution over the step is unitary and given by U R (θ), which maps neighbourhoods of the ground state manifold of one sector to states of the other sector and vice-versa.For a particular θ * , U R (θ * ) maps states from {|G + } exactly onto states from {|G − } and vice versa.
• In the second step ("cooling phase"), the system is governed by a Hamiltonian for which {|G ± } are ground state manifolds as well as by the same dissipation processes as in the first step.
The dissipative processes during both steps cool the system down so that under their sole influence all states in the H + sector would be driven to {|G + } and likewise for H − .Under the combined action of the unitary and dissipative terms, they take the system towards the |G ± of the sector it instantaneously finds itself in.
A quantum system driven with such a dissipative Floquet protocol, initialised in either of the ground state manifolds, shows a time crystalline response trivially if θ = θ * , as the expectation value of M oscillates stroboscopically between that in {|G + } and {|G − } with a period twice that of the Floquet drive, provided the dissipation is inactive in the rotation cycle.A certain robustness of the temporal order to deviations of θ from θ * is expected to be induced via the dissipation: if the unitary rotation does not take states from the ground state manifold of one sector (say {|G + }) entirely to that of the other sector ({|G − }), but admixes nearby excited states in the other sector, the cooling step of the drive can push the weight back onto the {|G − } manifold.Simply put, cooling kills off the excitations left behind by the imperfect rotation, stabilising a DTC.
. . .and how it can die: In a broad sense, the dissipative processes have three effects: • Hindering rotation during rotation phase Recall that {|G + } and {|G − } are not connected via local operators.U then naturally has the form of a global rotation of the degrees of freedom.If the rotation phase is not instantaneous, the state goes through excited states at intermediate times.However, the dissipation cools the system, opposing this creation of excitations hence making the rotation process less effective.Therefore the overall rotation with dissipation is less than without, trapping (part of) the weight in the wrong sector.This is unfavourable to the presence of a stable DTC.
• Correcting error caused by imperfect rotation during cooling phase Imperfect rotation potentially leaves the state in the correct sector, but not in the |G ± manifolds; dissipation corrects this, favouring the DTC.
• Broadening the distribution during both phases: The rotation and the cooling acting in conjunction can, for short times, increase the width of the distribution of the state's overlaps with the excited states such that the resulting state is spread over both the sectors.In the following cooling cycle of the Floquet drive, the weights in each sector can get pushed to their respective ground state manifolds, resulting in a finite weight in the wrong sector (see later discussion and Fig. 1).This is generally fatal to the DTC.
Of the three, the third (broadening) invariably causes the DTC signal to decay eventually.In its absence, when the probability distribution of the observable remains sharp over time, a stable DTC phase is possible with the first two mechanisms determining the parameter regime of the stability.Let us also note that while dissipation is favourable for the temporal order in the cooling cycle, it is detrimental in the rotation cycle, and it is a priori not obvious whether increasing the strength of the dissipation from some finite value favours or disfavours the temporal order.
In what follows, we introduce an explicit microscopic model and give three examples of dissipative processes.First we show that dissipation that cleanly separates the two sectors but broadens the magnetization distribution leads to a decay of the oscillation.We then introduce spatially local dissipation processes and show that: 1) in 1-d they fail to separate the two sectors, cause broadening, and lead to a decaying oscillation.In 2-d, they may cleanly separate the sectors and in addition do not result in broadening, so that in this case a stable DTC appears.
Quantum spin systems.To analyse the above ideas in a concrete setting, we consider a system of spins-1/2, first in 1-d.Using a basis constituted by the products states of σ z (which we henceforth denote as {|α }), the two sectors H ± can be taken to be the set of product states, . This is a natural choice for a system described by a ferromagnetic Ising Hamiltonian as in the limit of g → 0, {|α } is a possible set of eigenstates.Moreover, for g = 0, |G ± are adiabatically connected to the |⇑ (all-up) and |⇓ (all-down) states as long as the Hamiltonian is in the ferromagnetic phase, |g| < 1.Note that, defining the two sectors and the corresponding ground state manifolds in the fashion we do, also allows us to label the basis states and the sectors with the magnetisation density m = M /N (N being the system size).The unitary operator U R which in the thermodynamic limit maps states {|G + } ↔ {|α − } is given by U R (θ) = exp [−iθ σ x ] with θ ∈ (π/4, π/2] and is produced by the action of the Hamiltonian H R = θ t R j σ x j over time t R .It easily follows that for θ * = π/2, U R (θ * ) precisely maps the all-up state to the all-down state.In fact, since the ground state of the system breaks the Z 2 symmetry of the Hamiltonian spontaneously, U R (θ * ) connects the two ground states exactly throughout the ferromagnetic phase.As anticipated, the rotation U R (θ) is manifestly a non-local operation.In the thermodynamic limit, a product state with definite magnetization m is mapped to a (in the σ z -basis, non-product) state with definite magnetization m cos(2θ), so that the rotation operation does not result in broadening.
Depending on the dissipative processes involved, this model can show decaying (bottom of Fig. 1) or persistent (top) subharmonic oscillations depending on whether broadening occurs or not.In the following, we use three explicit Markovian dissipative processes to show how the absence (presence) of broadening due to them favours (disfavours) the persistence of the temporal order.
Lindblad dynamics.Focussing on Markovian dissipative processes, the equation of motion for the density matrix of the system is governed by the Lindblad equation where H(t) is the time-dependent (in our case, timeperiodic) Hamiltonian and {L i } is the set of timeindependent quantum jump operators which arise due to the coupling to the dissipative environment.Our binary Floquet protocol with period We focus mostly on the g → 0 limit of ĤTFIM .Direct jump operators.To demonstrate the deleterious effects of broadening we begin by considering a set of jump operators {L α } where m α denotes the magnetisation of the product state |α .These jump operators take the weight from any product state and transfer it directly to the ground state of the corresponding sector as well as causing exponential decay of offdiagonal elements of the density matrix in the product state basis.They therefore provide very efficient cooling.However, they lead to broadening of the distribution leading to the decay of the oscillatory signal.To show this explicitly, we study the time-dependent magnetisation of the system starting from the |⇑ state, making use of a simplification due to translational invariance to access very large system sizes [22].Fixing one of γt R and γt C and varying the other leads to the results shown in the two leftmost columns of Fig. 2: for either finite γt C , finite γt R , or both finite, the timecrystalline response of m decays exponentially with t while the lifetime of the subharmonic oscillations grows exponentially with γt C .On the other hand it decays at least polynomially with γt R , so overall stronger dissipation has opposite effects during each part of the driving.Nevertheless broadening of the distribution means the oscillations always decay via the mechanism shown in the lower panel of Fig. 1.
The magnetisation vanishes with time due to the state strobscopically being in a mixture of both ground state manifolds, with opposite magnetisations.
One then expects that an observable finite and equal in both the ground states will remain finite in a statistical mixture of the two, such as is obtained at long times in the present case.Such an observable is the correlator   5) is shown in the third column.For a fixed value of γtR (= 0.1 in (c)), the amplitude of the oscillations increases with γtC as dissipation corrects the error induced by imperfect rotationms.By contrast, for a fixed γtC (= 2 in (f)), the amplitude decreases with γtR as dissipation hinders the rotation.Here N = 51 and θ = 0.4π.
The persistent oscillations of C(t) are shown in the rightmost column of Fig. 2.This is a fundamental difference between this dissipative Floquet phase and the πspin glass, where persistent oscillations of C imply those of an initially finite m.The amplitude of the oscillations decreases with decreasing γt C because the signal is strongest in the two fixed point states |⇑ / ⇓ and γt C controls how well the system is cooled into the two ground states.In this sense, this order is also induced by dissipation.Increasing γt R causes C to remain close to unity (its value in the ground states) and resist rotation, consistent with the earlier general arguments.
The direct jump operators demonstrate that the broadening of the distribution in magnetisation is fatal to the time-crystalline order, even when the dissipation cleanly separates the two sectors.
Domain-wall annihilating jump operators.We now introduce a set of jump operators which avoid broadening in a natural way.These operators cause dynamics that only move domain walls (DWs): a free standing DW can move but not disappear.However two DWs can move into each other and annihilate.Such dynamics are fundamentally different in the 1-and 2-dimensional cases.
Denoting the neighbours of a site by {r }, to each product state |α and site there corresponds a jump operator where M {r },α is the net magnetisation of the spins in {r } and P i,s α i is a projector onto the spin at site i in spinstate s α i =↑ / ↓.The corresponding Lindblad dynamics along the diagonal is governed by a Pauli master equation with the γ determined according to the rules above while the off-diagonals decay exponentially.In 1-d, the dynamics along the diagonal amounts to the DWs executing a random walk, i.e. diffusing.The probability distribution of magnetisation starting from a sharp value m broadens at short times (Fig. 3 top left), and at long times becomes bimodal with two peaks at ±1 of height such that −p(−1) + p(1) = m (Fig. 3 top right), as found by solving Eq. ( 7) using a classical kinetic Monte Carlo approach.The resulting destruction of the DTC is shown in the bottom two panels.
For d ≥ 2, this type of domain wall dynamics eventually eliminates the minority phase by effectively causing a line (or surface) tension, tending to minimise the area of the interface between two non-conserved phases (see Fig. 4 for two examples of allowed transitions and Sup.Mat.[22] for a demonstration of how the dynamics minimises the interface length) so that the dissipation cleanly separates the two sectors.In general, local dissipative processes lowering the energy of (ferromagnetic) Ising-type hamiltonians in d ≥ 2 will behave in a qualitatively similar way.Less obviously, the dissipative dynamics does not broaden the magnetization distribution starting from a state sharp in magnetization, Fig. 5.It then follows that a DTC phase may be stable; this is supported by the lower panels of Fig. 5 where a rapid rotation is shown to result in persistent subharmonic oscillations while a slow rotation in a ferromagnetic phase in which the magnetization never changes sign.
In order to show this explicitly on finite-sized systems a numerical solution of the full Lindblad equation is desirable.However note that Eq.( 6) dictates that the number of jump operators is exponentially large in N , making a numerical solution all but impossible.We surmount this by noting that one can have another set of jump operators which lead to identical (to those of Eq.( 6)) dynamics for the diagonal elements of the density matrix but whose number grows polynomially with N .
Each 'majority rule' operator, for a given site, consists of a product of projectors onto the site's neigbouring spins (the set is over all possible configurations of P ↑/↓ , and the jump operators (6) shows an exponential decay of the timecrystalline order with polynomially decreasing lifetime with γtR.For the numerical solutions, θ = π/2 and γtC = 100. (a) Examples of the local jump operators leading to a persistent DTC in two dimensions: In general and for individual product states, these operators tend to either decrease the length of domain walls, such as the transition shown on the left, or cause minority regions to vanish, as on the right.
for the neighbours) multiplied by the spin raising (lowering) operator for the site if the projector configuration has fewer P ↓(↑) compared to P ↑(↓) , see e.g.Fig. 4. The resulting dynamics is displayed in Fig. 5 .Conclusions and outlook.We have discussed general mechanisms leading to dissipative (de)stabilisation of DTCs in disorder-free dissipative Floquet systems.Our example of a dissipation-stabilised DTC is completely distinct, relying on a fundamentally different mechanism, from the π-spin glass introduced in [1], which in turn is unstable to dissipation [16].We also uncover other, non- DTC but still dissipation-induced phases.Phenomenologically, a crucial difference between the dissipative Floquet system and the π-SG [1] is that in the former the oscillations in the magnetisation can decay even though its correlation function, C(t), can synchronise and oscillate persistently, while in the latter one implies the other.
While the mechanism of obtaining period-doubling from periodic switching between distinct sectors of Hilbert space is intuitively transparent, our analysis of its failure modes we believe also sheds light on recent work finding stable subharmonic oscillations [17][18][19].In these works the system Hamiltonian is fully connected.This typically leads to a stochastic description in which the noise vanishes with diverging system size, so that the master equation for the density matrix results in no broadening over the time evolution [23,24].
We believe that generally, treatments for short-range models based on approximate mean-field and other analyses involving only a few effective degrees of freedom may erroneously find stable time-crystalline behaviour in the absence of such a noise suppression mechanism.Our proposal is that the role provided by long-range interactions can, however, be replaced by the effectively macroscopic rigidity of the ordered component of a symmetry-broken system such as the Ising magnet subjected to the dissipative processes discussed here.
Finally, we have only considered Markovian dissipation.An open and interesting problem is to understand whether the physics unveiled here is changed qualitatively in the non-Markovian case, and whether there are non-fine-tuned non-markovian environments that lead to interesting new examples of oscillatory dynamics in quantum systems.

DOMAIN-WALL ANNIHILATING JUMP OPERATORS
In this section, we discuss further details of the domain-wall annihilating operators.The operators in Eq. ( 6) (main text) locally shift domain walls such that two of them can annihilate each other upon meeting.However the set of operators are rather inconvenient from the point of view of a numerical solution of the corresponding Lindblad equation as there are exponentially (in N ) many jump operators.We surmount this problem by considering a different set of jump operators which have the same effect on basis product states as the ones in Eq. ( 6) (main text) but which contains only polynomially in N of them.In particular, this new set contains 2 Z × N jump operators, (2 Z for each spin) where Z is the coordination number of the lattice and N is the number of spins.As mentioned in the main text, each operator in the set consists of a spin-flip (or lowering or raising) term for a given site, say , multiplied to projectors onto each configuration of its neighbours (hence 2 Z of them).If a neighbour-configuration has a net magnetisation which is positive (negative), the corresponding projector is multiplied to σ + (σ − ) and if the neighbour-configuration has equal number of up and down spins, the projector is multiplied to σ x , the latter with much smaller rate.
For a square lattice in 2-d, for a given site at position (x, y), there are three classes of jump operators: • All four of the neighbours are aligned with each other but anti-aligned with the spin at .There are two jump operators for this case, P ↑ x+1,y+1 P ↑ x+1,y−1 P ↑ x−1,y+1 P ↑ x−1,y−1 σ + x,y and the same with ↑→↓ and + → −.See Fig. 4 (right) [main text] for a visual example.
• Three of the neighbours are anti-aligned with the spin at and one is aligned.There are four configurations each for a the net magnetisation of the neighbours being positive or negative.So there are a total of eight operators in this class.As example, one of the operators in this class is P , and the projector on the down spin could be on any of the four neighbours which generates the other three operators.Similarly ↑→↓ and + → − generates the other four operators in this class.
• Two of the neighbours are up and two of them are down.
In this case, we always flip the spin at site , but the operator has a much smaller rate.An example operator in this class is P ↑ x+1,y+1 P ↑ x+1,y−1 P ↓ x−1,y+1 P ↓ x−1,y−1 σ x x,y There are six operators in this class (one for each of the twoup two-down configurations).
Note that the first two classes of the jump operators tries to reduce the length of the domain wall and favour the majority phase.However, they can get stuck if they encounter a straight domain wall.The third class of the jump operators serve to unfreeze such potentially frozen domain walls.The effect of these jump operators on the diagonal elements of the density matrix can be understood from a classical Monte Carlo simulation using Glauber dynamics with the local energy cost function for a spin at site given by σ z r∈{r } σ z r .While the Monte Carlo at zero temperature would try to generically take the system to the majority phase, there is a technical subtlety.The zero temperature classical Monte Carlo which tries to reduce the string tension of the domain wall can get stuck if it encounters a domain wall which straight; hence we run the Monte Carlo at a finite but small temperature.This essentially is the manifestation of the third class of the jump operators discussed above.
Fig. S2 shows a specific Monte Carlo trajectory with each panel showing the spin-configuration on the lattice at specific Monte Carlo times; green denotes regions of down-spin and yellow, up-spin.The evolution of the configuration clearly shows the shrinking of the domain-wall lengths and corresponds to a macroscopic coarse-grained version of the processes shown in Fig. 4 (main text).
As a final remark, we mention that the full solution of the Lindblad equation in the 2-d case (results of Fig. 5 (main text)) was obtained using the Monte Carlo wavefunction method, see [Mølmer et al., J. Opt.Soc.Am.B 10, 524 (1993)] for details.

FIG. 2 .
FIG. 2. Direct jump operators:The first column shows the instability of the subharmonic oscillations of the magnetisation density, m, for the jump operators (S4) for (a) γtR = 0.1 and different values of γtC and (d) γtC = 10 and different values of γtR.The respective lifetimes, τ , are shown in the second column, revealing an exponentially increasing and at least polynomially decreasing lifetime with γtC and γtR respectively.The data is for N = 51, and θ = 0.45π and θ = π/2 for the top and bottom rows respectively.Stability of the persistent oscillations of the correlation function C, Eq. (5) is shown in the third column.For a fixed value of γtR (= 0.1 in (c)), the amplitude of the oscillations increases with γtC as dissipation corrects the error induced by imperfect rotationms.By contrast, for a fixed γtC (= 2 in (f)), the amplitude decreases with γtR as dissipation hinders the rotation.Here N = 51 and θ = 0.4π.

FIG. 3 .
FIG. 3. Instability of the DTC to domain wall annihilating operators in 1D: (a) Results obtained from classical Monte Carlo in one dimension show that the distribution of the state in magnetisation for the jump operators in Eq. (6) broadens with time.(b) The mean magnetisation (dashed lines) stays constant over Monte Carlo times whereas the standard deviation (solid lines) grows, indicating broadening of the distribution.(c)-(d) Numerically solving the Lindblad equation with the time-periodic Hamiltonian, (3), and the jump operators(6) shows an exponential decay of the timecrystalline order with polynomially decreasing lifetime with γtR.For the numerical solutions, θ = π/2 and γtC = 100.

FIG. 5 .
FIG. 5. Stability of the DTC domain wall annihilating operators in 2D: (a)Results obtained from classical Monte Carlo in two dimensions show that the distribution of the state in magnetisation for the jump operators in Eq. (6) stays sharp.(b) The mean magnetisation (dashed lines) saturates to unity whereas the standard deviation of the system (solid lines) systematically goes down system size.(c)-(d) Numerically solving the Lindblad equation for square lattices (of size Nx × Ny) shows a persistent time-crystalline response of the magnetisation of for low γtR whereas an oscillating ferromagnet at high γtR For the numerical solutions, θ = π/2, γtC = 100, and γ = γ/10.

t
FIG. S2.Shrinking of domain walls in classical MonteCarlo: Evolution of the spin configuration on a 100×100 lattice obtained using classical Monte Carlo dynamics (Glauber dynamics) at inverse temperature, β = 2 visually shows the shrinking of the domain walls and eventual approach to a configuration completely taken over by the majority phase.Green regions denote spin-down and yellow regions, spin-up.