Period-$n$ discrete time crystals and quasicrystals with ultracold bosons

TL;DR

This paper explores novel period-$n$ discrete time crystals and quasicrystals in ultracold bosons on a ring lattice, revealing complex dynamical phases including time quasicrystals, with potential experimental realization.

Contribution

It introduces a Floquet-driven model of interacting bosons exhibiting period-$n$ time crystals and quasicrystals, including a new incommensurate timescale and a detailed phase diagram.

Findings

Discovery of period-$n$ discrete time crystals in bosonic ring lattices.

Identification of a 'time quasicrystal' phase with incommensurate timescale.

Rich phase diagram including thermal and oscillatory phases.

Abstract

We investigate the out-of-equilibrium properties of a system of interacting bosons in a ring lattice. We present a Floquet driving that induces clockwise (counterclockwise) circulation of the particles among the odd (even) sites of the ring which can be mapped to a fully connected model of clocks of two counter-rotating species. The clock-like motion of the particles is at the core of a period-$n$ discrete time crystal where $L=2n$ is the number of lattice sites. In presence of a "staircase-like" on-site potential, we report the emergence of a second characteristic timescale in addition to the period $n$-tupling. This new timescale depends on the microscopic parameters of the Hamiltonian and is incommensurate with the Floquet period, underpinning a dynamical phase we call 'time quasicrystal'. The rich dynamical phase diagram also features a thermal phase and an oscillatory phase, all of which we investigate and characterize. Our simple, yet rich model can be realized with state-of-the-art ultracold atoms experiments.

Equations

