Fate of a discrete time crystal in an open system

Following the recent realization 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.

Figure 1

Magnitude of matrix element $\left.〈\omega \right|{\sigma }_1^x\left|{\omega }^{\prime }\right.〉$ for a single disorder realization and a system size $L=10$ of the system of Eq. (1) for bond- and site-disorder amplitudes $\delta J=w=0.1J$ vs the quasienergies of the Floquet states between which it is taken. Darker color indicates larger magnitude. The weight is mostly concentrated at elements between Floquet states differing by quasienergy $\pi /T$, implying that the dynamics of the operator ${\sigma }_1^x$ will display a strong oscillation at frequency $\pi /T$ or period $2T$ as in Eq. (11).

Figure 2

Left: Real-time plot of $〈{\sigma }_1^x〉\left(t\right)$ using the Lindblad dephasing operators of Eq. (24). The average $h_i$, denoted by $\overline{h}$, is chosen such that $\overline{h}{T}_z=\pi$, while each site has a different, random $h_i$. The blue (triangular) and red (circular) markers indicate the choice $T_x=T_z=1$ $\left(T_x=T_z=4\right)$, and the time axis is rescaled by $T=T_x+T_z$ to allow comparison of the two sets. The main feature is that as $T_z$ is the longer, the decay of the spatiotemporal order is faster, as explained in Sec. 5a. The data displayed are for size $L=6$. Note that the decay of the DTC occurs during the part of the period where the spins are rotating, while there is no decay during the part where the plateaus are in contrast to the results of Fig. 4. Right: Fourier transform of the time evolution, normalized so that the peaks have the same height to facilitate comparison of the width. The important features are that (a) both cases have a peak at half the driving frequency and (b) slower driving results in a broader frequency (corresponding to a faster-decaying oscillation in real time, as on the left panel). We note that the decay due to the finite-size splitting of the members of each eigenstate pair is not visible on these time scales.

Figure 3

Main panel: Decay rates vs $T_z=\pi /\overline{h}$, the time over which the $\pi$ rotation is effected, for sizes $L=4,5$ in blue, red respectively (the data points coincide so the blue are not visible). Here, $J=T_x=1,\delta J=0.2,\mathrm{\Gamma }=0.2$ and the data have been averaged over 100 (400) disorder realizations for $L=4\left(5\right)$. Inset: Same as left, but with $\mathrm{\Gamma }{T}_z$ kept constant (by varying $\mathrm{\Gamma })$. The decay rate is constant, indicating that the decay rate is a function of the product $\mathrm{\Gamma }{T}_z$ only. This is consistent with the results of B. The data have been averaged over 200 disorder realisations for both $L=4,5$.

Figure 4

Left: Real-time plot of $〈{\sigma }_1^x〉\left(t\right)$ using the Lindblad dephasing operators of Eq. (26), for times [see Eq. (1)] $T_x=T_z=1$. The blue (triangular) and red (circular) markers indicate different dephasing rates. The data displayed are for size $L=5$. The main qualitative difference from Fig. 2 is that the decay of the oscillations occurs also during the part of the period where the spins are aligned with the $x$ axis (indicated by the sloped plateaus in this figure and the flat plateaus in Fig. 2). Right: Fourier transform of the data in the left panel. The frequency $\mathrm{\Omega }$ is scaled by the driving frequency $\omega$ and the vertical axis is scaled so that the highest value of each trace is 1 to make the broadening easier to see. Note that, first, the peak is at $\omega /2$ or half the driving frequency, and second, the spectrum broadens for stronger dephasing.

Figure 5

Instantaneous expectation value of average Hamiltonian $〈H_0〉$ for the case of dephasing (left panel; see Sec. 5a2) and thermal (right panel; see Sec. 5b) Lindblad operators. The green lines indicate the average of all the energy eigenvalues, which is the expectation value of $H_0$ in the infinite-temperature Gibbs state (or the fully mixed state). In the dephasing case, the final state is the fully mixed, featureless state. In the thermal case, the steady-state is a time-periodic state with energy below that of the infinite temperature state. The blue and red traces correspond to bath inverse temperatures $\beta =0.5$ and 0.05, respectively. The inset shows the expectation value over two cycles of this steady state with a linear time axis (this is obtained by direct calculation at a very large time of $t/T\sim 4000$ but the results are the same for larger times: This is a periodic steady state. For dephasing operators in the $x$ direction, the results are very similar to those shown in the left panel.

Figure 6

Left: Real-time plot of $〈{\sigma }_1^x〉\left(t\right)$ using the thermal Lindblad operators of Sec. 5b, for times [see Eq. (1)] $T_x=T_z=1$. Such operators lead to thermal (Gibbs) steady states in the static case at some inverse temperature $\beta$. The blue (triangular) and red (circular) markers indicate different environment temperatures, and the system-bath interaction does not couple eigenstates of different parity. The data displayed are for size $L=5$ and the system-bath coupling strength $\gamma =0.1$. Right: Fourier transform for the same parameters, showing the broadening due to the finite-temperature bath. The main insight to be gained from this figure is that generic system-bath couplings, even at finite temperatures (see Fig. 5), destroy the DTC.