The fate of a discrete time crystal in an open system

TL;DR

This paper investigates the stability of discrete time crystals in disordered one-dimensional quantum systems when coupled to an external environment, showing that such coupling destroys the DTC phase over time.

Contribution

It provides analytical and numerical evidence that environmental coupling destabilizes discrete time crystals in one-dimensional disordered systems.

Findings

DTCs are destroyed at long times due to environmental coupling.

The destruction occurs even when the system is protected from heating.

Disorder does not prevent the loss of DTC order in open systems.

Abstract

Following the recent realisation that periodically driven quantum matter can support new types of spatiotemporal order, now known as discrete time crystals (DTCs), we consider the stability of this phenomenon. Motivated by its conceptual importance as well as its experimental relevance we consider the effect of coupling to an external environment. We use this to argue, both analytically and numerically, that the DTC in disordered one-dimensional systems is destroyed at long times by any such natural coupling. This holds true even in the case where the coupling is such that the system is prevented from heating up by an external thermal bath.

Equations

$$H(t)=\{\begin{array}{ll}{H}_z& \text{if }0\le t<T_z\\ H_x& \text{if }{T}_x\le t<T=T_z+T_x\end{array}$$

$$H_z$$

$$H_x$$

$$U:=U(T,0)=\exp \left(-i{H}_x{T}_x\right)\exp \left(-i{H}_z{T}_z\right)$$

$$|{\omega }_{\alpha };t⟩=U(t,0)|{\omega }_{\alpha }⟩$$

$$P=\prod _j{\sigma }_j^z.$$

$$U|{\omega }_{\alpha },\pm ⟩=\pm \exp (-i{\omega }_{\alpha }T)|{\omega }_{\alpha },\pm ⟩.$$

$$⟨\omega ,+|{\sigma }_i^x{\sigma }_j^x|\omega ,+⟩=⟨\omega ,-|{\sigma }_i^x{\sigma }_j^x|\omega ,-⟩\ne 0$$

$$|\uparrow /\downarrow ⟩=\left(|+⟩\pm |-⟩\right)/\sqrt{2}.$$

$$U|\omega ,\uparrow /\downarrow ⟩=\exp \left(-i\omega T\right)|\omega ,\downarrow /\uparrow ⟩$$

$$P|\omega ,\uparrow /\downarrow ⟩=|\omega ,\downarrow /\uparrow ⟩$$

$$\begin{array}{rl}⟨{\psi }_{2n}|{\sigma }_j^x|{\psi }_{2n}⟩& =⟨{\psi }_0|{\sigma }_j^x|{\psi }_0⟩\\ ⟨{\psi }_{2n+1}|{\sigma }_j^x|{\psi }_{2n+1}⟩& =-⟨{\psi }_0|{\sigma }_j^x|{\psi }_0⟩\end{array}$$

$${\rho }_0=\sum _{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}{\rho }_{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}|{\omega }_{\alpha },{\sigma }_{\alpha }⟩⟨{\omega }_{\beta },{\sigma }_{\beta }|$$

$$\rho_{0}=\begin{array}[]{c|c|c|c|c}\ddots&&&\\ \hline\cr&\mathbf{\tilde{\rho}}_{\alpha-1,\alpha-1}&\mathbf{\tilde{\rho}}_{\alpha-1,\alpha}&\mathbf{\tilde{\rho}}_{\alpha-1,\alpha+1}\\ \hline\cr&\mathbf{\tilde{\rho}}_{\alpha,\alpha-1}&\mathbf{\tilde{\rho}}_{\alpha,\alpha}&\mathbf{\tilde{\rho}}_{\alpha,\alpha+1}\\ \hline\cr&\mathbf{\tilde{\rho}}_{\alpha+1,\alpha-1}&\mathbf{\tilde{\rho}}_{\alpha+1,\alpha}&\mathbf{\tilde{\rho}}_{\alpha+1,\alpha+1}\\ \hline\cr&&&&\ddots\end{array}$$

$$\mathbf{\tilde{\rho}}_{\alpha,\beta}=\left(\begin{array}[]{cc}\rho_{\omega_{\alpha,\omega_{\beta}};\uparrow,\uparrow}&\rho_{\omega_{\alpha,}\omega_{\beta};\uparrow,\downarrow}\\ \rho_{\omega_{\alpha,}\omega_{\beta};\downarrow,\uparrow}&\rho_{\omega_{\alpha,}\omega_{\beta};\downarrow,\downarrow}\end{array}\right).$$

$${\rho }_n=\sum _{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}{\rho }_{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}{\mathrm{e}}^{\left(-i({\omega }_{\alpha }-{\omega }_{\beta })nT\right)}|{\omega }_{\alpha },(-1{)}^n{\sigma }_{\alpha }⟩⟨{\omega }_{\beta },(-1{)}^n{\sigma }_{\beta }|.$$

$$\mathrm{tr}\left({\rho }_n{\sigma }_j^x\right)={\left(-1\right)}^n\sum _{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}{\rho }_{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}⟨{\omega }_{\alpha },{\sigma }_{\alpha }|{\sigma }_j^x|{\omega }_{\beta },{\sigma }_{\beta }⟩\exp \left(-i\left({\omega }_{\alpha }-{\omega }_{\beta }\right)nT\right).$$

$${\rho }^{BDE}=\sum _{\begin{array}{c}\omega \\ {\sigma }_{\alpha },{\sigma }_{\beta }\end{array}}{\rho }_{\omega ;{\sigma }_{\alpha }{\sigma }_{\beta }}|\omega ,{\sigma }_{\alpha }⟩⟨\omega ,{\sigma }_{\beta }|$$

$$\rho^{\mathrm{BDE}}=\begin{array}[]{c|c|c|c|c}\ddots&&&\\ \hline\cr&\tilde{\mathbf{\rho}}_{\alpha-1,\alpha-1}&0&0\\ \hline\cr&0&\mathbf{\tilde{\rho}}_{\alpha,\alpha}&0\\ \hline\cr&0&0&\mathbf{\tilde{\rho}}_{\alpha+1,\alpha+1}\\ \hline\cr&&&&\ddots\end{array}$$

$$U{\rho }^{BDE}{U}^{†}=\sum _{\begin{array}{c}\omega \\ {\sigma }_{\alpha },{\sigma }_{\beta }\end{array}}{\rho }_{\omega ;{\sigma }_{\alpha }{\sigma }_{\beta }}|\omega ,-{\sigma }_{\alpha }⟩⟨\omega ,-{\sigma }_{\beta }|$$

$$\mathrm{tr}({\rho }_n^{BDE}{\sigma }_j^x)={\left(-1\right)}^n\mathrm{tr}({\rho }_0^{BDE}{\sigma }_j^x)$$

$${\partial }_t\rho =\mathcal{L}\rho$$

$$\mathcal{L}\rho =-i\left[H,\rho \right]+\sum _a\left(L_a\rho L_a^{†}-\frac{1}{2}{\left[L_a^{†}{L}_a,\rho \right]}_{+}\right)$$

$$\mathcal{L}(t)\rho =-i\left[H(t),\rho \right]+\mathcal{D}(\rho )$$

$$\mathcal{D}(\rho )=\sum _{\omega ,{\omega }^{\prime }}\left(L_{\omega ,{\omega }^{\prime }}\rho L_{\omega ,{\omega }^{\prime }}^{†}-\frac{1}{2}{\left[L_{\omega ,{\omega }^{\prime }}^{†}{L}_{\omega ,{\omega }^{\prime }},\rho \right]}_{+}\right).$$

$$\begin{array}{rl}{L}_{\omega ,{\omega }^{\prime }}(t)& =\sqrt{{\gamma }_{\omega ,{\omega }^{\prime }}}\sum _{\sigma }|\omega ,\sigma ;t⟩⟨{\omega }^{\prime },\sigma ;t|\\ & =\sqrt{{\gamma }_{\omega ,{\omega }^{\prime }}}|\omega ;t⟩⟨{\omega }^{\prime };t|\otimes \mathbb{I}.\end{array}$$

$${\hat{\rho }}_{{\omega }_{\alpha },\omega \beta }=\sum _{{\sigma }_{\alpha },{\sigma }_{\beta }}{\rho }_{\begin{array}{c}{\omega }_{\alpha },{\sigma }_{\alpha }\\ {\omega }_{\beta },{\sigma }_{\beta }\end{array}}|{\sigma }_{\alpha }⟩⟨{\sigma }_{\beta }|,$$

$${\partial }_t{\hat{\rho }}_{{\omega }_{\alpha },{\omega }_{\beta }}=-i\left[H(t),{\hat{\rho }}_{{\omega }_{\alpha },{\omega }_{\beta }}\right]+{\mathcal{D}}_{{\omega }_{\alpha },{\omega }_{\beta }}$$

$$\begin{array}{rl}{\mathcal{D}}_{{\omega }_{\alpha },{\omega }_{\beta }}& ={\delta }_{\alpha ,\beta }\sum _{\delta }{\gamma }_{{\omega }_{\alpha },{\omega }_{\delta }}{\hat{\rho }}_{{\omega }_{\delta },{\omega }_{\delta }}\\ & -\frac{1}{2}{\hat{\rho }}_{{\omega }_{\alpha },{\omega }_{\beta }}\sum _{\delta }\left({\gamma }_{{\omega }_{\delta },{\omega }_{\alpha }}+{\gamma }_{{\omega }_{\delta },{\omega }_{\beta }}\right)\end{array}$$

$${\partial }_t{\hat{\rho }}_{{\omega }_{\alpha }}=-i\left[H(t),{\hat{\rho }}_{{\omega }_{\alpha }}\right]+\sum _{ϵ}\left({\hat{\rho }}_{{\omega }_{ϵ}}{\gamma }_{{\omega }_{\alpha },{\omega }_{ϵ}}-{\hat{\rho }}_{{\omega }_{\alpha }}{\gamma }_{{\omega }_{ϵ},{\omega }_{\alpha }}\right)$$

$$L_j=\sqrt{\Gamma }{\sigma }_j^x$$

Full text

I Introduction

The field of non-equilibrium quantum many-body dynamics has seen increasing interest and rapid progress in recent years. Periodically-driven, or Floquet, many-body systems have been among the most rapidly progressing topics because of the promise of non-trivial long-time behaviour and the resulting recent experimental activity.Bordia et al. (2016); Zhang et al. (2016); Choi et al. (2016)

The main obstacle to observing interesting physics at long times in Floquet matter is that a general ergodic system simply heats up maximally under driving, as its entropy increases due to the non-adiabatic nature of the perturbation.D’Alessio and Rigol (2014); Lazarides et al. (2014a); Ponte et al. (2015a) This process may be frustrated either by tuning to certain integrable points for which the system does not heat up and its long-time behaviour is described by a so-called “periodic Gibbs ensemble” (PGE),Lazarides et al. (2014b) or more robustly by introducing disorder which leads to a many-body localised (MBL) phase in static systems.Basko et al. (2006); Oganesyan and Huse (2007); Pal and Huse (2010) Driving an MBL system leads to finite-energy-density long-time statesLazarides et al. (2015); Ponte et al. (2015b); Bordia et al. (2016) which may display non-equilibrium phases based on the notion of eigenstate order.Huse et al. (2013); Pekker et al. (2014)

One particularly intriguing example is the $\pi-$spin glass,Khemani et al. (2016) since also known as discrete time crystal (DTC) Else et al. (2016); Von Keyserlingk et al. (2016); Yao et al. (2017). The $d+1$-dimensional spatiotemporal order characterising DTCs manifests itself as a subharmonic response of the system to the driving. Like Bragg peaks signalling an increased magnetic unit cell compared to the structural one upon the onset of antiferromagnetic Néel order, this is encoded in the temporal Fourier transform.

This prediction has immediately sparked activity aimed at the experimental confirmation of this latest addition to quantum statistical mechanics.Zhang et al. (2016); Choi et al. (2016) The experiments found a temporally decaying order parameter. This immediately poses the question of the stability of the spatiotemporal order in realistic environmental conditions. In addition, in this setting there is the obvious question about the role of decoherence – a concept of perennial interest in quantum physics, playing a role in aspects as fundamental as the measurement process, and applications as important as quantum computing.Hone et al. (2009); Vorberg et al. (2013); Dehghani et al. (2014); Seetharam et al. (2015); Shirai et al. (2016)

We address these issues by studying the evolution of the density matrix via a Lindblad equation describing coupling to a Markovian environment. We find that at least in the one dimensional disordered systems recently investigated, the spatiotemporal DTC order is destroyed by any realistic model of environmental effects. We investigate two generic models of environmental decoherence as well as a model of an external thermal bath, finding that while the latter still destroys the DTC it nevertheless does not lead to an infinite-temperature state.

In detail, our results are the following. First, we show that DTCs can be described, in the long time limit, by a form of the density matrix that we call block-diagonal ensemble (BDE), which has a block-diagonal form in a basis we define related to the Floquet basis and which is manifestly $2T$ periodic. For any local operator to be sensitive to this periodicity, disorder must be present and the initial states must themselves break the $\mathbb{Z}_{2}$ symmetry of the drive. We then discuss which properties lead to a Lindblad-type equation which has this form as a steady state. Because DTC is a property of the Floquet states, not detectable by spatially and temporally local measurements, we are led to conclude that Lindblad operators preserving DTC are most naturally expressed in the Floquet basis and have the property either of not coupling different Floquet states to each other at all, or of coupling Floquet states only to others selected based on the local operator that displays oscillations. This leads to the conclusion that coupling to physical environments typically destroys DTC, since in general the environment will not have this property. Thus while it is possible to write down Lindblad operators preserving the DTC, these do not appear to correspond to any natural physical processes.

The remainder of this paper is organised as follows. We first review DTCs in isolated Floquet systems in order to define the problem and fix notation. We then set up our analysis for the presence of dephasing non-unitary dynamics. We present our results on the timescales over which the DTCs persist as well as our general conclusions on the stability of the DTC in open one dimensional systems, based on an analysis of different types of coupling models to the environment. We conclude with an outlook.

II Discrete time crystals: Unitary evolution

We set up the problem by taking the first and arguably simplest model of a DTC as our object of study. For a brief review see Ref. Moessner and Sondhi, .

II.1 Model and phenomenology in the $\pi$-SG phase

The Floquet dynamics is provided by a binary drive, in which the period $T$ is subdivided into two parts during each of which a time-independent Hamiltonian generates unitary evolution:

[TABLE]

The model describing a chain of spins-1/2 represented by Pauli matrices $\sigma_{i}^{\gamma}$ isKhemani et al. (2016)

[TABLE]

Here $J_{i}=1+\delta J_{i}$ with $\delta J_{i}$ randomly drawn from a uniform distribution between $-\delta J/2$ and $\delta J/2$, while $h_{i}=\overline{h}+h_{i}$ with $h_{i}$ randomly drawn form a uniform distribution between $-w/2$ and $w/2$.

As a technical aside, we note that in the numerics that follows we set $J_{z}=0$ unless otherwise indicated as it has been shownKhemani et al. (2016); Von Keyserlingk et al. (2016) that the DTC phase is stable to the introduction of interactions. However the numerically accessible range of $J_{z}$ can be quite small especially as our exact approach only allows for small system sizes for which a flip of a single spin’s $z$ component with concomitant energy change $J_{z}$ can already amount to considerable energy “density.”

The unitary operator

[TABLE]

propagates the system over one period. Its eigenvectors $\ket{\omega_{\alpha}}$, $U\ket{\omega_{\alpha}}=\exp\left(-i\omega_{\alpha}T\right)\ket{\omega_{\alpha}},$ are the Floquet states which for stroboscopic dynamics play a role analogous to that played by energy eigenstates for static systems. The eigenvalues are of the form $\exp\left(-i\omega_{\alpha}T\right)$ with $T=T_{x}+T_{z}$ the period and the quasienergies $\omega_{\alpha}$ real. For later convenience let us also define

[TABLE]

which is useful in describing in-period dynamics.

The model of Eq. (1) supports a number of phases; in this paper we shall be exclusively concerned with its time-crystalline $\pi$-SG phase, in which the eigenstates of $U$ have the following salient properties:Thakurathi et al. (2013); Benito et al. (2014); Khemani et al. (2016); Von Keyserlingk et al. (2016)

Are eigenstates of the parity operator

[TABLE] 2. 2.

Come in pairs of opposite parity, separated by quasienergy $\pi/T$ (up to corrections exponentially small in the system size) and forming a doublet for each $\omega$. They can be labelled as $\ket{\omega_{\alpha},\pm}$ with

[TABLE]

The doublet can also be thought of as a pseudospin-1/2 for each $\omega$ block. 3. 3.

Display spatial long-range order: the expectation value of the operator $\sigma^{x}_{i}\sigma^{x}_{j}$ is finite in all eigenstates for arbitrarily large $|i-j|$ and is the same for both members of each doublet.

[TABLE]

It is however different for different $\omega$, and does not depend smoothly on $\omega$. 4. 4.

The matrix elements of the operator $\sigma^{x}_{j}$ are finite between the two members of a doublet as shown in Fig. 1.

The last point above indicates that the dynamics of this operator will exhibit subharmonic oscillations. To bring this out more clearly, let us rotate the doublet basis $\ket{\pm}$ to a broken symmetry basis,111Note that $\uparrow/\downarrow$ do not refer to physical spins; they simply label the two many-body states.

[TABLE]

These two basis states are eigenstates neither of the parity operator $P$ nor of the unitary operator $U$ and have the following properties:

The expectation value of $\sigma^{x}_{j}$ is generically finite in each of these. 2. 2.

From Eq. (6) it follows that the action of $U$ is to flip the spin and produce a phase:

[TABLE] 3. 3.

From Eq. (8) it follows that the action of the parity operator is to flip the spin:

[TABLE] 4. 4.

From the two last points one concludes that the broken symmetry states break spatiotemporal symmetry, with a subharmonic response of period $2T$: if $\ket{\psi_{0}}=\ket{\omega_{\alpha},\uparrow}$ and $\ket{\psi_{m}}=U^{m}\ket{\psi_{0}}$ for integer $m$ then

[TABLE]

III The long-time limit and the block-diagonal ensemble

In this section we use the above properties to show that the long-time state of the system is well described by a density matrix of a special form that we call the “block-diagonal ensemble” (BDE). This form follows naturally from the doublet structure and the form of the matrix elements of operators breaking the $\mathbb{Z}_{2}$ symmetry (Fig. 1) and generalises the diagonal ensemble occuring in systems with no special spectral structure.Russomanno et al. (2012); Lazarides et al. (2014b); Moessner and Sondhi

Let us now consider the form of the BDE matrix in detail. An initial density operator with matrix elements $\rho_{\alpha\beta}=\bra{\omega_{\alpha}}\rho\ket{\omega_{\beta}}$ after $n$ periods becomes $\rho(nT)=\sum_{\alpha,\beta}\rho_{\alpha\beta}\exp\left(-i\left(\omega_{\alpha}-\omega_{\beta}\right)nT\right)\ket{\omega_{\alpha}}\bra{\omega_{\beta}}$. For a generic (ergodic, non-MBL, non-spectrally-paired) system, a) the eigenvalues $\omega_{\alpha}$ are continuously distributed with no special spectral structure apart from repelling each other due to the ergodic nature of the system and b) local operators have matrix elements that are maximal near the diagonal $\omega_{\alpha}=\omega_{\beta}$, as for static systems.Beugeling et al. (2015) This leads to such operators having synchronised states which are described by the so-called diagonal ensembleRigol et al. (2008) in the floquet basis, in which terms with $\alpha\neq\beta$ do not contribute at long times. The result is that the long-time steady-state becomes independent of the period $n$, thus synchronised with the driving.Russomanno et al. (2012); D’Alessio and Rigol (2014); Lazarides et al. (2014b); Ponte et al. (2015a)

