Critical Time Crystals in Dipolar Systems

TL;DR

This paper investigates the stability and lifetime of discrete time crystalline order in 3D dipolar systems, revealing a critical regime stabilized by slow dynamics rather than localization, with implications for quantum matter probing.

Contribution

It introduces a perturbative approach to evaluate DTC lifetime in dipolar systems and predicts a sharp crossover indicating a novel critical DTC regime.

Findings

Decay of DTC order is parametrically slow in 3D dipolar systems.

A sharp crossover from stable DTC to loss of order is predicted.

Experimental results align with the theoretical predictions of a critical DTC regime.

Abstract

We analyze the quantum dynamics of periodically driven, disordered systems in the presence of long-range interactions. Focusing on the stability of discrete time crystalline (DTC) order in such systems, we use a perturbative procedure to evaluate its lifetime. For 3D systems with dipolar interactions, we show that the corresponding decay is parametrically slow, implying that robust, long-lived DTC order can be obtained. We further predict a sharp crossover from the stable DTC regime into a regime where DTC order is lost, reminiscent of a phase transition. These results are in good agreement with the recent experiments utilizing a dense, dipolar spin ensemble in diamond [Nature 543, 221-225 (2017)]. They demonstrate the existence of a novel, critical DTC regime that is stabilized not by many-body localization but rather by slow, critical dynamics. Our analysis shows that the DTC response can be used as a sensitive probe of nonequilibrium quantum matter.

Equations

$$H_0=\sum _i\Omega S_i^x+\sum _{i,j}\frac{J_{ij}}{r_{ij}^3}{S}_i^x{S}_j^x,$$

$$U_F=\exp \left[-i\sum _i(\pi +ϵ)S_i^y\right]\exp \left[-i{H}_0{\tau }_1\right],$$

$$H(t)$$

$${\Delta }_{IJ}\approx |h_I+h_J-m^{\ast }{\omega }_0|<J_{IJ}^{\text{eff}},$$

$$J_{IJ}^{\text{eff}}\sim {\left(\frac{ϵ}{{\tau }_1}\right)}^2\frac{J_{IJ}}{r_{IJ}^3}\left(\frac{1}{{\delta }_I^2}+\frac{1}{{\delta }_J^2}\right),$$