$$H=\{\begin{array}{ll}-3J_1{\sum }_{j=1}^n(a_{2j}^{†}{a}_{2j-1}+h.c.),& 0<t<\frac{1}{3}\\ -3J_2{\sum }_{j=1}^n(a_{2j}^{†}{a}_{2j+1}+h.c.),& \frac{1}{3}<t<\frac{2}{3}\\ 3{\sum }_{j=1}^{2n}\left[\frac{U}{2N}{n}_j(n_j-1)+h_j{n}_j\right],& \frac{2}{3}<t<1\end{array}$$

$$\mathcal{I}(t)=\sum _{j=0}^{n-1}{e}^{i\frac{2\pi }{n}j}|{\psi }_{2j+1}(t){|}^2,$$

$$\tilde{\mathcal{I}}(\omega )=\underset{M\to \infty }{\lim }\frac{1}{M}\sum _{t=0}^{M-1}\mathcal{I}(t)e^{-i\omega t}.$$

$$Z(t)=e^{-i\frac{2\pi }{n}t}\mathcal{I}(t),$$

$$d^2(t)=\sum _{j=1}^L(|{\psi }_j(t){|}^2-|{\psi }_j^{\prime }(t){|}^2{)}^2.$$

$${\omega }_n^{1,2}=\frac{2\pi }{n}\pm \frac{2\arccos p_n(ϵ)}{n},$$

$$\frac{d{a}_j}{dt}=\{\begin{array}{ll}i3J_1{a}_{j\pm 1}& 0<t(\mathrm{mod}T)<\frac{1}{3}\\ i3J_2{a}_{j\mp 1}& \frac{1}{3}<t(\mathrm{mod}T)<\frac{2}{3}\\ -i3\left[\frac{U}{N}{n}_j+h_j\right]& \frac{2}{3}<t(\mathrm{mod}T)<1,\end{array}$$

$$\left(\begin{array}{c}{a}_{2j}\\ a_{2j-1}\end{array}\right)\left(t+1/3\right)$$

$$\left(\begin{array}{c}{a}_{2j}\\ a_{2j+1}\end{array}\right)\left(t+2/3\right)$$

$$\left(\begin{array}{c}{a}_{2j}\\ a_{2j-1}\end{array}\right)\left(t+1/3\right)$$

$$\left(\begin{array}{c}{a}_{2j}\\ a_{2j+1}\end{array}\right)\left(t+2/3\right)$$

$$n_{2j}(t+2/3)$$

$$n_{2j+1}(t+2/3)$$

$$n_j(t)=n_{j\mp 2t}(0),$$

$$F$$

$$={\Pi }_h\left[\sin (ϵ{)}^2{1}_{2n}-\cos (ϵ{)}^2{K}_2{K}_1-i\sin (ϵ)\cos (ϵ)(K_1+K_2)\right],$$

$$K_1=\left(\begin{array}{ccccccc}0& 1& & & & & \\ 1& 0& & & & & \\ & & 0& 1& & & \\ & & 1& 0& & & \\ & & & & \ddots & & \\ & & & & & 0& 1\\ & & & & & 1& 0\end{array}\right),K_2=\left(\begin{array}{ccccccc}0& & & & & & 1\\ & 0& 1& & & & \\ & 1& 0& & & & \\ & & & 0& 1& & \\ & & & 1& 0& & \\ & & & & & \ddots & \\ 1& & & & & & 0\end{array}\right),$$

$$({\tau }_e{)}_{2i,2j+1}=({\tau }_e{)}_{2i+1,2j}=0,$$

$$({\tau }_o{)}_{2i,2j+1}=({\tau }_o{)}_{2i+1,2j}=0,$$

$$({\tau }_e{)}_{2i,2j}=({\tau }_o{)}_{2i+1,2j+1}=\tau ,$$

$$\tau =\left(\begin{array}{cccccc}0& & & & 1& \\ 1& 0& & & & \\ & \ddots & \ddots & & & \\ & & 1& 0& & \\ & & & 1& 0& \end{array}\right).$$

$$F=\left(\begin{array}{cccccccccc}{\alpha }_1& {\gamma }_1& & & & & & & {\beta }_1& {\gamma }_1\\ {\gamma }_1& {\alpha }_1& {\gamma }_1& {\beta }_1& & & & & & \\ {\beta }_2& {\gamma }_2& {\alpha }_2& {\gamma }_2& & & & & & \\ & & {\gamma }_2& {\alpha }_2& {\gamma }_2& {\beta }_2& & & & \\ & & {\beta }_3& {\gamma }_3& {\alpha }_3& {\gamma }_3& & & & \\ & & & & & \ddots & & & & \\ & & & & {\beta }_{n-1}& {\gamma }_{n-1}& {\alpha }_{n-1}& {\gamma }_{n-1}& & \\ & & & & & & {\gamma }_{n-1}& {\alpha }_{n-1}& {\gamma }_{n-1}& {\beta }_{n-1}\\ & & & & & & {\beta }_n& {\gamma }_n& {\alpha }_n& {\gamma }_n\\ {\gamma }_n& {\beta }_n& & & & & & & {\gamma }_n& {\alpha }_n\end{array}\right).$$

$$F\vec{y}=\lambda \vec{y},$$

$$\begin{cases}(F\vec{y})_{2j}&=\omega_{j}\big{(}\gamma y_{j}^{o}+\alpha y_{j}^{e}+\gamma y_{j+1}^{o}+\beta y_{j+1}^{e}\big{)}=\lambda y_{j}^{e}\\ (F\vec{y})_{2j-1}&=\omega_{j}\big{(}\gamma y_{j}^{e}+\alpha y_{j}^{o}+\gamma y_{j-1}^{e}+\beta y_{j-1}^{o}\big{)}=\lambda y_{j}^{o}.\\ \end{cases}$$

$$\{\begin{array}{l}\gamma {\tilde{y}}_k^o+\alpha {\tilde{y}}_k^e+\frac{\gamma }{{\omega }_k}{\tilde{y}}_k^o+\frac{\beta }{{\omega }_k}{\tilde{y}}_k^e=\lambda {\tilde{y}}_{k-1}^e\\ \gamma {\tilde{y}}_k^e+\alpha {\tilde{y}}_k^o+\gamma {\omega }_k{\tilde{y}}_k^e+\beta {\omega }_k{\tilde{y}}_k^o=\lambda {\tilde{y}}_{k-1}^o,\end{array}$$

$$\lambda \left(\begin{array}{c}{\tilde{y}}_k^e\\ {\tilde{y}}_k^e\end{array}\right)=M_k\left(\begin{array}{c}{\tilde{y}}_{k-1}^e\\ {\tilde{y}}_{k-1}^e\end{array}\right),$$

$$M_k=\left(\begin{array}{cc}\sin (ϵ{)}^2-\frac{\cos (ϵ{)}^2}{{\omega }_k}& -i\sin (ϵ)\cos (ϵ)(1+\frac{1}{{\omega }_k})\\ -i\sin (ϵ)\cos (ϵ)(1+{\omega }_k)& \sin (ϵ{)}^2-\cos (ϵ{)}^2{\omega }_k\end{array}\right).$$

$$\left(\begin{array}{c}{\tilde{y}}_k^e\\ {\tilde{y}}_k^e\end{array}\right)={\lambda }^n{M}_{k+1}^{-1}{M}_{k+2}^{-1}{M}_{k+3}^{-1}\dots M_{k+n-1}^{-1}{M}_k^{-1}\left(\begin{array}{c}{\tilde{y}}_k^e\\ {\tilde{y}}_k^e\end{array}\right),$$

$${\lambda }_j^{\pm }=\exp \left[i\frac{\pi +2\pi j\pm \arccos p_n(ϵ)}{n}\right],$$

$$p_n(ϵ)=-\frac{1}{2}\sum _{j_1,j_2,\dots ,j_n=1}^2(M_n{)}_{j_1,j_2}(M_{n-1}{)}_{j_2,j_3}\dots (M_1{)}_{j_n,j_1}$$

$$M_1=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),M_2=\left(\begin{array}{cc}\sin (ϵ{)}^2-\cos (ϵ{)}^2& -i2\sin (ϵ)\cos (ϵ)\\ -i2\sin (ϵ)\cos (ϵ)& \sin (ϵ{)}^2-\cos (ϵ{)}^2\end{array}\right),$$

$$p_2(ϵ)=$$

$$p_3(ϵ)=$$

$$p_4(ϵ)=$$

Full text

Introduction.— Symmetries pervade most fields of modern physics, ranging from nuclear physics to relativity and condensed matter physics. The spontaneous breaking of symmetries is the essential mechanism at the core of equilibrium phase transitions. Quite surprisingly, the possibility of systems breaking time-translational symmetry has been put forward only recently by Wilczek Shapere and Wilczek (2012); Wilczek (2012) who dubbed them ’time crystals’ and triggered an extraordinary amount of excitement Sacha and Zakrzewski (2017); Moessner and Sondhi (2017); Else et al. (2019). Since time crystals have been shown to be impossible in quantum ground states Bruno (2013), research has focussed on out-of-equilibrium conditions. Particularly successful has been the setting of Floquet systems, characterized by a time-periodic Hamiltonian $H(t)=H(t+T)$, for which the notion of a discrete time crystal (DTC) has been introduced in Refs. Sacha (2015a); Else et al. (2016); von Keyserlingk and Sondhi (2016); Khemani et al. (2016). A DTC is a system that, for some physical observable $O$ with expectation value $f(t)=\expectationvalue{O}{\psi(t)}$, in the thermodynamic limit and for a set of initial conditions $\ket{\psi(0)}$, features three properties Russomanno et al. (2017): (I) discrete time-translational symmetry breaking: $f(t)\neq f(t+T)$; (II) rigid subharmonic response: without fine-tuning $f(t)$ shows oscillations with a period $nT$, i.e. period-$n$-tupling, that is the dynamics features a characteristic frequency $\frac{2\pi}{n}$ with $n$ integer $\geq 2$; (III) persistence: the subharmonic oscillations extend up to infinite time.

Among the plethora of proposed DTCs, most feature period-doubling $n=2$ Sacha (2015a); Else et al. (2016); von Keyserlingk and Sondhi (2016); Khemani et al. (2016); Moessner and Sondhi (2017); Ho et al. (2017); Else et al. (2017); Iemini et al. (2018); Zhu et al. (2019); Russomanno et al. (2017); Gong et al. (2018); Yao et al. (2017); Zhang et al. (2017); Choi et al. (2017); Rovny et al. (2018); Smits et al. (2018). The possibility of a ”period-$n$ DTC ” ($n>2$), which has been discussed for $n$-hands clock models Sreejith et al. (2016); Surace et al. (2018), still lacks a physical implementation Giergiel et al. (2018a). On the other hand, investigations of DTCs for ultracold bosons have for the most part remained restricted to the context of a vibrating mirror in presence of a gravitational field Huang et al. (2018a); Sacha (2015a, b); Giergiel et al. (2018b, a). In light of the terrific experimental progress and control with cold atoms in optical lattices Bloch et al. (2008); Smits et al. (2018), proposals based on this platform would be very desirable.

In this letter we study the dynamical properties of bosons in a ring lattice. For a simple Floquet driving protocol of nearest-neighbor hopping and local interaction, we find period-$n$ DTCs ($n\geq 2$), thermal, and oscillatory phases. Surprisingly, in the presence of a staircase-like on-site potential, a new dynamical phase emerges characterized by three incommensurate frequencies (Fig. 1). The dynamics is strictly aperiodic, but long-time ordered: a discrete time quasicrystal (DTQC) Lifshitz (2003, 2011).

The idea of time quasi-crystallinity has recently appeared in the literature with various connotations and in various contexts such as dissipative classical systems Flicker (2018), finite-size systems Huang et al. (2018b), and quasiperiodically driven systems Dumitrescu et al. (2018); Peng and Refael (2018); Zhao et al. (2019). The term ’time quasicrystal’ has also been adopted in an experiment with magnons in which one of the system’s characteristic frequencies was incommensurate with the driving frequency Autti et al. (2018) and in a model where a periodic repetition in time of a (possibly large, but finite) Fibonacci word was observed Giergiel et al. (2018b).

In striking contrast to the cases above, we consider a quantum, macroscopic, and periodically driven system and we find a DTQC phase that features a characteristic subharmonic frequency $\frac{2\pi}{n}$ and two other, generally incommensurate, intrinsic frequencies $\omega_{n}^{1,2}$. The first frequency is locked to a submultiple of the driving frequency, breaking discrete time-translational symmetry and robust to perturbations, similar to a standard DTC, whereas the latter two frequencies depend on the microscopic parameters.

The remainder of this letter is organized as follows. First, we introduce our model and study its solvable limits. We then propose a set of suitable diagnostic dynamical order parameters to characterize the system. We study the dynamics first in absence of an on-site potential and find a period-$n$ DTC, a thermal, and an oscillatory phase. In the presence of a staircase-like on-site potential, we additionally find a DTQC phase. We demonstrate its rigidity, and find analytical expressions for the incommensurate frequencies $\frac{2\pi}{n}$ and $\omega_{n}^{1,2}$. Finally, we summarize our results and outline possible experimental implementations using state-of-the-art ultracold atom setups.

Model and integrable limits.— We consider a system of $N$ bosons on a one-dimensional ring lattice with $L=2n$ sites governed by the Floquet Hamiltonian with period $T=1$ Mizuta et al. (2018)

[TABLE]

where we set $J_{1}=J_{2}=\frac{\pi}{2}+\epsilon$ and $h_{2j-1}=h_{2j}=\frac{2\pi}{n}j\delta$, that is a staircase-like on-site potential. The ring structure implicitly corresponds to periodic boundary conditions $a_{j}=a_{j+L}$. The Floquet Hamiltonian (1) alternates hopping between odd nearest-neighbor links, hopping between even nearest-neighbor links, and on-site potential and two-body interaction. In the following $t=0,1,2,\dots$ denotes stroboscopic times.

We gain some intuition in the dynamical properties of the system by solving the Heisenberg equation of motion $\frac{da_{j}}{dt}=i[H(t),a_{j}(t)]$ ($\hbar=1$) in the two cases $\epsilon=0$ and $U=0$. The case $\epsilon=0$ is at the core of a period-$n$ DTC and is schematically illustrated in Fig. 1. Consider an initial product state in the site basis with $n_{j}$ bosons in an odd (even) site $j$. During the first part of the driving these particles move to site $j+1$ ($j-1$), during the second they move to site $j+2$ ($j-2$), and during the third the number of particles in each site is conserved (irrespective of $U$ and $\{h_{j}\}$). We thus find $n_{j}(t)=n_{j\pm 2t}(0)$ with $+$ and $-$ for even and odd sites $j$, respectively. The bosons in the ring lattice can be shown to be equivalent to a fully connected model of clocks with $n$-hands and of two species: clockwise and counterclockwise rotating (see SM). The clocks tick at every Floquet period, so that after a time $t=n$ the system returns to its initial condition, a mechanism at the core of a period-$n$ DTC.

In the non-interacting limit ($U=0$), the Heisenberg equation of motion is linear and easily solved as $\vec{a}(t+1)=F\vec{a}(t)$, where $\vec{a}=(a_{1},a_{2}\dots,a_{2n})^{T}$ and where $F$ is a $2n\times 2n$ dimensional matrix. The dynamics is thus characterized by oscillations at the $2n$ frequencies corresponding to the phases of the $2n$ eigenvalues $\{\lambda_{j}\}$ of $F$. For instance, for $\epsilon=0$ we find doubly degenerate eigenvalues of the form $\lambda_{j}^{1,2}=e^{i\phi+i\frac{2\pi}{n}j}$ (with $\phi$ some non-relevant phase), which correctly signals the period-$n$ DTC with clock-like clockwise (counterclockwise) rotation of particles.

Dynamical order parameters.— We are now interested in the dynamics away from the solvable limits $\epsilon=0$ and $U=0$. To this end, we solve a semiclassical Gross-Pitaevskii equation of motion, which is obtained from the Heisenberg equation upon replacing the bosonic operators with $c$-numbers $\vec{a}\rightarrow\sqrt{N}\vec{\psi}$. In the thermodynamic limit of macroscopic occupation ($N\rightarrow\infty$) for a fixed finite ring size $L$ Sacha and Zakrzewski (2017), considering a symmetry broken initial state with all the particles in site $j=1$ (i.e. $\psi_{j}(0)=\delta_{j,1}$), the semiclassical dynamics can either be chaotic, signaling quantum thermalization Cosme and Fialko (2014) or not, in which case we expect it to become exact Polkovnikov (2003). We indeed confirm this for our model explicitly using exact diagonalization and finite-size scaling (see SM).

It is now crucial to introduce suitable dynamical order parameters. To track the clock-like circulation of the particles we introduce a generalized imbalance on the odd sites

[TABLE]

and consider its Fourier transform

[TABLE]

A measure of the time crystallinity is given by

[TABLE]

whereas the chaoticness of the semiclassical dynamics is quantified by the distance $d^{2}$ between two slightly different initial states $\vec{\psi}(0)=(1,0,0,\dots,0)^{T}$ and $\vec{\psi}^{\prime}(0)=(1-\Delta,\Delta,0,\dots,0)^{T}$ with an arbitrary small $\Delta=10^{-10}$

[TABLE]

A growth of $d^{2}(t)$ to a finite value $\sim 1$ corresponds to classical sensitivity to initial conditions and signals quantum thermalization Cosme and Fialko (2014). Finally, we consider the infinite time averages $\langle d^{2}\rangle_{t}=\lim\limits_{M\rightarrow\infty}\frac{1}{M}\sum_{t=0}^{M-1}d^{2}(t)$ and $\langle Z\rangle_{t}=\tilde{\mathcal{I}}\left(\frac{2\pi}{n}\right)$.

In the integrable limit $\epsilon=0$ we have a DTC, the phase of $\mathcal{I}$ grows linearly in time, $\tilde{\mathcal{I}}(\omega)$ is peaked at $\omega=\frac{2\pi}{n}$, $\langle Z(t)\rangle_{t}=1$ and $\langle d^{2}\rangle_{t}\approx 0$. In the following, we will use these quantities to characterize the dynamical phases of the system.

Dynamical phases without on-site potential.— We start by considering a vanishing on-site potential ($\delta=0$). In Fig. 2 we focus on the example $L=6$, but similar results can be obtained for any $L=2n$ (see SM). We find that three dynamical phases are possible: (I) period-$n$ DTC: for small $\epsilon\approx 0$ the clock-like circulation of the particles is rigidified by the interaction $U\neq 0$ and signaled by $Z\approx 1$. The characteristic frequency of the imbalance $\mathcal{I}$ is locked to $\omega=\frac{2\pi}{n}$, i.e. the frequency of the peak of $\tilde{\mathcal{I}}$ is robust against perturbations. This is the quintessential property of a DTCs: the interactions make the system subharmonic response rigid to mistakes in the driving, that is the range of $\epsilon$ corresponding to the DTC grows with $U$. We note that, mathematically, a robust subharmonic response is possible only if the semiclassical dynamical equation is nonlinear, which occurs for non-vanishing interactions $U\neq 0$ and confirms the many-body nature of the DTC. (II) Oscillatory phase: for a large $\epsilon$ but small interaction $U$ the system exhibits non-ergodic oscillations with a few characteristic frequencies signaled by sharp peaks in $\tilde{\mathcal{I}}(\omega)$. In contrast to the case of a DTC, these frequencies are not locked, but rather depend on $\epsilon$. (III) Thermal phase: for large $\epsilon$ and $U$ classical chaos emerges, i.e. $\langle d^{2}\rangle_{t}\sim 1$. The three phases touch in the tricritical point $\epsilon=U=0$.

We emphasize that the observation of the DTC is not exclusive of the initial condition we considered here but is rather valid for a generic initial Fock state featuring a macroscopic imbalance, thanks to the underlying clock-like particle circulation. Also, we notice that different classes of initial conditions are expected to correspond to DTCs with different periods. For instance, an initial state where one site every $2m$ is equally occupied would lead to period $m$-tupling (independently of $L$). In this case, a finite-density occupation could still realize a DTC for $L\rightarrow\infty$ Mizuta et al. (2018). Moreover, in general our model does not require any disorder since in a fully connected system (to which our model can be mapped) disorder is not a necessary ingredient to prevent thermalization under Floquet driving Russomanno et al. (2017); Surace et al. (2018). This is due to an extensive number of integrals of motion as explained in detail in the SM. Finally, we note that the effective full-connectivity of the equivalent clock model ultimately emerges from the symmetry of the multi-particle bosonic wavefunction of the considered physical system of cold atoms, and is in this sense an intrinsic property of the system which cannot be broken by perturbations to the bosonic Hamiltonian.

The realization of period-$n$ DTC in a simple bosonic model in a ring lattice is the first major finding of this work.

Tilted lattice and discrete time quasicrystal.— Next, we show that tilting the lattice with a staircase-like potential ($\delta\approx 1$) favors a new DTQC phase. In presence of such a potential, we gain a clear intuition of the dynamics solving exactly the limit of $U=0$ and $\delta=1$ (with $\epsilon$ generally $\neq 0$). In such a limit, the system’s characteristic frequencies are linked to the eigenvalues of $F$ and turn out to be (details in the SM) $\frac{2\pi}{n}$, $\omega_{n}^{1}$ and $\omega_{n}^{2}$. Only the first frequency is locked to a submultiple of the driving frequency, whereas the other two are generally incommensurate with it and are given by

[TABLE]

where $p_{n}(\epsilon)$ is a trigonometric polynomial in $\epsilon$ with $p_{n}(0)=1$, e.g. $p_{2}(\epsilon)=\cos 2\epsilon$ and $\omega_{2}^{1,2}=\pi\pm 2\epsilon$ for $L=4$. In the limit $\epsilon=0$ we find $\omega_{n}^{1,2}=\frac{2\pi}{n}$, which consistently links back to the period-$n$ DTC. In the thermodynamic limit – in our case corresponding to a macroscopic number of bosons $N\to\infty$ but finite $L$ – this regime represents a proper dynamical phase, which we will refer to as a DTQC Lifshitz (2003, 2011). The discovery of a DTQC and of its characteristic frequencies (Eq. (6)) represents the second major finding of this work. We emphasize that, differently from previous studies Dumitrescu et al. (2018); Peng and Refael (2018); Zhao et al. (2019), the driving that we consider is perfectly periodic and characterized by a single frequency whereas the system response features a subharmonic frequency $\frac{2\pi}{n}$ and two intrinsic frequencies $\omega_{n}^{1,2}$ that depend on the microscopic parameters of the system.

We show that these results are not an artifact of fine-tuning, but rather underline a proper dynamical phase. A non-zero interaction ($U\neq 0$) indeed rigidifies the DTQC. This means that, even if the frequencies $\omega_{n}^{1,2}$ depend on the model parameters (particularly on $\epsilon$), the characteristic subharmonic frequency $\frac{2\pi}{n}$ does not shift. For $\delta\approx 1$ we observe the new DTQC phase, which is characterized by three main sharp peaks in $\tilde{\mathcal{I}}$ (Fig. 3a). One peak is locked to $\frac{2\pi}{n}$ and does not change for $\epsilon$ and $\delta$ within a certain range, whereas the other two are well approximated by the non-interacting prediction of Eq. (6). We emphasize that, differently from previous works Autti et al. (2018), we believe the presence of two incommensurate frequencies to be a necessary but not sufficient ingredient for the definition of a DTQC 111A definition relying solely on the presence of two incommensurate frequencies would in fact classify a single, undamped, driven harmonic oscillator as a DTQC whenever the driving frequency $\omega_{d}$ is incommensurate with the oscillator’s characteristic frequency $\omega_{c}$.. Rather, in analogy with the more established notion of DTC, we judge essential the presence of the robust subharmonic peak, which is a genuine many-body feature. Beyond the DTQC, we also still observe the DTC, oscillatory, and chaotic phases, which are reflected in $\tilde{\mathcal{I}}(\omega)$ featuring a single peak locked to $\frac{2\pi}{n}$, few sharp peaks and no sharp peaks at all, respectively. In Fig. 3b we sketch the various phases as a function of $\delta$ and $U$ for $\epsilon=0.1$ by looking at $\langle Z\rangle_{t}=\tilde{\mathcal{I}}(\frac{2\pi}{n})$ and $\langle d^{2}\rangle_{t}$. The interaction $U$ rigidifies the DTQC, which expands from $\delta=1$ for increasing $U$ and is signaled by $\langle Z\rangle\approx 0.5$. For increasing interactions the system enters a thermal and/or a DTC phase. In Fig. 3c we characterize the system as a function of $\epsilon$ and $U$ for $\delta=0.9$. We see that the interaction necessary to enter the DTC grows with $\epsilon$. Eventually, if $\epsilon$ is too large, no DTC is possible.

Experimental implementations.— An appealing feature of our proposal lies in its ultimate simplicity making it very natural for experimental implementation. Indeed, the Hamiltonian (1) only requires a number of basic ingredients such as nearest-neighbor hopping, local interaction and on-site potential, which are all readily available for ring-shaped optical lattices Franke-Arnold et al. (2007); Houston et al. (2008); Vyas et al. (2013); Amico et al. (2014). The main difficulty is the discrete time dependent switch of tunnelling between alternating bonds. For small systems this is easily achievable if alternating bonds point in different real space directions. Alternatively, the even and odd site labels could be imprinted onto two different spin states. This would enable the discrete switch also on larger ring lattices. More concretely, the dynamics of the first and of the second third of the Floquet Hamiltonian would in this case correspond to a transverse field-induced spin-flip and to a intra-well barrier lowering, respectively. Furthermore, we have also checked that our results are robust against small undesired hoppings at odds with the alternating Floquet protocol.

Conclusions.— In conclusion, we have studied the non-equilibrium properties of interacting bosons in a ring lattice. Considering a Floquet driving which alternates hopping on even and odd nearest-neighbor links, we induced a clock-like circulation of the particles whose direction depends on the parity of sites. Introducing a suitable set of dynamical order parameters and solving the exact dynamics in two integrable limits and a Gross-Pitaevskii equation across the whole parameter space, we identified several dynamical phases: a period-$n$ DTC, an oscillatory phase, a thermal phase, and a new DTQC. The DTQC phase is characterized by a frequency which is locked to a submultiple of the Floquet frequency and by a second frequency which is incommensurate with it. Our work demonstrates a wide spectrum of non-equilibrium, non-trivial dynamical phases of bosons, which represent a unique opportunity for experimental investigation.

As natural for bosons in a finite-size lattice, the proposed system can effectively be mapped onto a fully-connected model. A natural question for future investigation regards then the fate of the various dynamical phases in physical (rather than effective) systems of clocks where such a full-connectivity is broken, e.g. in presence of long-range power-law interaction. Moreover, our work adds a totally new perspective to the multiple connotations of the still developing concept of time quasi-crystallinity, which certainly deserves further study. For instance, an intriguing question concerns the possibility for a system driven with two incommensurate frequencies $\omega_{1}$ and $\omega_{2}$ to respond at subharmonic frequencies $\omega_{1}/n_{1}$ and $\omega_{2}/n_{2}$.

I I) Dynamical equations and solvable limits

We derive explicitly the quantum dynamical equation and solve it in the limit cases $\epsilon=0$ and $U=0$. The Heisenberg dynamical equations for the annihilation operator $a_{j}$ at site $j=1,2,\dots,L=2n$ is obtained with a straightforward computation of the commutators of the Hamiltonian (1) with $a_{j}$ and reads ($\hbar=1$)

[TABLE]

where the upper and the lower signs are for odd and even $j$, respectively. In the following we implicitly assume a stroboscopic time $t=0,1,2,\dots$, unless differently specified. Recalling $J_{1}=J_{2}=\frac{\pi}{2}+\epsilon$, we integrate the dynamics over the first two fractions of the Floquet driving to obtain

[TABLE]

whereas the third fraction of the Floquet driving can be integrated only in the two limit cases $\epsilon=0$ and $U=0$.

I.1 Limit $\epsilon=0$

For a fine-tuned hopping strength with $\epsilon=0$, we get

[TABLE]

and thus

[TABLE]

Moreover, since $[H,n_{j}]=0$ for $2<t<3$, we finally get

[TABLE]

where again the upper and lower sign refer to odd and even $j$, respectively. We thus obtain that $n_{j}(t=n)=n_{j}(t=0)$, that is the solvable limit $\epsilon=0$ is at the core of a period-$n$ DTC.

I.2 Limit $U=0$

In the non-interacting limit ($U=0$), also the dynamics associated to the third fraction of the Floquet period becomes linear, and can be trivially solved by $a_{j}(t+1)=e^{-ih_{j}}a_{j}(t+2/3)$. The evolution of the bosonic operators $\vec{a}=(a_{1},a_{2},\dots,a_{2n})^{T}$ over one period can be compactly written as $\vec{a}(t+1)=F\vec{a}(t)$ with

[TABLE]

where $\mathbb{1}_{2n}$, $K_{1}$, $K_{2}$ and $\Pi_{h}$ are $2n\times 2n$-dimensional matrices. In particular, $\Pi_{h}$ is diagonal and with entries $(\Pi_{h})_{j,j}=e^{-ih_{j}}$, $\mathbb{1}_{2n}$ is the identity and $K_{1}$ and $K_{2}$ are $2n\times 2n$-dimensional involution matrices swapping sites $(1,2),(3,4),\dots,(L-1,L)$ and sites $(2,3),(4,5),\dots,(L,1)$, respectively

[TABLE]

It is easy to check that $K_{2}K_{1}=\tau_{e}^{\dagger}+\tau_{o}$ where $\tau_{e}$ and $\tau_{o}$ are defined by

[TABLE]

with $\tau$ the following $n\times n$-dimensional matrix

[TABLE]

The matrix $K_{2}K_{1}$ thus generates a clockwise (counterclockwise) circulation of the particles within the odd (even) sites.

In general, the system dynamics is characterized by $2n$ frequencies which are linked to the eigenvalues of $F$. Considering the staircase-like potential $h_{2j}=h_{2j-1}=\frac{2\pi}{n}j\delta$, these eigenvalues take a particular form in the limit $\delta=1$, which we are about to show and which is at the core of the DTQC phase.

Let us say $\omega_{j}=e^{i\frac{2\pi}{n}j}$, $\alpha_{j}=\omega_{j}\sin(\epsilon)^{2}$, $\beta_{j}=-\omega_{j}\cos(\epsilon)^{2}$, $\gamma_{j}=-i\omega_{j}\cos(\epsilon)\sin(\epsilon)$ and write $F$ as

[TABLE]

The eigenvalue problem reads

[TABLE]

where $\vec{y}=(y_{1}^{o},y_{1}^{e},y_{2}^{o},y_{2}^{e},\dots,y_{n}^{o},y_{n}^{e})^{T}$ is a $2n$-dimensional column vector. We rewrite the eigenvalue problem (15) component by component as

[TABLE]

Dividing both members of Eqs. (16) by $\omega_{j}$ and introducing the Fourier transform $\tilde{y}_{k}^{e/o}=\sum_{j=1}^{n}\omega_{j}y_{j}^{e/o}$ (with $k=1,2,\dots,n$) we get

[TABLE]

that we rewrite in a compact form as

[TABLE]

where $M_{k}$ reads

[TABLE]

It is easy to show that $\det(M_{k})=1$ for every $k$, so that $M_{k}$ admits an inverse $M_{k}^{-1}$. Inverting Eq. (18), iterating it and exploiting periodicity in momentum space ($\tilde{y}_{k}^{e/o}=\tilde{y}_{k\pm n}^{e/o}$) we get

[TABLE]

which is itself a $2\times 2$-dimensional eigenvalue problem for the eigenvalue $\lambda^{n}$. The matrix $M_{n}M_{n-1}\dots M_{1}$ has determinant $1$, and its eigenvalues can therefore be written as $-e^{\pm i\acos(-\frac{\Tr}{2})}$ where $\Tr$ is its trace. For $\delta=1$, we thus finally get the eigenvalues of $F$ to be

[TABLE]

where

[TABLE]

is in general a trigonometric polynomial in $\epsilon$, of order $2n$ and such that $p_{n}(0)=1$. For instance, for $L=2n=4$ we obtain

[TABLE]

so that $\Tr=2\sin(\epsilon)^{2}-2\cos(\epsilon)^{2}$, i.e. $p_{2}(\epsilon)=\cos(2\epsilon)$. We plot $p_{n}(\epsilon)$ for $n=2,3,\dots,10$ in Fig. S1a, and report here its explicit expression for $n=2,3,\dots,6$

[TABLE]

As explained in the main text, the eigenvalues (21) are at the core of a DTQC dynamical phase. The characteristic frequencies $\omega_{n}^{1,2}$ are plotted in Fig. S1b. In the limit $\epsilon\rightarrow 0$ we get $\omega_{n}^{1,2}\rightarrow\frac{2\pi}{n}$, consistently recovering the DTC.

II II) Map to a fully connected clock model