In the $\pi-$SG case, by contrast, a) the eigenvalues are continuously distributed except for the doublet pairing structure, and b) there exist operators the matrix elements of which are appreciable between pairs of states separated by quasienergy $\omega/2=\pi/T$. If the operators of interest possess the latter property, terms off-diagonal in the Floquet basis and differing by quasienergy $\pi/T$ are the only ones that survive leading to what we call the “Block Diagonal Ensemble” (BDE). Density matrices of this form describe a steady state which is periodic with period twice that of the driving (rather than with the same period as the driving).

Concretely, in terms of the broken-symmetry states a general initial density matrix may be written as

[TABLE]

which may be visually represented as

[TABLE]

with the matrices

[TABLE]

Time evolved (see Eq. (6)) this becomes $\rho_{n}=U^{n}\rho(U^{\dagger})^{n}$,

[TABLE]

The expectation value of $\sigma_{j}^{x}$ after $n$ periods is

[TABLE]

That is, time evolution affects parts of this sum corresponding to the diagonal blocks ($\alpha=\beta$) of the representation in Eq. 13 only via a period-$2T$ flipping of the pseudospin, while all the other blocks acquire an additional time-dependent phase.

These phases are not correlated with each other, nor are they commensurate with the driving. Noting that all diagonal blocks give a $2T$-periodic contribution to the sum, while off-diagonal blocks also contribute a “random” frequency (incommensurate with that of the driving and the other blocks), we conclude (as earlier) that the off-diagonal blocks may be set to zero in the long-time limit for the same reason the diagonal ensemble holds in the absence of the doublet structure. This then leads to the BDE,