$$\displaystyle C\approx\left\{\begin{array}[]{cc}(\epsilon/\tau_{1}W)^{2}&\textrm{for }\omega_{0}\gg W\\ \epsilon^{2}&\textrm{for }\omega_{0}\ll W\end{array}\right.,$$

$$dP=\left(J^{\text{eff}}(R)/W^{\text{eff}}\right)n_{0}4\pi R^{2}dR,$$

$$\displaystyle\Gamma\sim\left\{\begin{array}[]{cc}\frac{J_{0}\epsilon^{2}}{a_{0}^{3}\tau_{1}W^{2}}\exp{\left[-\frac{3W^{3}\tau_{1}^{2}}{4\pi J_{0}n_{0}\epsilon^{2}}\right]}&\textrm{for }\omega_{0}\gg W\\ \frac{J_{0}\epsilon^{2}\tau_{1}}{a_{0}^{3}}\exp{\left[-\frac{3}{2J_{0}n_{0}\epsilon^{2}\tau_{1}}\right]}&\textrm{for }\omega_{0}\ll W.\end{array}\right.$$

$$H^{\prime }(t)$$

$$\frac{V}{{\tau }_1}-[{\Omega }^{(n)},H_0]+n{\omega }_0{\Omega }^{(n)}=0.$$

$$⟨s^{\prime }|{\Omega }^{(n)}|s⟩=\frac{⟨s^{\prime }|V|s⟩}{(E_s-E_{s^{\prime }}-n{\omega }_0){\tau }_1}.$$

$$H^{\prime }(t)=H_0-\frac{{ϵ}^2}{2}[\Omega (t),V]\sum _n\delta (t-n^{-}{\tau }_1)$$

$$H^{\prime }(t)=H_0+\sum _{IJ}\frac{A_{IJ}{J}_{IJ}}{r_{IJ}^3}\left(S_I^{+}{S}_J^{+}+\text{H.c.}\right)\sum _n\delta (t-n^{-}{\tau }_1),$$

$$A_{IJ}\approx -2S_I^x(s)S_J^x(s){\left(\frac{ϵ}{{\tau }_1}\right)}^2\left(\frac{1}{{\tilde{\delta }}_I^2}+\frac{1}{{\tilde{\delta }}_J^2}\right),$$

$$H(t)=\sum _{ij}\frac{J_{ij}}{r_{ij}^3}{S}_i^x{S}_j^x+ϵ\sum _i{S}_i^y\sum _n\delta (t-n^{-}{\tau }_1),$$

$$U(t)=Q(t)\tilde{U}(t)Q^{†}(0),$$

$$\tilde{U}(t)=\mathcal{T}{e}^{-i{\int }_0^t{H}^{\prime }(t)dt},$$

$$H^{\prime }(t)$$

$$H_0=\sum _{ij}\frac{J_{ij}}{r_{ij}^3}{S}_i^x{S}_j^x,$$

$$V=\sum _i{S}_i^y.$$

$$⟨\psi |S_i^x(n{\tau }_1)|\psi ⟩=⟨\psi |U^{†}(n{\tau }_1)S_i^{x}U(n{\tau }_1)|\psi ⟩$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}=⟨\tilde{\psi }|{\tilde{U}}^{†}(n{\tau }_1){\tilde{S}}_i^x\tilde{U}(n{\tau }_1)|\tilde{\psi }⟩,$$

$$|\tilde{\psi }⟩=Q^{†}(0)|\psi ⟩$$

$${\tilde{S}}_i^x=Q^{†}(0)S_i^{x}Q(0),$$

$$⟨\psi |S_i^x(n{\tau }_1)|\psi ⟩\approx ⟨\psi |{\tilde{U}}^{†}(n{\tau }_1)S_i^x\tilde{U}(n{\tau }_1)|\psi ⟩.$$

$$H^{\prime }(t)$$

$$+{ϵ}^2\left(\frac{1}{2}[\Omega ,[\Omega ,H_0]]-[\Omega ,V]\sum _n\delta (t-n^{-}{\tau }_1)+\frac{i}{2}[\Omega ,{\partial }_t\Omega ]\right)$$

$$\cdots .$$

$$\frac{V}{{\tau }_1}-[{\Omega }^{(n)},H_0]+{\omega }_{0}n{\Omega }^{(n)}=0.$$

$$⟨s^{\prime }|{\Omega }^{(n)}|s⟩=\frac{⟨s^{\prime }|V|s⟩}{(E_s-E_{s^{\prime }}-n{\omega }_0){\tau }_1}.$$

$$E_s-E_{s^{\prime }}=2\sum _{j\ne I}\frac{J_{jI}}{r_{jI}^3}{S}_j^x(s)S_I^x(s),$$

$$\frac{ϵ}{{\tau }_1}\ll \underset{m\in \mathbb{Z}}{\min }\left|\left(2\sum _{i\ne J}\frac{J_{iJ}}{r_{iJ}^3}{S}_i^x(s)S_J^x(s)-m{\omega }_0\right)\right|.$$

$$ϵ<\min (W,{\omega }_0){\tau }_1,$$

Full text

Introduction.—Understanding and controlling non-equilibrium quantum matter is an exciting frontier of physical science. While periodic driving has long been used to control the properties of quantum systems, it was only recently realized that periodically driven (Floquet) systems can also host new states of matter that are not possible in equilibrium. In particular, this is possible in a class of driven disordered systems exhibiting many-body localization (MBL) Nandkishore and Huse (2015), called Floquet-MBL systems, which avoid unbounded heating to infinite temperature Ponte et al. (2015a); Lazarides et al. (2015); Abanin et al. (2016). The latter is generally expected to befall all ergodic isolated systems due to external driving Lazarides et al. (2014); Ponte et al. (2015b); D’Alessio and Rigol (2014).

One remarkable example of a novel, nonequilibrium phase of matter is the recently introduced discrete time crystal (DTC) Khemani et al. (2016); Else et al. (2016); von Keyserlingk et al. (2016); Yao et al. (2017), which is characterized by a spontaneously broken discrete time-translation symmetry of the underlying drive. In such systems, physical observables exhibit robust oscillations with a period that is an integer multiple of the underlying driving period $T$. Indeed, key signatures of such robust DTC order have been observed in two recent experiments Zhang et al. (2017); Choi et al. (2017). In particular, one of these realizations involved a disordered ensemble of $\sim 10^{6}$ spins associated with nitrogen-vacancy (NV) centers in diamond, which interact between themselves via dipolar couplings Choi et al. (2017). The origin of apparent robustness of the observed DTC order in such a system Choi et al. (2017), however, has not been fully understood. Although this system is disordered due to the random positions of the NV centers in 3D, the long-range dipolar interactions are believed to preclude localization Anderson (1958a); Levitov (1990); Burin (2006); Yao et al. (2014); Gutman et al. (2016). Moreover, a prethermal regime of the DTC Else et al. (2017), was also ruled out Choi et al. (2017) since in the experiment the initial polarized state is effectively at an infinite temperature with respect to the effective Hamiltonian due to the randomly varying signs of the dipolar interactions. Since neither localization nor prethermalization are likely the mechanisms that stabilize the DTC order, this raises important questions about the origin of the observed robust DTC response.

This Letter develops a theoretical treatment of DTC order in systems with long-range interactions. We utilize a perturbative approach to analyze the interplay of long-range interactions, periodic driving, and positional disorder of spins. Focusing on dipolar systems in 3D, we show that although DTC order is only transient, it can persist for asymptotically long times with a strongly suppressed thermalization rate. This behavior is intrinsically connected to the slow thermalization dynamics of disordered dipolar systems in 3D, which has been previously shown to be consistent with the so-called critical regime Anderson (1958a); Kucsko et al. (2016) without a periodic drive. As a function of experimental parameters, we find that the relaxation time shows a sharp crossover between a regime where the DTC response is robust and a regime where it decays rapidly. This crossover is reminiscent of a phase transition, thereby allowing us to obtain the effective phase diagram of the DTC which is in good agreement with experimental results. Thus, our work provides an explanation of the recent experimental observations Choi et al. (2017), and also demonstrates the possibility of the DTC in systems with critical dynamics, a regime which we refer to as ‘critical time crystals’. Furthermore, our perturbative approach can be used to study the nonequilibrium properties in other driven disordered systems with long-range interactions.

Our key results can be understood by considering a simple spin model that describes an ensemble of dipolar interacting NV centers, used in the experiments of Ref. Choi et al. (2017). Using strong microwave excitations, the effective Ising interactions between spins were engineered, described by the following Hamiltonian:

[TABLE]

where $\vec{S}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ are Pauli spin-1/2 operators, $\Omega$ the strong microwave driving along $\hat{x}$, $J_{ij}$ the orientation-dependent coefficient of dipolar interactions with typical strength $J_{0}$, and $r_{ij}$ the distance between spins $i$ and $j$. We assume that the spin-1/2 particles are randomly distributed in three-dimensional space with density $n_{0}$ and neglect coupling to the environment [Fig. 1(a)]. DTC order was observed by interrupting the evolution under Hamiltonian (1) with rapid, global pulses that rotate the spin ensemble along the $\hat{y}$ axis by an angle $\pi+\epsilon$. The corresponding Floquet unitary is given by

[TABLE]

where $\tau_{1}$ is the period for which the spins are allowed to interact for. In the experiment, the period is chosen such that $\Omega\tau_{1}=2\pi n$, and therefore $\Omega$ can effectively be taken to be [math] in Eq. (1). When the system is initialized in a state where all spins are polarized along the $+\hat{x}$ direction, a nontrivial temporal response may be revealed by measuring the average polarization $P(n\tau_{1})$ of the ensemble along $\hat{x}$ after $n$ Floquet cycles, or equivalently, $q(n)\equiv(-1)^{n}P(n\tau_{1})$, which serves as an order parameter for the DTC phase. The stability of the DTC order can be ascertained by studying the decay rate of $q(n)$ for large number of cycles as a function of $\tau_{1}$ and $\epsilon$.

In order to describe the dynamics of $q(n)$, we move into a so-called toggling frame, which rotates by $P_{\pi}\equiv\prod_{j}\exp{\left[-i\pi S_{j}^{y}\right]}$ each time a global pulse is applied to the system. Since $P_{\pi}S_{i}^{x}(P_{\pi})^{-1}=-S_{i}^{x}$, the $2\tau_{1}$-periodic oscillation in $P(n\tau_{1})$ naturally appears as a time-independent spin polarization in this new frame. The dynamics of the system is then described by the Floquet unitary $\bar{U}_{F}=\exp{\left[-i\sum_{i}\epsilon S_{i}^{y}\right]}\exp{\left[-iH_{0}\tau_{1}\right]}$, or, equivalently, by an effective time-dependent Hamiltonian

[TABLE]

Thus, our problem reduces to studying the depolarization dynamics of an initialized polarized spin ensemble under the time evolution of $H(t)$.

Physical picture.—The essence of our analysis is to study resonant spin dynamics that lead to depolarization perturbatively in $\epsilon$, while accounting for energy exchanges provided by the external drive. In particular, since $\sum_{n}\delta(t-n^{-}\tau_{1})=\frac{1}{\tau_{1}}\sum_{m}e^{im\omega_{0}t}$, the pulsed periodic spin rotations can be viewed as spin excitation with harmonics of the fundamental frequency $\omega_{0}\equiv 2\pi/\tau_{1}$ and fixed magnitude $\epsilon/\tau_{1}$. While this driving allows energy absorption and emission in integer multiples of $\omega_{0}$, the interplay of strong interactions and positional disorder suppresses direct energy exchanges such that typical spins depolarize only via indirect higher-order processes in $\epsilon$.

Let us first consider the case without perturbations, i.e. $\epsilon=0$. Then the polarization of each spin along $\hat{x}$ is conserved. When all spins are initially polarized, each spin therefore experiences a mean-field potential $h_{i}\equiv\sum_{j\neq i}(J_{ij}/r_{ij}^{3})\langle S_{j}^{x}\rangle$. Because of the random positioning of spins, the strength of $h_{i}$ is also random with zero mean and variance $W^{2}=\langle\frac{1}{4}(\sum_{j\neq i}J_{ij}/r_{ij}^{3})^{2}\rangle$, where $\langle\cdot\rangle$ denotes averaging over different positions.

When $\epsilon\neq 0$, there is depolarization due to spins experiencing a time-varying on-site field along the $\hat{y}$ axis. Let us therefore consider the first-order process where spins individually flip due to the action of this field. If a spin experiences a a strong mean-field potential $h_{i}$ compared to the applied field, that is, if $h_{i}\tau_{1}\gg\epsilon$, then it does not flip – it experiences an effective field that is approximately pointing along the $\hat{x}$ axis and therefore precesses around it without significant depolarization. On the other hand, if $h_{i}\tau_{1}$ is close to an integer multiple of $2\pi$, then the spin rotates along the $\hat{y}$ axis and depolarizes. Physically, this corresponds to an effectively resonant excitation of (individual) spins that arises when one of the driving harmonics is tuned close to their energy: $|h_{i}-m^{*}\omega_{0}|<\epsilon/\tau_{1}$ for some optimal integer $m^{*}$ [see Fig. 1(b)]. Such resonances occur with a small probability in the limit of $\epsilon\ll W\tau_{1}$, and amount to a reduction of the total polarization by a constant factor proportional to $\epsilon/({\rm min}(W,\omega_{0})\tau_{1})$. However, if $\epsilon\sim W\tau_{1}$, a substantial fraction of spins rapidly depolarize due to resonant processes shown in Fig. 1(b). Note that the phenomenological phase boundary extracted in Ref. Choi et al. (2017), based on the existence of self-consistent closed spin trajectories, is consistent with the perturbative condition $\epsilon\ll W\tau_{1}$.

We next focus on the second-order process illustrated in Fig. 1(c) and (d) in which a pair of spins $I$ and $J$ simultaneously flip their polarizations while exchanging their energies with each other and with the external drive. Such processes are resonant when

[TABLE]

where $h_{I}$ and $h_{J}$ are effective on-site potential for spins $I$ and $J$, respectively, $m^{*}$ is the optimal harmonic number that minimizes the energy difference, and $J^{\textrm{eff}}_{IJ}$ is the effective amplitude of the pair-flip process. The amplitude $J^{\textrm{eff}}_{IJ}$

[TABLE]

can be estimated from the interference of two paths in the second-order perturbation theory, as illustrated in Figs. 1(c) and 1(d). Here $\delta_{I(J)}\equiv{\rm min}_{\ell^{*}}(h_{I(J)}-\ell^{*}\omega_{0})$ is the energy difference between initial or final states and intermediate virtual states, up to extra energy provided by a driving harmonic $\ell^{*}$. We find that $J^{\textrm{eff}}_{IJ}$ is an effective long-range interaction decaying as $\sim 1/r_{IJ}^{3}$ allowing the flipping of remote spin pairs.

The resonance condition (4) is sensitive to $\omega_{0}$ and behaves qualitatively differently in two limiting cases: (i) $\omega_{0}\gg W$ and (ii) $\omega_{0}\ll W$ [see Fig. 1(c,d)]. In the former case, the optimal choice is $m^{*}=\ell^{*}=0$ since spins cannot absorb or emit such a large energy $\omega_{0}$. In the latter case, effective energy differences (both $\Delta_{IJ}$ and $\delta_{I},\delta_{J}$) are bounded by $\omega_{0}$ as the external drive can always compensate energy in units of $\omega_{0}$. These considerations yield the scaling $J^{\textrm{eff}}(r)\sim CJ_{0}/r^{3}$ with

[TABLE]

and the effective range $W^{\textrm{eff}}$ of the energy differences $\Delta_{IJ}$ becomes $W^{\textrm{eff}}\sim W$ for $\omega_{0}\gg W$ and $W^{\textrm{eff}}\sim\omega_{0}$ for $\omega_{0}\ll W$.

We now estimate the probability that a given spin finds a resonant partner within a ball of radius $R$. This is obtained by integrating the probability of finding such a partner in a shell $R$ and $R+dR$

[TABLE]

from a short distance cutoff $a_{0}$ to $R$, which gives $P(R)\sim\log(R/a_{0})$. Here the first factor in Eq. (9) is the probability of satisfying Eq. (4), and the second factor is the average number of spins within a shell of size $R$ with the density $n_{0}$. As this probability diverges, it implies that pairwise spin flips prevail, and the system thermalizes, with the DTC order slowly decaying over time. We can extract the time scale associated with these pair-spin-flip processes using the typical distances $R_{*}$ of resonant spin pairs. Solving $P(R_{*})\sim 1$ gives $R_{*}\approx a\exp{[W^{\textrm{eff}}/4\pi CJ_{0}n_{0}]}$. Finally, the effective depolarization rate is estimated from the interaction strengths of typical pairs, i.e., $\tilde{\Gamma}\sim J^{\textrm{eff}}(R_{*})$, leading to the decay rate per Floquet cycle $\Gamma\equiv\tilde{\Gamma}\tau_{1}$:

[TABLE]

This exponentially slow in $1/\epsilon^{2}$ decay of the DTC order is a central result of the present Letter and is a direct consequence of critically slow thermalization of dipolar systems in 3D Anderson (1958b); Kucsko et al. (2016). Interestingly, the depolarization is exponentially sensitive to the parameters $\tau_{1}$ and $\epsilon$ in two distinct ways: In regime (i) $\Gamma$ is a function of $\tau_{1}^{2}/\epsilon^{2}$ while in regime (ii) it only depends on $1/\epsilon^{2}\tau_{1}$. These considerations allow us to identify an effective phase boundary using the criteria $\tau_{1}^{2}/\epsilon^{2}=A$ or $1/\epsilon^{2}\tau_{1}=B$ with some constants $A$ and $B$. Remarkably, this boundary illustrated in Fig. 2 captures the key features observed in the experiment Choi et al. (2017): the linear growth of $\epsilon$ for short $\tau_{1}$ and slow diminishing of $\epsilon$ at longer $\tau_{1}$ Cor .

Technical procedure.—We now outline the technical procedure that formalizes the above discussion (see SOM for details). The key idea is to identify a time-dependent unitary transformation of the Hamiltonian $H(t)$ such that nonresonant single spin flips are essentially “integrated out” and only residual two-spin-flip processes become dominant terms in the effective Hamiltonian $H^{\prime}(t)$. More specifically, we start from the Hamiltonian (3) with $H_{0}$ representing the Ising interactions and $V$ the applied field, and perform a time-periodic unitary transformation $Q(t+\tau_{1})=Q(t)$, which gives rise to

[TABLE]

Our goal is to eliminate terms that are linear in $\epsilon$ from $H^{\prime}(t)$. Following Ref. Abanin et al. (2016), we look for $Q(t)$ of the form $Q(t)=e^{\epsilon\Omega(t)}$ with anti-Hermitian operator $\Omega(t)=\sum_{n}\Omega^{(n)}e^{in\omega_{0}t}$. Expanding Eq. (22) in powers of $\epsilon$, and requiring that the $O(\epsilon)$ term equals 0 gives an equation for the $n$th Fourier mode $\Omega^{(n)}$:

[TABLE]

The matrix elements of the operator $\hat{\Omega}^{(n)}$ can be computed in the eigenstate basis $|s\rangle$ of $H_{0}$ (which is a product state basis in $S_{i}^{x}$ operators):

[TABLE]

Noting that ${V}=\sum_{i}S^{y}_{i}$, the operator $\Omega^{(n)}$ has nonzero matrix elements only between spin configurations $s$ and $s^{\prime}$ that differ by one spin flip. If $|s\rangle$ and $|s^{\prime}\rangle$ differ by the value of spin $I$, $E_{s}-E_{s^{\prime}}=2\sum_{j\neq I}\frac{J_{jI}}{r_{jI}^{3}}S_{j}^{x}(s)S_{I}^{x}(s)=2h_{I}S_{I}^{x}.$ We assume that the on-site field $h_{I}$ is random (due to positional disorder and orientation dependence of $J_{ij}$) and sufficiently strong such that resonances are rare, i.e. the denominator in (15) typically does not diverge and the procedure controlled. Then, the rotated Hamiltonian to second order becomes

[TABLE]

A straightforward calculation SOM ; fli using expression (15) gives an effective Hamiltonian of the following form:

[TABLE]

where $S_{I}^{+}\equiv(S_{I}^{z}+iS_{I}^{y})/\sqrt{2}$ is the spin raising operator in $S_{I}^{x}$ basis for the spin $I$, and $A_{IJ}$ is the coefficient

[TABLE]

where we introduced the notation $\frac{1}{\tilde{\delta}_{J}^{2}}=\sum_{\ell}\frac{1}{\left(h_{J}-\ell\omega_{0}\right)^{2}}.$

The effective Hamiltonian (17) contains the larger disordered part $H_{0}$, and long-range terms which can flip pairs of spins; the latter are suppressed proportional to $\epsilon^{2}$, leading to slow relaxation. From Eq. (18) it is evident that the amplitudes for flipping a pair of spins depend on $h_{I},h_{J}$, which in turn are determined by the positions of the spins. Assuming that $h_{I},h_{J}$ take typical values of the order $W$, and taking the contribution of the harmonic $\ell^{*}$ for which $h_{J}-\ell\omega_{0}$ is minimized (this gives the leading contribution to $\tilde{\delta}_{J}$), the expression (18) for the two-spin-flip amplitude reduces to the estimate (5) above.

We emphasize that the above unitary transformation is distinct from the rotating frame transformations employed to derive effective Hamiltonians in the high-frequency limit Abanin et al. (2017a, b). Rather, it utilizes the randomness of our Hamiltonian in order to effectively integrate out non-resonant single-spin-flip processes.

Phase diagram.—Using the effective Hamiltonian approach described above, we obtain the phase diagram of the critical DTC. To improve upon the estimates for $\Gamma(\epsilon,\tau_{1})$, we take into account the fact that the distribution of the potential $h_{i}$ stems from the positional randomness of spins, and numerically sample $h_{i}$ from a distribution of $2000$ spins in a 3D region with density $9.26\times 10^{-3}\text{\leavevmode\nobreak\ nm}^{-3}$ with a short distance cutoff of $3$ nm Kucsko et al. (2016).

While Eq. (12) already provides analytical predictions for the decay rate $\Gamma$ by estimating the typical distance $R_{*}$ of resonant spin pairs, in numerics we find it more amenable to estimate $\Gamma$ from an explicit depolarization in time profile; the counting arguments in Eq. (9) predicts a power-law decay of polarization $q(n)$, from which the decay time scale $1/\Gamma$ is extracted by equating $q(n)$ to a small threshold SOM . The phase boundary is then identified from a criterion $\Gamma(\epsilon,\tau_{1})=\Gamma_{*}=1/100$.

This approach yields the phase diagram illustrated in Fig. 2, which is in very good agreement with the experimental observations Choi et al. (2017); Per . At high driving frequency, the boundary approximately follows a relation $\tau_{1}\propto|\epsilon|$ (also obtainable using a semi-classical argument), while at low frequency $\tau_{1}\propto 1/\epsilon^{2}$, which indicates that DTC order becomes less stable as $\tau_{1}$ is increased, due to the fact that multiphoton processes lead to faster depolarization. The DTC phase is most robust in the crossover regime, where $\omega_{0}\sim W$.

We also note that, strictly speaking, DTC order has finite relaxation rate at any $\epsilon\neq 0,\tau_{1}\neq 0$. However, we find that the relaxation rate $\Gamma$ increases very sharply at a certain value of $\epsilon$, as illustrated in Fig. 3, which matches the experimental observations and is reminiscent of a phase transition. Note, however, that unlike for a true phase transition, this increase does not become infinitely sharp even in the thermodynamic limit.

Summary and discussion.—We described a new approach to analyze the dynamics of periodically driven spin systems with long-ranged interactions and applied it to explain the recently observed surprising stability of DTC in dipolar spin system. The results of our analysis are in very good agreement with experimental observations. They demonstrate that these observations correspond to a novel, critical regime of the DTC order.

Furthermore, our general approach can be applied to analyze the interplay of long-range interactions, randomness, and periodic driving in a broad class of experimental systems. The present analysis focused on the experimentally relevant case of critical interactions, decaying as $1/r^{\alpha}$, where $\alpha$ coincides with the dimensionality of the system, $\alpha=d=3$. This leads to direct relaxation processes of spin pairs. It is interesting to extend the analysis to the case $\alpha>d$ (e.g. $\alpha=3,d=2$), where resonant spin-pair-flip processes are rare and presumably do not provide the main relaxation channel. Experimentally, such a situation can be realized by reducing the dimensionality of the dipolar spins systems. In the static case, relaxation is expected to occur via multispin processes: in essence, a sparse resonant network may form, which can act as a heat bath that mediates the relaxation of other spins Burin (2006); Gutman et al. (2016). We expect that future experiments on DTC in reduced dimensions will allow one to probe such a delicate interplay of various relaxation mechanisms in driven systems with long-range interactions. Our theoretical approach is well suited for analyzing such systems. Finally, apart from these specific realizations, our analysis demonstrates that the DTC response to periodic perturbations can be used as a sensitive probe of nonequilibrium quantum states and phases of matter.

I Supplemental material: Critical time crystals in dipolar systems

II A. Details on technical procedure

In this section, we present the details on the technical procedure used to rotate the depolarization inducing Hamiltonian

[TABLE]

into an effective Hamiltonian $H^{\prime}(t)$.

As mentioned in the main text, because of the disorder in the interactions, single spin-flip processes, effected by the action of a single operator $S^{y}_{i}$ (i.e. the off-diagonal $O(\epsilon)$ term in Eq. (19)), are typically not resonant and do not induce significant depolarization. One has to consider other channels for depolarization which are of higher order in $\epsilon$, such as two, three, spin-flip processes, and inquire if they are resonant processes – if so, then the dominant channel of decay is the one that governs the asymptotic behavior of depolarization dynamics.

The purpose of the transformation we employ is thus to extract the terms that give rise to dominant depolarization processes: we will rotate the original Hamiltonian in such a way that off-diagonal terms that generate non-resonant processes are ‘integrated out’, giving a resulting effective Hamiltonian $H^{\prime}(t)$ whose leading order off-diagonal terms generate resonant proceses. In our system, it will turn out that the dominant decay channel is given by two spin-flip processes.

II.1 Rotating the Hamiltonian

To that end, let us perform this transformation in detail. We write the unitary time evolution operator $U(t)=\mathcal{T}e^{-i\int_{0}^{t}H(t^{\prime})dt^{\prime}}$, generated by Eq. (19), as

[TABLE]

where we have yet to define the unitary $Q(t)$. With this decomposition, $\tilde{U}(t)$ is given by

[TABLE]

i.e. it is generated by a rotated time-dependent Hamiltonian (via the Schrodinger equation):

[TABLE]

where

[TABLE]

Here $J_{ij}=J_{0}q_{ij}$ where $q_{ij}=-(1-3(\hat{z}\cdot\hat{r}_{ij})^{2})$ encodes the angular dependence of the interactions between spins $(i,j)$, as in Ref. Choi et al. (2017).

We will pick $Q(t)$ to be time-periodic; then because $Q(n\tau_{1})=Q(0)$, the expected value of observables as a function of time, such as the polarization (of one site) $S_{i}^{x}(n\tau_{1})$, is given by

[TABLE]

where $|\psi\rangle$ is the initial state which we take to be polarized in the $x$-direction. In other words, $Q(0)$ is just some static rotation that rotates both the state and observable.

If in addition $Q(0)$ is a ‘small’ rotation (as we will choose, and to be made precise below), then both the state and observable are close to the unrotated ones,$|\tilde{\psi}\rangle\approx|\psi\rangle$ and $\tilde{S}^{x}_{i}\approx S^{x}_{i}$, and one can conclude that

[TABLE]

That is, the time dependence (and consequently, depolarization) is completely captured in $\tilde{U}(t)$ and hence, $H^{\prime}(t)$, the rotated effective Hamiltonian.

II.2 Choosing the rotation $Q(t)$

Now let us construct $Q(t)$, which we write as $Q(t)=e^{\epsilon\Omega(t)}$, where $\epsilon$ is explicitly the small parameter and $\Omega(t)$ a periodic anti-Hermitian operator. Expanding Eq. (22) we get

[TABLE]

Utilizing a procedure similar to Ref. Abanin et al. (2016), we equate the order $\epsilon$ piece to [math] with the constraint that $\Omega(t)$ is time-periodic. Note that this transformation is distinct from the transformations employed in Refs. Abanin et al. (2017a, b) to generate effective Hamiltonians in high-frequency driven systems. There, the small parameter was the inverse of the driving frequency $\omega$, but here, in anticipation that we will take into account the disorder in the interactions, the small parameter is served by $\epsilon$, the strength of the off-diagonal perturbation to $H_{0}$. Decomposing $\Omega(t)$ in terms of its Fourier modes $\Omega(t)=\sum_{n}\Omega^{(n)}e^{in\omega_{0}t}$, where $\omega_{0}=2\pi/\tau_{1}$, and using the fact the Fourier transform of the Dirac comb $\sum_{n}\delta(t-n^{-}\tau_{1})$ is $\frac{1}{\tau_{1}}\sum_{n}e^{-in\omega_{0}t}$, the equation that the $n$-th Fourier mode has to obey is

[TABLE]

In the basis of the eigenstates $|s\rangle$ of $H_{0}$ which are product states in the $x$-direction, we can therefore write the solution as

[TABLE]

Note that since $V=\sum_{i}S^{y}_{i}$, the eigenstates $s$ and $s^{\prime}$ representing spin configurations in the matrix element of $\Omega^{(n)}$ can only differ by one spin-flip. If $|s\rangle$ and $|s^{\prime}\rangle$ differ at the $I$th spin,

[TABLE]

where $S_{j}^{x}(s)$ represents the $S^{x}$ component of the $j$th spin for the configuration $|s\rangle$, which for the starting state is just $+\frac{1}{2}$.

Because of the disorder in $J_{ij}/r_{ij}^{3}$ (accorded for by the random positions and relative angles between the spins), resonances (i.e. terms where the denominator $\approx 0$) are controlled, provided

[TABLE]

This gives us a condition that our perturbative procedure should work only for

[TABLE]

where $W^{2}$ is the variance of the interactions,

[TABLE]

The angular brackets represent averaging over different spins $i$, and we have used the fact that the mean is [math]. In that case, the corresponding fraction of spins that are resonant is then small, and goes as $\sim\frac{\epsilon}{\min(W,\omega)\tau_{1}}$.

II.3 Effective Hamiltonian

The rotated Hamiltonian (27) then becomes

[TABLE]

Let us concentrate on the second term and look at its matrix elements. This is

[TABLE]

Now $\Omega^{(m)}$ and $V$ are sums of terms which each individually flip a single spin, so $|s\rangle$ and $|s^{\prime}\rangle$ can differ by either only zero spin flips (i.e. $|s^{\prime}\rangle=|s\rangle$) or two spin flips ($|s^{\prime}\rangle=|s_{IJ}\rangle$) where $s_{IJ}$ stands for the spin configuration $s$ with the $I$-th and $J$-th spins flipped. The diagonal process (zero spin flips) leads to a renormalization of the energy of $H_{0}$, while the off-diagonal process where $V$ flips spin $I$ and $\Omega^{(m)}$ flips spin $J$ results in the matrix element

[TABLE]

The term in the parenthesis is

[TABLE]

where $\omega_{0}=2\pi/\tau_{1}$.

Using this, we can thus write the rotated Hamiltonian (34) in operator form:

[TABLE]

where $S_{I}^{+}\equiv S_{I}^{y}+iS_{I}^{z}$ is the spin-raising operator in the $S_{I}^{x}$ basis for spin $I$, $A_{IJ}$ is the coefficient of the interactions between two spins $(I,J)$, given by

[TABLE]

This effective Hamiltonian gives Eqs. (5, 13) and (14) in the main text. Note that technically speaking, Eqs. (37) when written in operator form would also give flip-flop terms $S^{+}_{I}S^{-}_{J}+h.c.$. However, because we are ultimately interested in the dynamics of depolarization beginning from an initial polarized state in the $\hat{x}$ direction, and moreover such flip-flop terms are polarization conserving, they do not contribute to depolarization and so we drop them in the effective Hamiltonan.

Depending on the two limits of high or low frequencies ($\omega_{0}\gg W$ or $\omega_{0}\ll W$ respectively), $A_{IJ}$ has qualitatively different behaviors. Its scaling behavior with respect to $\epsilon$ and $\tau_{1}$ can be extracted in those limits, which is Eq. (6) in the main text. Thus, this gives rise to different behaviors for the depolarization rate in the two limiting cases, and can ultimately be used to determine the phase boundary of the critical time crystal, as was done in the main text.

III B. Numerical procedure to extract decay rate and phase boundary

In this section, we describe the numerical procedure used to 1) determine the decay rate $\Gamma$ as a function of the perturbation $\epsilon$, and 2) obtain the phase diagram of the critical DTC shown in the main text.