Our bosonic model is equivalent to a model of fully connected (i.e. infinite-range interacting) clock variables of two counter-rotating species, as we show here in the limit $\epsilon=0$. For simplicity we adopt the inverse route: we present the clock model and map it back to the bosonic one, following in spirit the mapping that Ref. Surace et al. (2018) carried out for clocks of a single species.

Consider a system of $N$ clock variables with $n$-hands, which generalizes the spin case corresponding to $n=2$. To each clock we further associate a binary label $s=e,o$ indicating the species of the clock. We will refer to clocks of species $e,o$ as to even and odd clocks, respectively. A natural basis of the Hilbert space is the one of states $\ket{\{s_{i},j_{i}\}}=\bigotimes_{i=1}^{N}\ket{s_{i},j_{i}}_{i}$ with the $i$-th clock being of species $s_{i}$ and in the $j_{i}$-th hand. We introduce the local operators $\sigma_{s,i}$, $\eta_{s,i}$ and $\tau_{s,i}$ acting on the $i$-th clock and defined as

[TABLE]

where $\bar{s}=e,o$ for $s=o,e$, respectively, $\omega=e^{i\frac{2\pi}{n}}$, and where $\{h_{s,j}\}$ are so far unspecified numbers. The operator $\tau_{s,i}$ acts on the $i$-th clock moving its hand one step forward, if of species $s$, or annihilating it, is of species $\bar{s}$. Consider a Floquet operator $U_{F}=e^{-iH}U_{K}$, where $H$ is the Hamiltonian