[TABLE]

visually represented as

[TABLE]

The action of $U$ on this is only to flip $\uparrow/\downarrow$ to $\downarrow/\uparrow$, that is,

[TABLE]

Hence, provided that $\rho_{\omega_{\alpha},\omega_{\alpha};\downarrow,\uparrow}\neq\rho_{\omega_{\alpha},\omega_{\alpha,};\uparrow,\downarrow}$ or $\rho_{\omega_{\alpha},\omega_{\alpha};\downarrow,\downarrow}\neq\rho_{\omega_{\alpha},\omega_{\alpha};\uparrow,\uparrow}$ (or both), a density matrix of the form Eq. 16 has period $2T$, since $U^{2n}\rho^{\mathrm{BDE}}(U^{\dagger})^{2n}=\rho^{\mathrm{BDE}}$ for any integer $n$ while $U^{2n+1}\rho^{\mathrm{BDE}}(U^{\dagger})^{2n+1}\neq\rho^{\mathrm{BDE}}$. In particular, we note that the expectation value $\sigma_{j}^{x}$ in the BDE is

[TABLE]

from which it follows that the spin oscillates with period $2T$ in the BDE.

For finite system sizes $L$, the two partner states are not separated by quasienergy exactly $\pi/T$ but rather by $\pi/T+\epsilon$ with $\epsilon$ exponentially small in the system size, and different for each pair of states. After a time $\sim\exp(L)$, which diverges in the thermodynamic limit, the oscillations are therefore washed out and the DTC disappears.Von Keyserlingk et al. (2016)