Assuming that the initial state is polarized in the $x$-direction, we first generated the probability distribution $P(h_{I})$ of $h_{I}$, which is the mean-field potential field felt by spin $I$, i.e. $h_{I}=\frac{1}{2}\sum_{j\neq I}\frac{J_{Ij}}{r_{Ij}^{3}}$. Here $J_{ij}=J_{0}q_{ij}$ and $q_{ij}=-(1-3(\hat{z}\cdot\hat{r}_{ij})^{2})$ encodes the angular dependence of two spins $(i,j)$. We modeled this by uniformly distributing $N=2000$ spins in a 3-torus of linear dimension $L=60$ nm, with a minimum (UV) distance cutoff $r_{\text{UV}}=3$ nm, and caculating for each spin $h_{i}$. Such a choice of parameters gives a particle density of $n_{0}=9.26\times 10^{-3}$ nm*-3*. Using $J_{0}=2\pi\times 52$ Mhz nm*-3* as in the experiments, we can characterize this distribution by its variance $W^{2}$, which gives $W=4.49$ Mhz.

From this distribution, we then drew two values $(h_{I},h_{J})$, provided that they are each not resonant. That is, we discarded values if

[TABLE]

where $\alpha_{1}$ is some fixed $O(1)$ number. We determined this number from the condition that along the line $W\tau_{1}=|\epsilon|$ where our perturbative analysis holds, at least $80\%$ of the spins are not resonant. In the region where $W\tau_{1}>\epsilon$, this implies that $>80\%$ of the spins are off-resonant. In our simulation, $\alpha_{1}=3.7$. Then, $(h_{I},h_{J}$) represent the mean-field potentials felt by spins $I$ and $J$ respectively.