[TABLE]

and $U_{K}$ is a local kicking operator of the form

[TABLE]

The first term of the Hamiltonian (33) plays the role of interaction among clocks of the same species and couples all sites, making the system fully connected. In the chosen basis, $H$ is diagonal whereas the operator $U_{K}$ acts moving one step forward (backward) the hands of the odd (even) clocks. After $n$ Floquet periods, the state of the system will therefore come back to the initial condition. This mechanism is at the basis of a period-$n$ DTC.

Being fully connected, the system can be described just counting the number of clocks of a certain species and with hands pointing in a certain direction. That is, it is possible to describe the system in terms of bosonic operators $b_{j},b_{j}^{\dagger}$ with $j=1,2,\dots,2n$ and fulfilling the standard bosonic commutation relations $[b_{j},b_{j^{\prime}}^{\dagger}]=\delta_{j,j^{\prime}}$ and $[b_{j},b_{j^{\prime}}]=0$ (for further details we refer to Appendix D of Surace et al. (2018)). We say

[TABLE]

where $P=\sum_{i_{1},i_{2},\dots,i_{N}}\Pi_{i_{1},i_{2},\dots,i_{N}}$ is the symmetrization operator, with sum running over the possible permutations of $1,2,\dots,N$ and with $\Pi_{i_{1},i_{2},\dots,i_{N}}$ being the corresponding permutation operator. The bosonic operators are such that