IV Lindblad operators and the BDE form

The most general equation for the possibly non-unitary evolution of a quantum density matrix is given by the Lindblad equation (see Appendix A)

[TABLE]

with the Lindblad operator

[TABLE]

where $\left[\cdot,\cdot\right]_{+}$ is the anticommutator and the $L_{a}$ arbitrary operators encoding the non-unitary part of the dynamics. In what follows, we use this form with the time-dependent Hamiltonian of Eq. (1) and choices of Lindblad operators $L_{a}$ appropriate to different types of environmental couplings.

We begin by discussing which features the operators $L_{a}$ must display in order to not destroy the BDE structure responsible for the subharmonic oscillations. As the DTC is a phenomenon that cannot be detected by local measurements at a single point in time, the Lindblad operators that preserve DTC are naturally written down in terms of Floquet operators. We will identify properties required for the DTC to survive and argue that these are unlikely to appear in real systems.

The Lindblad equations we consider are of the form $\partial_{t}\rho=\mathcal{L}\left(t\right)\rho$ with

[TABLE]

with

[TABLE]

We take the Lindblad operators to be of the general form

[TABLE]

Explicitly, defining

[TABLE]

we have

[TABLE]

with

[TABLE]

which manifestly has steady states of the BDE form of Eq. (15): Intra-block dynamics is only generated by the unitary term, which as discussed in Sec. III results in the subharmonic features, while inter-block dynamics are only generated by the $\mathcal{D}$ term. The time evolution due to the latter couples the diagonal blocks ($\alpha=\beta$) to each other and only to each other, while each off-diagonal block is not coupled ($\alpha\neq\beta$) to any other blocks so that each simply decays exponentially with time.