Next, we estimated the probability that the two spins $(I,J)$ a distance $r$ apart form a resonant pair. In principle this is obtained by comparing the effective hopping $J_{\rm eff}(r)\equiv A_{IJ}J_{IJ}/r^{3}$, against the minimum quasi-energy gap $\Delta$. Here $A_{IJ}J_{IJ}\sim-\left(\frac{\epsilon^{2}}{\tau_{1}^{2}}\right)J_{0}\overline{q_{IJ}}\sum_{m}\left(\frac{1}{(h_{I}-m\omega_{0})^{2}}+\frac{1}{(h_{J}-m\omega_{0})^{2}}\right)$ and $\Delta\approx\min_{m}\left|h_{I}+h_{J}-m\omega_{0}\right|.$ The bar over the angular dependence $q_{ij}$ represents the typical angular dependence; this can be derived analytically and gives $\overline{q_{ij}}=2/\sqrt{5}$. In other words, for a fixed $r$, if

[TABLE]

where $\alpha_{2}$ is another order one number, then this spin pair is resonant. We take $\alpha_{2}=3.4$.

In practice, we make a simplifying assumption that for large enough $r$, which we call $d$, the probability of finding two spins a distance $r$ apart which are resonant is small and is simply proportional to ${J_{\rm eff}(r)}/{\alpha_{2}\Delta}\propto 1/r^{3}$ (we have numerically checked this assumption holds).