[TABLE]

Since the Hamiltonian (33) is invariant under site permutations, it is easy to show that

[TABLE]

where $n_{j}=b_{j}^{\dagger}b_{j}$ is the number operator for the $j$-th bosonic mode. From Eq. (37), performing the sum over $m$, removing a constant term $U$ and saying $h_{2j-1}=h_{o,j}$ and $h_{2j}=h_{e,j}$, we finally obtain

[TABLE]

As well, we find that the kickting operator acts as

[TABLE]

Summing up, we have therefore shown that the fully-connected clock model can be mapped into the model of bosons in a ring with $2n$ sites. The effect of the kicking operator is to rotate the bosons in the odd (even) sites of two steps in the clockwise (countercockwise) direction, whereas the clock Hamiltonian (33) maps into a local potential and two-body interaction for bosons. In the main, the kicking operator is then realized thanks to the sequential action of first and the second fractions of the ternary Floquet Hamiltonian.

The effective full-connectivity of the bosonic model has an important implication: the system is invariant under the permutation of clocks $(1,2)$, $(2,3)$, $\dots,(N-1,N)$, which corresponds to an extensive number $\lfloor\frac{N}{2}\rfloor$ of integrals of motion. Thanks to the existence of these integrals of motion, a priori the system does not necessitate any disorder to escape the fate of thermalization, making non-trivial dynamical phases such as DTCs possible.