IV.1 No mixing between blocks

We begin with a DTC-preserving choice for the Lindblad operators of (20), namely, $\gamma_{\omega,\omega^{\prime}}=\gamma_{\omega}\delta_{\omega,\omega^{\prime}}$ so that $L_{\omega,\omega^{\prime}}(t)=\sqrt{\gamma_{\omega}}\delta_{\omega,\omega^{\prime}}\ket{\omega;t}\bra{\omega;t}\otimes\mathbb{I}$ is a projector onto a doublet. For this projector form $\mathcal{D}_{\omega_{\alpha},\omega_{\alpha}}=0$ while for $\alpha\neq\beta$ we have $\mathcal{D}_{\omega_{\alpha},\omega\beta}=-\frac{1}{2}\mathbf{\hat{\rho}}_{\omega_{\alpha},\omega_{\beta}}\left(\gamma_{\omega_{\alpha}}+\gamma_{\omega_{\beta}}\right)$. Therefore the entire time dependence is exponential damping of all off-diagonal blocks (with rates determined by the various decoherence channels) while the diagonal ones are completely unaffected by the evolution. Consequently, any initial density matrix with pseudospin imbalance (thus an initial finite expectation value of $\sigma^{x}_{j}$) will result in a BDE long-time density matrix displaying subharmonic oscillations.

While we have succeeded in constructing Lindblad operators leading to the BDE form and subharmonic oscillations, we still need to discuss how realistic such a choice of operators is. To answer this, it is useful to think of the process described by Eq. (18) as the system evolving unitarily, with the unitary evolution interrupted at random and with a typical rate set by the $\gamma$, by measurements described by the operators $L_{\omega,\omega^{\prime}}$.Dalibard et al. (1992); Carmichael (1993) For the environment to have the effect described here it must effectively measure the occupancies of doublets. This is not a natural operation for two reasons: Firstly, a measurement of a projector onto a doublet is impossible to achieve with local operations as it would require a simultaneous measurement at all points in space; secondly, the projector is time-dependent, in a way dictated by the unitary time evolution of the Floquet system itself. Such an external environment therefore seems very fine-tuned and unlikely to appear naturally, even if it is not physically forbidden.

IV.2 Mixing between blocks

One may ask whether it is possible to remove the restriction of no mixing between the blocks and have DTC still survive. Here we show that the answer is yes for certain fine-tuned choices of Lindblad operators; for general forms of the $\gamma_{\omega_{\alpha},\omega_{\beta}}$ the DTC is destroyed.

To lighten the notation we use the fact that the off-diagonal blocks evolving according to Eq. (22) simply decay exponentially, so the long-time BDE state may be obtained by solving the equation

[TABLE]