Then, one can simplify the counting of resonant pairs by just estimating the the probability $\mathcal{P}$ of finding two spins some fixed distance $d$ apart. We generate $4000$ pairs $(I,J)$, and count the number of pairs that satisfy Eqn. (41) for $r=d$; the fraction of such pairs is $\mathcal{P}$. To get the probability of resonant pairs at distances $R>d$, we can then just multiply $\mathcal{P}$ by the factor $(d/r)^{3}$. In the simulations, we obtain $d$ as the distance for which the last step of Eqn. 37 is justified, that is, if $W\sim J_{0}\overline{q_{ij}}/d^{3}$, which gives $d=4.02$ nm for the parameters we have used.

Having determined $\mathcal{P}$, we extracted the decay rate $\Gamma(\epsilon,\tau_{1})$. The survival of polarization probability is given by a power law $P(t)=(t/t_{0})^{-q}$ where $q\propto\mathcal{P}$. This can be derived as the product of probabilities of having no resonant spins at each distance $[r,r+dr]$ up to $R(t)=(J_{0}t)^{1/3}$ Choi et al. (2017):

[TABLE]

where

[TABLE]

The factor in the parenthesis $\left(\frac{\overline{A_{IJ}J_{IJ}}}{\alpha_{2}\Delta d^{3}}\right)$ is the probability $\mathcal{P}$, and hence we see $q\propto\mathcal{P}$.