III III) On the validity of mean field

Throughout the main text, numerical results are obtained solving the Gross-Pitaevskii Equation (GPE), which is often referred to as the mean field (MF) limit of the Heisenberg equation. Here we discuss the regimes in which this is legitimate, and how to interpret the results when it is not.

In the considered thermodynamic limit ($N\rightarrow\infty$ for a fixed $L$), the correct and well-established phase-space framework to deal with non-equilibrium dynamics is the Truncated Wigner Approximation (TWA) Polkovnikov (2003). This computes expectation values as averages over the classical GPE trajectories obtained for an ensemble of classical initial conditions, whose stochastic distribution is chosen consistently with the actual quantum initial condition. When the system is initialized in $\ket{\psi(0)}=\ket{N,0,\dots,0}$, the initial condition for the GPE reads $\psi_{j}(0)=\delta_{j,1}e^{i\theta_{j}}$ and, thanks to the gauge symmetry, can actually be considered without loss of generality to be $\psi_{j}(0)=\delta_{j,1}$. Since the initial condition no longer contains any degree of freedom, in this case the TWA is equivalent to a ”single-shot” Gross-Pitaevskii Equation (SSGPE). This observation sounds very promising from a computational point of view, because it allows avoiding running the GPE multiple times, but it comes with a drawback: the SSGPE is not guaranteed to be accurate since the initial condition features non-macroscopically occupied sites. In particular, we expect the SSGPE to be inaccurate when chaotic. In the framework of TWA, chaotic semiclassical equations are in fact expected to underline thermalization Cosme and Fialko (2014), whereas the SSGPE displays in general persistent (in fact, chaotic) oscillations.