with $\mathbf{\hat{\rho}}_{\omega_{\alpha}}=\mathbf{\hat{\rho}}_{\omega_{\alpha},\omega_{\alpha}}$ (see Eq. (21). As before, the unitary (first) term only produces intra-block dynamics by flipping the pseudospin periodically; the interblock dynamics is described by a master equation for the $\mathbf{\hat{\rho}}_{\omega}$. The steady state is determined by setting the sum on the right hand side to zero and is automatically of the BDE form. For general $\gamma$ that do not separate the state space into disjoint sets, such master equations generally result in a steady state in which all the blocks have a finite occupancy, with the occupancies determined by the details of the $\gamma$. This will in general result in vanishing expectation values for $\sigma^{x}_{j}$ due to the fact that its expectation value varies randomly between blocks and is not correlated with the pseudospin; this follows from the discussion around Eq. (7) in section II.1. We now discuss when this response might not vanish, finding that the conditions again correspond to environments with very unnatural properties.

First, there will be a subharmonic response if the steady state corresponds to a single occupied block such that there is pseudospin imbalance. A single-block steady state will occur if the $\gamma$ are such that there is a state into which there are transitions but out of which there are no transitions. This would require the effect of the environment to favour one particular floquet state over all others, corresponding again to an environment performing fine-tuned, time-dependent operations everywhere in space simultaneously.

Second, a subharmonic oscillation of the operator $\sigma^{x}_{j}$ for a given $j$ will also result if the $\gamma$ only connect blocks for which, for the same pseudospin orientation for each block, $\sigma^{x}_{j}$ has the same sign. In this case the sign of the spin will flip in time, while its magnitude may decrease but will not vanish. However, for a finite fraction, say $\nu$, of the spins to oscillate, the $\gamma$ must only connect a fraction $2^{-\nu L}$ of the blocks which vanishes in the thermodynamic limit. Thus the environment would need to be selective in which states it couples, which is the same problem as in our earlier attempts to construct DTC-preserving operators. In addition this scheme suffers from all the drawbacks mentioned earlier, namely, that the environment would need to be performing finely-tuned time-dependent and non-local measurements on the system.

V Physically motivated Lindblad operators

We now turn to the direct numerical solution of the Lindblad equation for three commonly studied and physically realizable types of Lindblad operators. Two of them are dephasing operators, corresponding to generic interactions with the environment that destroy quantum coherence and have no particular energy structure. The third are “thermal” operators, which in the static case lead to thermal (Gibbs) long-time states. In all three cases we find that DTC is destroyed, in agreement with the arguments of Sec. IV, since these operators do not have the special structure required to preserve the BDE form.

In what follows we set $h_{i}$ and $T_{z}$ such that the average $h_{i}$ is $\overline{h}$ with $\overline{h}T_{z}=\pi$. In all cases we take as initial state an equal superposition of the ground and first excited states of an Ising Hamiltonian in the ferromagnetic phase, with $J=1$ and $h=0.2$. This ensures that $\sigma_{1}^{x}$ has a finite expectation value (the site $i=1$ is arbitrarily selected–our results are independent of this choice).

V.1 Dephasing

We start by looking at dephasing operators, relevant to trapped ion and cold atom experiments as well as to the experiment reported in Ref. Choi et al., 2016. Physically, they model the effect of an external environment that performs projective measurements of the spin in some direction. Alternatively, they model environmental effects without a preferred energy scale. One example is the non-unitary dynamics generated by incoherent scattering of the lattice laser light in cold atom systems.Pichler et al. (2010)

The two cases we study differ in that one preserves the parity of the initial state while the other does not; however, as we will see, the DTC is destroyed in both cases.

V.1.1 Parity-violating dephasing: $x$ direction

The first type of dephasing operator we consider is

[TABLE]

for $j=1,2,\ldots,N$, aligned along the ferromagnetic direction of Eq. (1). Real-time (as opposed to stroboscopic) results are shown in Fig. 2. The figure shows $\sigma^{x}$ at an arbitrarily selected position for fixed values of $\Gamma$ and $\bar{h}T_{z}$ and varying $T_{x}$ and $T_{z}$.

The Lindblad operators of Eq. (24) cause decay of density matrix elements off-diagonal in the $\sigma^{x}$ basis and do not preserve parity. In the absence of a unitary part (ie for $H=0$), any product state of eigenstates of $\sigma^{x}_{j}$ is a steady state as is a statistical mixture of such states; while in the presence of a generic, time-dependent $H$ the single steady state is the identity matrix and the steady state is the fully mixed state $\rho\propto\mathbb{I}$.

The right panel of Fig. 2 shows the Fourier transform of the evolution, confirming that

• •

In both cases there is a peak at $1/2$ the driving frequency, and

• •

Slower driving does result in stronger damping, broadening the peak of the transform.

We quantify the decay rate by fitting the quantity $\left|\left<\sigma_{1}^{x}\right>(t)\right|$ with an exponentially decaying function, ie, by determining $\xi_{t}$ in $\left|\left<\sigma^{x_{1}}\right>(nT)\right|=\sigma\left(0\right)\exp(-\xi_{t}nT)$. The results of this fit are shown in Fig. 3: The left panel shows the dependence of $\xi_{t}$ for fixed $\Gamma$ and varying $T_{z}$, while the right panel shows the same quantity for varying $T_{z}$ but now keeping $\Gamma T_{z}$ fixed, demonstrating that it is the quantity $\Gamma T_{z}$ that determines the damping rate rather than $\Gamma$ or $T_{z}$ alone.

A simple picture for this is provided by a minimal model for this type of driving for a single spin presented in Appendix B. There it is shown that if the flipping (paramagnetic) part acts for time $\tau$, the Lindblad operator is proportional to some $\gamma$ and the transverse field is $h$ such that $h\tau=\pi$ (thus in the absence of the dephasing term would exactly flip the spin) then the $x$ component after time $\tau$ is $\left\langle\sigma^{x}\right\rangle\left(\tau\right)=\mathrm{tr}\left(\rho\left(\tau\right)\sigma^{x}\right)=\rho_{x}(\tau)$ with

[TABLE]

so that

• •

the spin loses polarisation; and

• •

the rate of loss of polarisation (damping rate) depends on the product $\gamma\tau$ rather than on each factor individually.

These properties we also find for our many-spin problem, where the rate of damping only depends on the product $\Gamma T_{z}$ as shown in the right panel of Fig. 3, suggesting that the toy model of the appendix correctly describes the relevant physics for this type of Lindblad operator. In particular this shows that the more rapid the flipping of the spins, the longer the lifetime of the DTC for a given damping rate $\Gamma$.

V.1.2 Parity-preserving dephasing: $z$ direction

The second type of Lindblad operator we consider is

[TABLE]

causing dephasing in a direction aligned with the transverse field. Under the action of this type of dynamics, parity is a conserved quantity: as follows from Eq. 30, $\mathrm{{tr}\left(P\rho\left(t\right)\right)}=\mathrm{{tr}\left(P\rho\left(0\right)\right)}$ for all $t$, with $P$ the parity operator of Eq. 5.

Real time and Fourier transformed results for this are shown in Fig. 4. The difference from the case of V.1.1 is that now the decay occurs throughout the evolution, rather than mostly during the “flipping” $L_{z}$ part, therefore there are no longer plateaus in the expectation value of $\sigma^{x}_{j}$. Otherwise the overall behaviour is very similar to the previous case: The DTC is destroyed. This is a concrete demonstration that, as follows from the arguments of IV.2, it is not enough that the symmetry is preserved for the survival of the DTC phase.

V.2 Finite-temperature bath

As the operators of both Eq. (24) and Eq. (26) are Hermitian, the fully-mixed density matrix $\rho=\mathbb{I}$ is a solution. This state corresponds to infinite temperature (see Fig. 5), and is not unexpected since these Lindblad operators cause transitions between states without any energy preference. To show that this is not what destroys the DTC, we now select a set of Lindblad operators that in the static limit would have a finite-temperature Gibbs state as the steady-state. Under periodic driving, there are two processes at play: the unitary evolution, which tends to increase the energy of the system, and the dissipative evolution which tends to take the system to some finite-energy state. While their interplay and resulting steady-state are complicated, we find that the long-time state is not infinite temperature but still does not display DTC.

In the static case, a thermal bath is often modelled by coupling the system of interest to an environment with a thermal energy distribution. This is done via a term in the system-bath Hamiltonian $H_{SB}=V\otimes B$ with $V$ acting on the system degrees of freedom and $B$ acting on the bath degrees of freedom. Under a set of broad assumptions, one arrives at a Markovian master equation for the system of the form of Eq. 18 with Lindblad operators of the form

[TABLE]

with $\ket{m}$ eigenstates of the instantaneous Hamiltonian, $H\ket{m}=\epsilon_{m}\ket{m}$ and $\Gamma_{mn}=2\pi\left|V_{mn}\right|g\left(\epsilon_{n}\right)$. Here $g\left(\epsilon_{n}\right)$ models everything about the external thermal bath and $V_{mn}=\bra{m}V\ket{n}$ is a matrix element of the system part of the bath-system coupling operator in the instantaneous eigenstates. If this operator commutes with the parity operator, $\left[V,P\right]=0$ and the Hamiltonian also satisfies $\left[H,P\right]=0$ then the matrix elements between eigenstates of different parity vanish so that the parity a conserved quantity. The choice $g\left(\epsilon\right)\propto\exp\left(-\beta\epsilon\right)$ ensures that the steady-state of this static system is the Gibbs distribution $\rho_{S}=\sum_{n}\exp\left(-\beta\epsilon_{n}\right)\ket{n}\bra{n}$ while the time evolution up to the steady-state depends on the specific form of the operator $V$, that is, on the details of the system-bath coupling.

We shall apply this approach to the Floquet case by taking the bath to act in the way just described throughout each of the two parts of the evolution; that is, during the $H_{x/z}$ part of the evolution, the Lindblad operators will be taken to be of the form of Eq. 27 with $\ket{m}$ the eigenstates of the Hamiltonian $H_{x/z}$. This is valid so long as the bath equilibration timescales are much shorter than the times $T_{x/z}$ over which each part of the Hamiltonian is acting.

The particular form of the coupling operator we shall take is

[TABLE]

which commutes with the parity and therefore does not couple the two sectors at all.

As already anticipated, there is a fundamental difference between this form of Lindblad operators and the dephasing operators we looked at earlier: this choice does not lead to a fully-mixed state but rather to some periodic steady-state. Physically this is because this type of bath is not energetically structureless and does not cause jumps between all pairs of eigenstates with equal probabilities; it rather favours a form of the density matrix that is a) diagonal in the instantaneous energy eigenstates and b) has a Gibbs energy distribution, that is, the preferred density matrix is a statistical mixture of eigenstates with Gibbs weights. In the Floquet case, the unitary dynamics rather favours the BDE form of Sec. III. This results in a competition between the two forms, leading to nontrivial long-time steady states. We show this explicitly by numerically solving the Lindblad equation and plotting the expectation value of the instantaneous Hamiltonian in the right panel of Fig. 5. The main figure shows the time evolution of this energy starting from the same initial state discussed in the previous sections and up to some late times for two different values of the inverse temperature of the external bath, $\beta$. The inset shows the periodic steady-state reached by the system at long times. Evidently, for both values of $\beta$, the long-time state is not fully mixed, is not an infinite-temperature state, and has non-trivial dynamics, in contrast to the two dephasing cases presented earlier.