The dimensionless decay rate per Floquet cycle $\Gamma$ can then be obtained from $\Gamma(\epsilon,\tau_{1})=1/n^{*}(\epsilon,\tau_{1})$ where $n^{*}(\epsilon,\tau_{1})$ is the number of Floquet cycles such that $P(t)$ drops to some fixed threshold $1/A$. Solving $P(n^{*}\tau_{1})=1/A$ for $n^{*}$ yields

[TABLE]

where

[TABLE]

We take $A=10$; physically this corresponds to the situation where the polarization drops to $10\%$ of its starting value.

Lastly, to get the phase diagram, we estimated the phase boundary of the time crystal as the contour in the $\epsilon$-$\tau_{1}$ plane satisfying $\Gamma(\epsilon,\tau_{1})=\Gamma_{*}=1/100$. In other words, the phase boundary demarcates the regions in $\epsilon$-$\tau_{1}$ parameter space having significant decay after $100$ Floquet cycles or not. The choice of $100$ Floquet cycles was picked to match the experimental observations. Referring to Fig. 2 in the main text, one sees the linear ($\tau_{1}\propto|\epsilon|$) phase boundary for small $\tau_{1}$, and also the ‘closing up’ of the phase boundary, $\tau_{1}\propto 1/\epsilon^{2}$, for large $\tau_{1}$, as predicted by our theory.