III.1 Exact diagonalization

To support the previous claims we compare the SSGPE with exact diagonalization (ED) computations for $L=4$ sites, finite $N\leq 28$ and vanishing on-site potential ($\delta=0$). For various values of $\epsilon$ and $U$ corresponding to the DTC, oscillatory and thermal phases, in Fig. S2a1-a3 we plot $Z(t)$. The DTC is characterized by $Z\approx 1$ already for a relatively small $N$, whereas in the oscillatory phase the accuracy of the MF solution is worse and deteriorates in time. Nevertheless, in both cases, for increasing $N$ the accuracy of MF improves, and in the thermodynamic limit $N\rightarrow\infty$ (for a fixed, finite $L$) we expect MF to be exact up to arbitrarily large times. Conversely, in the thermal phase the MF solution is inaccurate irrespective of $N$. In this phase, $Z_{MF}$ fluctuates chaotically whereas $Z_{ED}$ saturates to a steady value, i.e. thermalizes, as expected Cosme and Fialko (2014).

To better interpret the results, we check the finite-size scaling considering the parameter

[TABLE]

which serves as a measure of the MF error with respect to ED for a given $N$. In Fig. S2b we show that $\mathcal{E}$ decays as a power of $N$ in the DTC and oscillatory phases, whereas it does not decay in the thermal phase. This supports the idea for which, in the non-chaotic regime and in the thermodynamic limit, the SSGPE is exact. On the other hand, a failure of the SSGPE in the chaotic regime corresponds to the onset of quantum thermalization Cosme and Fialko (2014).

IV IV) Larger rings

In the main paper we reported numerical results for $L=4,6$. However, our findings are valid for an arbitrary even number of sites $L=2n\geq 4$. To support this claim, here we corroborate our results studying a system with $L=8$ sites, for which the time crystalline phases are characterized by period $4$-tupling. To sketch the dynamical phases we look as usual at the height of the Fourier subharmonic peak $\tilde{\mathcal{I}}(\frac{2\pi}{n})=\langle Z\rangle_{t}$ and at the distance $\langle d^{2}\rangle_{t}$ between two initially very close copies of the system, which are shown in Fig. S3. The results that we find are completely analogue to the ones of the main text.