Having established the nontriviality of the long-time steady-state, in Fig. 6 we show that nevertheless the DTC order decays under the influence of this external parity-preserving bath. In both panels the two traces indicate different temperatures (and the same bath-system strength coupling); both cases results in a decay, leading to a broadening of the peak in the Fourier transform shown in the right panel.

VI Conclusions

In this work we have presented an analysis of the long-time periodic steady-state of Floquet systems with nontrivial spectral structure. These systems include the $\pi$-spin glass, a.k.a. discrete time crystal, in which all states appear as members of a doublet separated by a fixed quasienergy, leading to a subharmonic response at a doubled period $2T$. The resulting long-time density matrix has the form of a block diagonal ensemble (BDE) Eq. (15), rather than the diagonal ensemble of systems with no special spectral structure. Such a structure is necessary but not sufficient for the characteristic subharmonic oscillations that characterise DTCs; a further requirement is the presence of spin order in the eigenstates, which in one dimension requires the presence of disorder.

We have argued in general and supported with concrete numerical simulations that the subharmonic response is destroyed by any reasonable model of an environment, even if the symmetry underpinning the spectral doublets is not destroyed. This is a direct result of its reliance on disorder, for even static eigenstate order would be destroyed by the effect of broad-band environmental coupling. This holds true even if the external environment has a preferred energy scale (temperature) such that the system goes to a nontrivial long-time periodic steady state and preserves the symmetry, as for our “thermal” Lindblad operators in Sec. V.2, for which a steady state at finite energy density (thus, not infinite temperature) is obtained but nevertheless the DTC is destroyed.

This study leads to a number of interesting further topics. The first is a detailed analysis of the dependence of these phenomena on gross system parameters – in particular, studying the influence of system dimensionality is a natural follow-on from recent studies of many-body localisation in higher dimension.Chandran et al. (2016) This question is all the more pressing as one of the recent works studying the presence of time crystallinity was in fact undertaken on a three-dimensional system.Choi et al. (2016) This ties in naturally with a second major question, namely how one can hope to stabilise time-crystalline behaviour for a long “intermediate” time window, even if at asymptotically long times such order is completely lost. In particular, can one use the coupling to an environment to prolong, rather than curtail, the lifetime of subharmonic responses of the isolated system? There clearly remains much scope for exciting discoveries in both theory and experiment.

VII Acknowledgements

We thank G. J. Sreejith and M. Heyl for discussions as well as A. Das, V. Khemani and S. Sondhi for collaboration on earlier work.

Appendix A The Lindblad equation

We now explain how Markovian non-unitary evolution can be described for open quantum systems via a Lindblad master equation.

A.1 Setup

As mentioned in the text, the Lindblad equation (18) is the most general formulation of the differential time evolution of a system undergoing non-unitary Markovian dynamics; the operators $L$ describe the non-unitary part of the evolution and are completely arbitrary. Its formal solution is

[TABLE]

One intuitive interpretation of this evolution is that the system evolves unitarily, while at random intervals undergoing a measurement induced by the operator $L_{\alpha}$. The result of the measurement is discarded, leaving the system in a mixed state, and the evolution continues. As we are taking the non-unitary evolution to be due to the interaction with some macroscopically large external environment, this interpretation is intuitively appealing.

In general, the Lindblad equation has at least one fixed point for which $\partial_{t}\rho=0$; more generally, the fixed point is the eigenvector of $\mathcal{L}$ corresponding to the eigenvalue with 0 real part. As we shall see below, the presence of symmetries ensures the existence of multiple fixed points.

Finally, let us note that if all the $L_{a}$ are Hermitian, $L_{a}=L_{a}^{\dagger}$ then a fixed point is the identity $\rho=\mathbb{I}$. Thus for Hermitian Lindblad operators typically the steady state is the fully-mixed, infinite temperature state $\rho=\mathbb{I}$.

A.2 Symmetries

Analogously to the Heisenberg picture in unitary quantum mechanics, one may define the adjoint Lindblad operator by its action on some observable $A$ as follows:Ribeiro and Vieira (2015)

[TABLE]

The Lindblad operator and its adjoint are related by

[TABLE]

from which it follows that the time-dependent expectation value is

[TABLE]

We now concentrate on the case where the Hamiltonian commutes with some operator, taking for concreteness the parity operator, Eq. 5. Noticing that one may write the adjoint operator as

[TABLE]

we immediately conclude that if $\left[L_{a},P\right]$ for all $a$ then the expectation value of the operator $P$ is time independent; it therefore remains a conserved quantity.

As the operator $P$ has two possible values, this implies that there must be two fixed points.

A.3 Time-dependent systems

The Lindblad equation is first order in time and linear in $\rho$. This implies that its solution for a time-depenent $\mathcal{L}\left(t\right)$ can be obtained in a very similar way to that for the Schrodinger equation. In particular, if the Lindblad operators $L_{a}$ are time-independent and $H\left(t\right)$ as in Eq. 1 then the operator propagating $\rho$ over a period $T$ is

[TABLE]

with $\mathcal{L}_{\alpha}\rho=-i\left[H_{\alpha},\rho\right]+\sum_{a}\left(L_{a}\rho L_{a}^{\dagger}-\frac{1}{2}\left[L_{a}^{\dagger}L_{a},\rho\right]_{+}\right)$ for $\alpha=x,z$ so that time evolution is generated by $\rho\left(T\right)=\mathcal{K}\rho\left(0\right)$. This defines a quantum map. We note that $\mathcal{K}$ need not be the exponential of a time-dependent Lindblad operator.

Appendix B Analytical approach to dephasing in the $x$ direction

Because of the special feature of this model that decay occurs only during the rotation phase when $H_{z}$ is acting, and because $H_{z}$ is a sum of single-spin terms, we can write down a single-spin simplified version of the problem which reproduces the observed dependence of the rates on the parameters.

During the part of the period over which the ferromagnetic term $H_{x}$ is acting and the system is evolving according to $\mathcal{L}_{x}$, the spins align along the $x$ direction. While so aligned, the dephasing term has no effect as any state with all spins (anti)aligned with the $x$ axis is a steady state of $\mathcal{L}_{x}$ with the Lindblad operators of Eq. (24). On the other hand, the paramagnetic term $H_{z}$ rotates the spin at position $i$ by $h_{i}T_{z}$ which is on average $\pi$; thus the spins approximately flip. This flipping does not have to be exact, due to the action of the $H_{x}$. However during the process of flipping each spin is obviously not in a steady state of $\mathcal{L}_{z}$ so the dephasing term is effective.

We now abstract this process into a single-spin evolving under single-spin versions of $H_{z}$ and the Lindblad operators, taking

[TABLE]

and a dephasing Lindblad operator

[TABLE]

with the density matrix evolving as in Eq. 18. We will use this to propagate our spin forward in time by $\tau=h/\pi$, so that in the unitary case $\gamma=0$ the spin’s $x$ component would exactly flip.

The density matrix can be decomposed as

[TABLE]

and Eq. 18 results in

[TABLE]

The $\rho_{0}$ is constant, preserving the trace, while $\rho_{z}$ decays exponentially if it is not zero. The $x$ and $y$ components perform damped oscillations, as can be seen by writing $\omega_{0}=h$ and $\eta=\gamma/h$ whence

[TABLE]

which is a damped harmonic oscillator: for $\eta>1$ or $\gamma>h$ it is overdamped and the spin decays exponentially (without any oscillations) to $\rho_{x}=\rho_{y}=\rho_{z}=0$. On the other hand, for $\eta<1$ or $\gamma<h$ it is underdamped, decaying to the same point but oscillating on the way with frequency $\sqrt{h^{2}-\gamma^{2}}$ and decay constant $\gamma$.

Integrating this for time $\tau=\pi/h$, which in the absence of damping would result in the spin flipping, the condition for being underdamped is

[TABLE]

For the initial condition $\rho_{y}=\rho_{z}=0$ and $\rho_{x}\neq 0$ we obtain the result displayed in Eq. (25).