Observation of tunable discrete time crystalline phases

TL;DR

This paper reports the experimental observation of multiple discrete time crystalline phases in a nanoelectromechanical system, demonstrating their many-body nature, long-range order, and tunability via mechanical strain, advancing understanding of non-equilibrium phases.

Contribution

It presents the first observation of multiple DTC phases in a NEMS device, showing their controllability and complex phase diagram, which broadens the scope of time crystal research.

Findings

Multiple DTC phases observed in NEMS system

Confirmation of many-body and long-range order

Controlled phase transitions via mechanical strain

Abstract

Discrete time crystals (DTCs) are emergent non-equilibrium phases of periodically driven many-body systems, with potential applications ranging from quantum computing to sensing and metrology. There has been significant recent interest in understanding mechanisms leading to DTC formation and a search for novel DTC phases beyond subharmonic entrainment. Here, we report observation of multiple DTC phases in a nanoelectromechanical system (NEMS) based on coupled graphene and silicon nitride membranes. We confirm the time-crystalline nature of these symmetry broken phases by establishing their many-body characters, long-range time and spatial order, and rigidity against parameter fluctuation or noise. Furthermore, we employ controlled mechanical strain to drive the transitions between phases with different symmetries, thereby mapping the emergent time-crystalline phase diagram. Overall, our work takes a step towards establishing time crystals as a system with complexity rivaling that of solid state crystals.

Equations

$$C_{i,j}(\tau )=|cov(Z_i,Z_j(\tau ))/{\sigma }_{Z_i}{\sigma }_{Z_j}(\tau )|.$$

$$\mathcal{H}=\frac{p_1^2}{2m_g}+\frac{p_2^2}{2m_g}+\sum _{i=3}^n\frac{p_i^2}{2m_s}+\frac{1}{2}{m}_g{\omega }_1^2{x}_1^2+\frac{1}{2}{m}_g{\omega }_2^2{x}_2^2+\frac{1}{2}{m}_s\sum _{i=3}^n{\omega }_i^2{x}_i^2+\frac{1}{4}\beta (x_1^4+x_2^4)-\sum _{i=3}^n{\alpha }_{1i}^{\prime }{x}_1{x}_i-\sum _{i=3}^n{\alpha }_{2i}^{\prime }{x}_2{x}_i.$$

$${\ddot{x}}_1+{\omega }_1^2{x}_1+\frac{\beta }{m_g}{x}_1^3-\sum _{i=3}^n\frac{{\alpha }_{1i}^{\prime }}{m_g}{x}_i=0,$$

$${\ddot{x}}_2+{\omega }_2^2{x}_2+\frac{\beta }{m_g}{x}_2^3-\sum _{i=3}^n\frac{{\alpha }_{2i}^{\prime }}{m_g}{x}_i=0,$$

$${\ddot{x}}_i+{\omega }_i^2{x}_i-\frac{{\alpha }_{1i}^{\prime }}{m_s}{x}_1-\frac{{\alpha }_{2i}^{\prime }}{m_s}{x}_2=0,\forall i=\{3,4,\cdots ,n\}.$$

$${\ddot{x}}_1+{\omega }_1^2{x}_1+\frac{\beta }{m_g}{x}_1^3-\frac{\alpha }{m_g}{x}_2=0,$$

$${\ddot{x}}_2+{\omega }_2^2{x}_2+\frac{\beta }{m_g}{x}_2^3-\frac{\alpha }{m_g}{x}_1=0.$$

$${\mathcal{H}}_{eff}=\frac{p_1^2}{2m_g}+\frac{p_2^2}{2m_g}+\frac{1}{2}{m}_g{\omega }_1^2{x}_1^2+\frac{1}{2}{m}_g{\omega }_2^2{x}_2^2+\frac{1}{4}\beta (x_1^4+x_2^4)-\alpha x_1{x}_2.$$

$$H=\frac{p_1^2}{2m_g}+\frac{p_2^2}{2m_g}+\frac{1}{2}{m}_g{\left({\omega }_0-\frac{\Delta }{2}\right)}^2{x}_1^2+\frac{1}{2}h\cos (\left\{2{\omega }_0+{\Omega }_p\right\}t+{\chi }_1)x_1^2$$

$$+\frac{1}{2}{m}_g{\left({\omega }_0+\frac{\Delta }{2}\right)}^2{x}_2^2+\frac{1}{2}h\cos (\left\{2{\omega }_0+{\Omega }_p\right\}t+{\chi }_2)x_2^2+\beta (x_1^4+x_2^4)-\alpha x_1{x}_2.$$

$$\left[\begin{array}{cc}{\left({\omega }_0-\frac{\Delta }{2}\right)}^2& -\frac{\alpha }{m_g}\\ -\frac{\alpha }{m_g}& {\left({\omega }_0+\frac{\Delta }{2}\right)}^2\end{array}\right].$$

$${{\omega }^2}_{\pm }=\left({{\omega }_0}^2+{(\Delta /2)}^2\right)\left[1\pm \frac{{\omega }_0^2{\Delta }^2+(\alpha /m_g)}{\left({{\omega }_0}^2+{(\Delta /2)}^2\right)}\right],$$

$$\delta \omega =\sqrt{\frac{{\omega }_0^2{\Delta }^2+(\alpha /m_g)}{{{\omega }_0}^2+{(\Delta /2)}^2}},$$

$$\alpha \simeq \frac{m_g{\Delta }^2}{2},$$

$${\ddot{x}}_1+\gamma {\dot{x}}_1+{\omega }_1^2{x}_1+\frac{h}{m_g}\cos \left\{({\omega }_1+{\omega }_2)t+{\chi }_1\right\}{x}_1+\frac{\eta }{m_g}{x}_1^2{\dot{x}}_1+\frac{\beta }{m_g}{x}_1^3-\frac{\alpha }{m_g}{x}_2=f_1(t),$$

$${\ddot{x}}_2+\gamma {\dot{x}}_2+{\omega }_2^2{x}_2+\frac{h}{m_g}\cos \left\{({\omega }_1+{\omega }_2)t+{\chi }_2\right\}{x}_2+\frac{\eta }{m_g}{x}_2^2{\dot{x}}_2+\frac{\beta }{m_g}{x}_2^3-\frac{\alpha }{m_g}{x}_1=f_2(t).$$

$$\frac{dx}{dt}=q(x,t)+g(x,t)\zeta (t),$$

$$x(t+h)=x(t)+{\int }_t^{t+h}q(x(s),s)ds+{\int }_t^{t+h}g(x(s),s)\zeta (s)ds.$$

$$x(t+h)=x(t)+hq(x(t),t)+q^{\prime }(x(t),t){\int }_t^{t+h}(x(s)-x(t))ds+w_h(t)g(x(t),t)$$

$$+g^{\prime }(x(t),t){\int }_t^{t+h}(x(s)-x(t))\zeta (s)ds+(h+w_h(t))\times O[(x(s)-x(t){)}^2],$$

$$g^{\prime }(x(t),t){\int }_t^{t+h}(x(s)-x(t))\zeta (s)ds=\frac{1}{2}{g}^{\prime }(x(t),t)g(x(t),t)w_h^2(t)+O[h^{3/2}].$$

$$x(t+h)=x(t)+hq(x(t),t)+g(x(t),t)\sqrt{h}u(t)+\frac{h}{2}g(x(t),t)g^{\prime }(x(t),t)u^2(t).$$

$$x(t+h)=x(t)+hq(x(t),t)+D\sqrt{h}u(t).$$

$${\ddot{x}}_1+ϵ\gamma {\dot{x}}_1+{\left(\omega -ϵ\frac{\delta {\omega }_1}{2}\right)}^2{x}_1+ϵ\frac{h}{m_g}\cos (2\omega t+{\chi }_1)x_1+ϵ\frac{\eta }{m_g}{x}_1^2{\dot{x}}_1+ϵ\frac{\beta }{m_g}{x}_1^3-ϵ\frac{\alpha }{m_g}{x}_2=f_1(t)$$

$${\ddot{x}}_2+ϵ\gamma {\dot{x}}_2+{\left(\omega -ϵ\frac{\delta {\omega }_2}{2}\right)}^2{x}_2+ϵ\frac{h}{m_g}\cos (2\omega t+{\chi }_2)x_2+ϵ\frac{\eta }{m_g}{x}_2^2{\dot{x}}_2+ϵ\frac{\beta }{m_g}{x}_2^3-ϵ\frac{\alpha }{m_g}{x}_1=f_2(t),$$

$$x_1(t,T)=x_1^0(t,T)+ϵx_1^1(t,T)+\cdots$$

$$x_2(t,T)=x_2^0(t,T)+ϵx_2^1(t,T)+\cdots$$

$$\frac{d}{dt}=\frac{\partial }{\partial t}+ϵ\frac{\partial }{\partial T}=D_t+ϵD_T,$$

$$\frac{d^2}{d{t}^2}=\left(\frac{\partial }{\partial t}+ϵ\frac{\partial }{\partial T}\right)\left(\frac{\partial }{\partial t}+ϵ\frac{\partial }{\partial T}\right)=D_{tt}+2ϵD_T{D}_t+{ϵ}^2{D}_{TT}.$$

$$D_{tt}{x}_1^0+{\omega }^2{x}_1^0=0,$$

$$D_{tt}{x}_2^0+{\omega }^2{x}_2^0=0,$$

$$D_{tt}{x}_1^1+{\omega }^2{x}_1^1=-2(D_t{D}_T{x}_1^0)-\gamma (D_t{x}_1^0)+\omega \delta {\omega }_1{x}_1^0-\frac{h}{m_g}\cos (2\omega t+{\chi }_1)x_1^0-\frac{\eta }{m_g}{(x_1^0)}^2(D_t{x}_1^0)$$

$$-\frac{\beta }{m_g}{(x_1^0)}^3+\frac{\alpha }{m_g}{x}_2^0,$$

$$D_{tt}{x}_2^1+{\omega }^2{x}_2^1=-2(D_t{D}_T{x}_2^0)-\gamma (D_t{x}_2^0)+\omega \delta {\omega }_2{x}_2^0-\frac{h}{m_g}\cos (2\omega t+{\chi }_2)x_2^0-\frac{\eta }{m_g}{(x_2^0)}^2(D_t{x}_2^0)$$

$$-\frac{\beta }{m_g}{(x_2^0)}^3+\frac{\alpha }{m_g}{x}_1^0.$$

$$x_1^0=\left(\frac{a(T)}{2}\exp (i\omega t)+\text{c.c.}\right),$$

$$x_2^0=\left(\frac{b(T)}{2}\exp (i\omega t)+\text{c.c.}\right).$$

$$\frac{d{a}_r}{dT}=$$

Full text

††preprint: XXX

In “More is different”, Anderson speculated on the possibility of symmetry breaking and emergence of stable collective phases in time, in close analogy with emergent crystalline phases in space [1]. Over last decade, renewed interest in possible temporal symmetry breaking resulted in the discovery of discrete time crystals (DTCs), which are emergent phases of interacting periodically driven systems [2, 3, 4, 5, 6, 7, 8, 9, 10]. In a DTC phase, multiple interacting modes of a driven system settle into a common period that is longer than the drive. Such collective breaking of discrete time translational symmetry results in a dynamical steady state, the DTC phase, that exhibits novel features reminiscent of equilibrium thermodynamic phases. Such phases have been suggested for applications in computation, metrology, and sensing  [11, 12, 13]. Experiments in quantum systems with ions, atoms and defects in solids observed signatures of DTC in the form of subharmonic entrainment [14, 15, 16]. Nevertheless, multiple questions regarding the origin and nature of DTCs remain unanswered. First, what mechanisms ensure stability of DTC phases against drive heating and fluctuations [17, 18]? Second, can a truly many-body DTC phase arise in a fully classical system [17]? Third, are there DTC phases with complex symmetry types beyond subharmonic entrainment [17, 19, 20, 21, 22]?

Here, we answer some of these questions by observing a range of new DTC phases as well as transitions between them in a classical NEMS. The phases satisfy three accepted criteria of genuine DTCs  [19, 23, 17, 12]:

(1) Interacting many-mode system: We fabricate a device in which multiple mechanical modes of a Silicon Nitride (SiNx) resonator interact with two spatially separated graphene resonators. Moreover, exceptional elastic and mechanical non-linear properties of graphene [24] make the interaction tunable, resulting in long-range interaction between otherwise non-interacting many (on average 16) SiNx modes (see SI, section IV).

(2) Collective breaking of discrete symmetry in time: When driven parametrically at twice of the average frequency of the hybrid modes, universal collective modes emerge sharply above a threshold drive power, breaking discrete time translational symmetry of the drive. In classical dynamical systems, such discrete time translation symmetry breaking is ubiquitous, ranging from parametrically driven swings to Faraday waves [25]. Though in principle such systems are composed of many modes, the symmetry breaking oscillation can nevertheless be described as parametric driving of a single normal mode without any distinct many-mode features  [26, 27, 28, 29]. In contrast, symmetry breaking in weakly coupled many-mode systems at finite temperature and fluctuations, a hall-mark of DTC phases, is non-trivial  [12, 11, 30, 23]. Even with two weakly coupled parametric oscillators, common overlap of dynamical instability regions of the resonator modes has lead to DTC-like features [19]. Remarkably, with many interacting modes we observe enhanced rigidity along with a rich palette of distinct DTC phases that can be distinguished by symmetry types including sub-harmonic, anharmonic and a novel biharmonic phase [15, 14, 16, 31, 32, 33]. Furthermore, the spectra, integrated over millions of oscillation periods and measured at differing spatial locations on the resonator surfaces, show evidence for emergence of long-range order in time and space due to interactions between modes of spatially distinct graphene and SiNx oscillators. We note that these DTC phases cannot be fully described within a mean-field model, as would be the case for a one or few-body system.

(3) Rigidity of DTC phases: In contrast to sensitivity of conventional phase synchronization to external parameters [34, 35], the observed DTC phases are surprisingly stable with respect to a broad range of parameter fluctuations. The rigidity to externally added noise is beyond predictions of effective two-mode mean-field models [36, 33]. With increased noise DTCs eventually melt to a non-crystalline phase, marked by broadened spectra and sharp phase boundary, a recently conjectured telltale signature of genuine many-body DTCs [17].

Analysis of a mean-field model indicate that for DTCs, akin to crystals in space, additional physical constraints lead to new emergent phases and such constraints are dynamical in nature, resulting in novel DTC phases due to bifurcations of fixed points and limit cycles with changing mode coupling and frequency detuning [36, 33, 23].

Device and experiment

The device consists of a large area ($320\times 320$ $\mu$m2) SiNx resonator with atomically thin graphene deposited onto two $20$ $\mu$m holes drilled in its surface, forming a second set of resonators $G_{1}$ and $G_{2}$ (Fig. 1A). The system is driven by applying an AC gate voltage via gate electrode (see SI, section I). For initial characterization, thermal noise or low-power measurements are used. Vibrational spectra are measured with a confocal microscope focused on either of the graphene drums ($X_{1}$ or $X_{2}$ in Fig. 1A) or SiNx ($X_{3}$ in Fig. 1A) [37]. When measured on SiNx, a low inbuilt tension and high quality factor result in individually resolvable but densely packed spectra of modes (Fig. 1B, below). Measurements on graphene reveal lower quality factors fundamental modes of the relevant drums hybridized with multiple SiNx modes (Fig. 1B, above). This hybridization induces effective nonlinear cross-talk between the SiNx modes (Fig. 1B) [38]. With increasing tension in the system, controlled by the DC component of gate voltage, the frequencies of flexible graphene drums up-shift (dotted lines in Fig. 1C) while the modes of rigid SiNx remain virtually unchanged (horizontal lines in Fig. 1C). Nevertheless, G1 and G2 remain spectral decoupled in the entire range of DC gate voltages and interact with differing sets of SiNx modes. No signature of direct interaction (e.g. avoid crossing vs. DC gate voltage [38]) between the two graphene resonators is observed.

Emergent DTC phases

To access time-crystalline behaviours, the device is parametrically driven at twice of the mean frequency of the clusters of hybrid modes around G1 and G2 (inset, Fig. 2A). With increasing drive strength, the spectral amplitude at $\omega$ first rises from the noise floor at $4$V and then undergoes a sharp transition to a high-amplitude steady state at around $8$V (Fig. 1D). In region I of Fig. 2A, the corresponding spectra contain distinct signatures of graphene-SiNx interactions of the separated clusters of thermal modes around G1 and G2. With increasing drive voltage, the two clusters merge sharply (around 4V, region II) to a single shared spectra, that nevertheless remain thermal-like with small amplitude and without any sharp oscillation frequencies. This suggests that a collective many-body state comprised of multiple modes in remote resonators develops satisfying criterion 1. Eventually, at around 8V, a sharp subharmonic peak, satisfying criterion 2, emerges at frequency $\omega$ with oscillation amplitude that are orders of magnitude higher compared to the thermal modes (region III). The DTC spectra is observed to be global in nature, with similar features measured on resonator surfaces at $X_{1}$, $X_{2}$ (blue and green respectively, in Fig. 2B) and $X_{3}$, in contrast to the local features observed below the threshold (Fig. 1B). Along with long-range order in space, the spectra remain stable for millions of oscillation periods over time scales of seconds.

Tuning the frequency separation of graphene and sets of SiNx with DC gate voltage, we observe two other distinct phases with similar long-range order properties in addition to the subharmonic phase shown in Fig. 2B. These are anharmonic (Fig. 2C) and biharmonic (Fig. 2D) phases, distinguished by even or odd numbers of frequency comb lines in the spectra. The subharmonic, anharmonic, and biharmonic phases manifest different types of time translational symmetry breaking, satisfying the criterion 2 of DTC. Indeed, several recent theoretical works predict the first two of these DTC phases [32, 39, 33]. To the best of our knowledge, the  biharmonic phase is novel. In the remainder of the manuscript, we confirm that assignment of these as time-crystalline phases, map their phase diagram, and explore transitions between them.

It can be noted that the anharmonic DTC phase provides an indirect means of quantifying correlations by noting that the two sharp spectral peaks at frequencies $\omega_{1}$ and $\omega_{2}$ are pinned around the thermal peak clusters corresponding to the two graphene drums G1 and G2 (Fig. 2C). The centers of these drums (points $X_{1}$ and $X_{2}$) are separated by 40 micrometers. Therefore, by comparing the correlation between the peaks at $\omega_{1}$ and $\omega_{2}$, we can gauge the extent of spatial correlations in the system. For a single mode corresponding to either G1 or G2, quadrature fluctuations, measured in frame co-evolving at either of the two peak frequencies, show uncorrelated noise reminiscent of thermal equilibrium (Figs. 2E, F). In striking contrast, cross-correlation measurements between the quadrature of two different modes (effectively probing inter-drum correlations) show strong correlation, with a $24$ $dB$ suppression along a specific phase compared to fluctuations in the orthogonal quadrature (Figs. 2G, H), providing additional evidence of emergent long-range order in higher-order correlations or fluctuations. As expected, this correlation in fluctuations is absent below the DTC threshold and thereby, can be ascribed to it’s emergent long-range order.

A mean-field model

Mean-field models of few coupled modes have been successful in providing insights into dynamical features of DTCs [19, 33]. To investigate the nature of the observed phases further, we thereby model the system as two coupled nonlinear graphene modes (with transverse displacement $x_{1}$ and $x_{2}$) with a common parametric drive and open to a thermal bath with damping and fluctuations. The coarse graining model replaces the zoo of many coupled SiNx modes with a direct coupling between the two otherwise decoupled graphene drums. In effect, the model is represented by two damped Mathieu-Van der Pol-Duffing oscillators with bi-linear coupling between them (see SI, section III and  [38]).

Numerical simulations (see SI, section III) indicate the same three DTC phases (subharmonic, anharmonic, and biharmonic; Figs. 3A, B, C) emerging in different parts of the phase diagram (Fig. 3D). This hints that few body classical models at least capture spectral features of genuine many-body DTC phases. A linear stability analysis of the model provides insights into the dynamical nature of the DTC phases. Collective oscillatory DTC phase are mathematically represented as solutions of sinusoidal with frequency $\omega$. Amplitudes of the solutions are slowly varying owing to the weak nonlinearity assumed to facilitate the corresponding rigorous perturbative analyses (see SI, section III). The phases then correspond to robust stable attracting solutions observed for a set of initial conditions of finite measure, while difference in DTC phases manifest as different asymptotic forms of the slowly varying amplitudes. For high enough detuning, stable fixed point, indicating non-crystalline phase (red region, Fig. 3D) appears at the origin (see inset). Reduction in detuning renders this stable solution unstable, eventually yielding limit cycle solutions symmetric about the origin (see inset) via the Hopf bifurcation (black curve, Fig. 3D). The corresponding steady-state oscillations have frequency shifted from $\omega$ resulting in anharmonic DTC phases (blue region, Fig. 3D) [32]. Remarkably, as the coupling strength is further decreased, the limit cycle too loses its stability. However, it bifurcates to give rise to two new stable limit cycles that are not symmetric about origin (see inset), leading to a novel unexplored scenario where DTC spectra bears both the spectral features at frequency $\omega$ as well as the frequencies symmetrically shifted from $\omega$. This suggest a new DTC phase: bi-harmonic (green region, Fig. 3D). The green region indicates existence of more complex solutions that need further detailed analysis. As the coupling strength is further reduced, the two limit cycles eventually disappear leaving only stable fixed points with constant amplitude, shifted away from the origin (see inset). This case corresponds to subharmonic entrainment (yellow region, Fig. 3D) with frequency at $\omega$  [15, 14, 16, 31]. Higher order effects of nonlinearities further complicate the spectra through a series side-peaks in the frequency spectra. Nevertheless, the mean-field model provides crucial insights into the dynamical nature of the DTC phases and dependence to a coarse-grained coupling and mode detuning.

Rigidity of DTC phases

Finally, we check the rigidity of the observed phases with respect to noise and parameter fluctuation, the final criterion of the true DTC phase. As an example, to test rigidity against fluctuations, we mix a white noise signal to the drive channel and record the resulting DTC spectra with increasing noise amplitude. Fig. 4A shows a plot of the peak of the subharmonic spectra that remains stable till the relative noise becomes comparable to the drive and eventually melts to a broad noisy spectrum (Fig. 4B). Remarkably, this transition is marked by a sharp phase boundary at specific noise threshold, reminiscent of a first-order thermodynamic phase transition (Fig. 4C, blue). In contrast, numerical simulation of mean field model shows rapid fluctuation between the DTC and a melted phase above a certain noise strength. When averaged over many runs, no sharp transition is observed (Fig. 4C, black). We interpret the stark difference between the experimental system and a mean-field model as a direct evidence for the many body nature of the emergent DTC phase.

The phases are stable in a large frequency window when the drive frequency is detuned. As a function of drive frequency (and at constant tension), we observe broad bands of DTC phases with sharp transitions between them (Fig. 4D). In a remarkably rich phase diagram, we observe additional phases, distinguished by their symmetry, in addition to previously seen subharmonic, anharmonic and biharmonic phases. The continuum phase characterized by a broadband spectrum, requires further investigation and can be a melted, chaotic or a coherent phase. The non-crystalline phase appears for certain pump detuning values which indicates that the drive strength is below threshold. This reflects the dependency of parametric threshold value on pump detuning and coupling of modes. Overall the phase diagram of our hybrid NEMS resonators with transitions between distinct DTC phases rivals the complexity of crystallographic phases of most solid state system.

Outlook

Our results establish existences of truly many-body DTCs in classical systems with distinctly different phases. Observation of rigidity against parameter variations and sharp phase boundaries are found to be emergent many-body features beyond mean-field paradigm. Furthermore, the hybrid NEMS, with exceptional nonlinear and elastic properties of graphene resonators are found to be essential in providing tunable interactions, nonlinearity and damping between many SiNx modes resulting in the stability of the DTC phases.

These results imply multiple directions to explore. Unexplained spectral features (Region II of Fig. 2A and continuum phase of Fig. 4D), can lead to understanding of liquid-like phases in time [21, 9, 8] and the nature of eventual thermalization in these non-equilibrium systems [21, 20]. Correspondence of these DTC phases with their quantum counterparts [20, 21] can lead to insights into interacting quantum systems. For NEMS, multiple graphene resonators in spatially ordered geometries can be extended to new space-time crystalline phases, exceptional spectral rigidity of DTCs can be applied as unique frequency markers in information technology [40, 41, 42] while tunability to phase boundaries can be used in sensors and metrology [43, 44, 45, 46, 47].

Acknowledgments

We thank Adhip Agarwala, Souvik Bandyopadhyay, Sayak Bhattacharjee, Diptarka Das, Amit Dutta, Alexander Eichler, Shilpi Gupta and V. Prasanna for illuminating discussions. This work is supported by SERB, Department of Science and Technology, India under DST grant no: SERB/PHY/2015404, ERC Grant No. 639739 and DFG, German Research Foundation, projects 449506295 and 328545488, CRC/TRR 227. A.S acknowledges CSIR, New Delhi and A.K.R acknowledges MHRD for financial support. J.A.M. acknowledges Prime Minister’s Research Fellows (PMRF) scheme of the Ministry of Human Resource Development, Govt. of India for financial support.

I Materials and Methods

I.1 Device

Low stress silicon-rich silicon nitride is deposited on both sides of a $300\times 300$ $(\mu m)^{2}$ area silicon membrane chip. Array of through holes of $10$, $15$ and $20$ $\mu m$ is patterned on the $300$ $nm$ thick silicon nitride surface. Metallic contact is deposited on the top surface for electrical actuation. High quality chemical vapor deposition (CVD) growth and wet transfer technique is used to make mono layers. Due to Van der Waals force, graphene membranes are clamped on the through holes as a suspended circular drum. The sample is subsequently annealed in an $Ar-H_{2}$ environment at $350^{o}C$.

I.2 Experimental Setup

Sample is placed in a capacitive arrangement, with a highly doped silicon surface used as a back electrode (see Fig. S.1). DC gate voltage is used to tune the graphene mode frequencies, while AC gate voltage actuates hybrid vibrational modes. The device is kept at room temperature and in vacuum (pressure $\sim 10$ mTorr ). Out of plane motion of the membranes are detected using a confocal microscope, which forms one arm of a Michelson interferometer. To stabilize the signal against external drifts and fluctuations, the two path lengths are actively locked. A frequency and amplitude stabilized external cavity diode laser (ECDL) is used to probe the system. Power of the input signal is kept very low ($\sim 400$ $\mu W$) such that it does not heat up the graphene sample. Detected signal is analyzed with a DAQ for imaging, spectrum analyzer to observe frequency response and a lock-in amplifier for homodyne detection  [37, 38].

II Evidence of correlated fluctuations and long-range order in anharmonic phase

To measure correlation between two spatially separated points on the surface, we study anharmonic response more closely. We have understood earlier that two different peaks of anharmonic response correspond to two different measurement points, one is on graphene drum 1 ($X_{1}$) and the other one is on drum 2 ($X_{2}$). Quadrature measurement using a network analyzer is done to detect fluctuations at $\omega_{1}$ and $\omega_{2}$. We define in phase and out of phase components of peak $\omega_{1}$ and $\omega_{2}$ as $X_{1},Y_{1},X_{2}$ and $Y_{2}$, respectively. When system is not driven, we observe thermal noise distribution (see Fig. S.2A and B ) in phase space due to room temperature.

Furthermore, to quantify the correlation between modes, absolute correlation coefficient is calculated for $X_{1},Y_{1},X_{2},Y_{2}$, defined as [48]

[TABLE]

Here, $cov(Z_{i},Z_{j}(\tau))=\sum_{k}(Z_{i,k}-\bar{Z_{i}})(Z_{j,k}(\tau)-\bar{Z_{j}})/N$, $N$ is total number of data points, $Z_{i}\in\{X_{1},Y_{1},X_{2},Y_{2}\}$ , $\sigma_{Z_{i}}$ is the variance and $\tau$ is time delay between two data points.

Fig. S.2C is showing absolute correlation coefficients at $\tau=0$ correspond to thermal peak. We observe that only auto correlation terms are nonzero because when system is not driven, correlation between quadrature does not exist. When parametric pumping is on, overall fluctuation increases. Looking at $X_{1}$ vs $Y_{1}$ Fig. (S.2D) and $X_{2}$ vs $Y_{2}$ (S.2E) one can conclude that these are same as quadrature of a thermal peak but cross quadrature show anti-correlation (Fig. 2G and H in main text). As a result, we observe non-zero off-diagonal elements of covariance matrix (Fig. S.2F).

Fig. S.3 shows linear variation of correlation coefficients with increasing pump amplitude but with different slope. Covariance $X_{1}-X_{1}$ (red) remains one for all the time when pump strength is increased on the other hand $X_{1}-Y_{1}$ (black) remains zero. $X_{1}-X_{2}$ (navy blue) and $X_{1}-Y_{2}$ (dark cyan) increases with pump amplitude linearly. Navy blue curve saturates around $0.8$ and starts decreasing again.

III Details of the mean-field model

III.1 Effective Hamiltonian of the system

Fig. 1B in main text shows that the fundamental mode of graphene drums can have resonance coupling with many densely packed high quality factor higher order $SiNx$ modes. Here, we have two graphene drums, and both of them are coupled to a few substrate modes on resonance but they are not coupled to each other directly. In a similar fashion, $SiNx$ modes are not directly coupled to each other but coupling exists through graphene modes.

Hamiltonian of this coupled system can be represented as

[TABLE]

We can write following differential equations for this system:

[TABLE]

Now, we know that $SiNx$ is much heavier than graphene membrane $(m_{s}\gg m_{g})$. So displacement of graphene drums are order of magnitude larger than $SiNx$. Observing thermal peaks at fundamental frequencies we take ansatz $x_{1}\sim\cos(\omega_{1}t)$ and $x_{2}\sim\cos(\omega_{2}t)$. In this situation, we will focus on $SiNx$ equations. Because secular terms decay down very fast, only coupling terms act as driving forces for this equation, implying $x_{i}\sim A_{i}\cos(\omega_{1}t)+B_{i}\cos(\omega_{2}t)\sim A^{\prime}_{i}x_{1}+B^{\prime}_{i}x_{2}$ ($A_{i}$, $B_{i}$, $A^{\prime}_{i}$, and $B^{\prime}_{i}$ are some slow variables of time) $\quad\forall i=\{3,4,\cdots,n\}$.

If we plug in this into equation S.3 a and b, we can rewrite equations,

[TABLE]

One of the two coupling terms will get absorbed as fundamental frequency correction terms, and finally we can think of an effective coupling between two graphene drums.

We write effective Hamiltonian for the system,

[TABLE]

In our experiments, we take $\omega_{1}=\omega_{0}-\Delta/2$ and $\omega_{2}=\omega_{0}+\Delta/2$, and parametrically drive the system with frequency $2\omega_{0}+\Omega_{p}$ such that the Hamiltonian becomes:

[TABLE]

Estimation of coupling value for simulation: Estimation of this effective coupling value will guide us to do numerical calculations. For $h=0$, we will try to solve Eq. (S.4) by solving eigenvalue problem of the following matrix:

[TABLE]

Eigenvalues of this matrix are,

[TABLE]

or,

[TABLE]

where $\delta\omega=\omega_{+}-\omega_{-}$.

Due to small coupling $\alpha$, we assume that the separation between the peaks at $\omega_{+}$ and $\omega_{-}$ is approximately equal to the detuning $\Delta$, we arrive at

[TABLE]

for $m_{g}\approx{10}^{-15}$ and $\Delta\approx{10}^{5}$ we get $\alpha\approx{10}^{-4}$.

Fine tuning of $\alpha$ can be done by looking at experimental data.

III.2 Numerical integration of stochastic differential equation

As this system is in contact with a thermal bath at the room temperature, we can introduce effective damping and thermal noise $f$ and write coupled second order differential equations:

[TABLE]

Any standard ordinary differential equation (ODE) solver will not give correct result due to stochastic terms. Also while checking robustness, another stochastic term will be added to the equation. So we decide to follow Milstein’s algorithm to solve the equations numerically.

III.2.1 Algorithm to solve stochastic differential equation (SDE)

Let us take a first order differential equation with deterministic and stochastic term  [49]:

[TABLE]

where $q(x,t)$ and $g(x,t)$ are differentiable functions and $\zeta(t)$ is Gaussian white noise ($<\zeta(t)>=0$ and

$<\zeta(t_{1})\zeta(t_{2})>=\delta(t_{1}-t_{2})$).

For $g(x)=0$, this equation becomes an ODE and can be solved using Euler algorithm or RK4.

After a discrete time interval $h$, updated value of $x$ can be expressed as

[TABLE]

Taylor series expansion of functions $q(x,t)$ and $g(x,t)$ give

[TABLE]

where $w_{h}(t)=\int_{t}^{t+h}\zeta(s)ds$.

Incorporating properties of Gaussian noise one can write $w_{h}(t)=\sqrt{h}u(t)$, where $u(t)$ is a set of independent Gaussian random variable with zero mean and variance $1$.

We want to consider all the terms up to first order of h. We can ignore third term of the equation because it is of the order of $h^{2}$. As $w_{h}(t)$ depends on $\sqrt{h}$, we need to evaluate integral of the fourth term. Using Stratonovich calculus we get,

[TABLE]

So equation III.2.1 becomes

[TABLE]

This recurrence relation is known as the Milshtein algorithm. To use this algorithm for additive noise one have to consider $g(x(t),t)$ is a constant. In that case this relation becomes Euler-Maruyama algorithm:

[TABLE]

Here D is a constant represents additive noise strength.

In equation S.11, thermal fluctuation term is additive, so Euler-Maruyama algorithm will work for us.

III.2.2 Numerical results for DTC phases

We solve equations S.11 using the above mentioned algorithm and we get different spectra for different set of parameter values. Comparing with experimental data, one may conclude that our mean field model is sufficient to explain the time crystalline phases but this results are lacking of essential many body features of the system.

Subharmonic: In experiment, at higher DC gate voltage fundamental modes of graphene drums are close enough to show subharmonic phase when driven parametrically. All the contributing modes start oscillating at half of the pump frequency (Fig. 2A, main text). Above a threshold value ($8V$) of pump voltage we observe subharmonic response. Using equations S.11 with suitable parameters we simulate similar behavior (see Fig. S.4).

Anharmonic: We observe transition from noncrystalline phase to anharmonic phase by changing pump voltage when system is driven at $\omega_{1}+\omega_{2}$. Frequency spectrum (see Fig. S.5A) on $X_{1}$ is recorded. We simulate similar result (Fig. S.5B) using coupled nonlinear equations (S.11).

Biharmonic: We observe transition from noncrystalline phase to anharmonic phase and anharmonic phase to biharmonic phase by changing pump voltage. First transition happens at threshold value $\sim 7V$ and threshold value for biharmonic is $9.5V$. We record frequency spectrum (see Fig. S.6A) and simulate similar result (Fig. S.6B) using coupled nonlinear equations (S.11).

Continuum: Apart from three different time crystalline phases we observe another phase above the parametric pump threshold value which is not discussed in details in this work. For certain parameters, we observe amplification of a wide spectrum around half of the pump frequency (see Fig. S.7A) which we identify as continuum. We observe same feature while simulating our coupled equation for very specific range of parameters (Fig. S.7B). We are still investigating this phase to understand in details.

We study these time crystalline phases experimentally by tuning pump frequency (Fig. 4D in main text). A numerical plot (Fig. S.7C) with two coupled graphene mode shows few features same as Fig. 4D (main text) as also suggestive from Fig. 3 (main text) but experimental data is much more richer than Fig. S.7C due to many-body nature of the system, variation of number of modes, coupling and detuning in experiment.

III.3 Perturbative analysis for a mean-field phase diagram

To understand analytically, we introduce a small book-keeping parameter $\epsilon\equiv 1$ to keep track the order of perturbative equation. We define, $\omega_{0}\equiv\omega-\epsilon\Omega_{p}/2$. Now the model equations take the form,

[TABLE]

where, $\delta\omega_{1}=\Omega_{p}+\Delta$ and $\delta\omega_{2}=\Omega_{p}-\Delta$.

We are now interested in the robust solutions of the above equations. The system is analytically tractable only through perturbative analysis. Considering only the deterministic part, we adopt the multiple time scales method (any equivalent perturbative techniques could be used) and start by writing the ansatz for the equation (S.18a) and equation (S.18b) as,

[TABLE]

As here we have two time scales $t$ and $T\equiv\epsilon t$, therefore

[TABLE]

Using the ansatz given by equations (S.19a)-(S.19b) in equations (S.18a)-(S.18b), the $0^{\textrm{th}}$ order and the $1^{\textrm{st}}$ order equations in $\epsilon$ are,

[TABLE]

The solution of the $0^{\textrm{th}}$ order equations for $x_{1}$ and $x_{2}$ are,

[TABLE]

Using the zeroth order solutions in Eq. (S.20) and Eq. (S.20) and equating the coefficients of $\exp(i\omega t)$ equal to zero, to remove the secular terms, in differential equations for $x_{1}$ and $x_{2}$, we get the evolution of the real and the imaginary parts of the amplitudes, given by,

[TABLE]

where, $a_{r}$ and $b_{r}$ are real parts of $a$ and $b$ respectively. $a_{i}$ and $b_{i}$ are the imaginary parts of $a$ and $b$ respectively.

Till now, we have discussed the theoretical model and corresponding equations of motion. Using the multiple timescales analysis, we find the first order differential equations for real and imaginary parts of amplitudes in equations (S.22d)-(S.22d). We are now all set to study the dynamics of the system in four dimensional phase space via linear stability analysis. In Fig. S.8, we show various dynamical phases of the system in two dimensional planes at different sets of parameter values.

Non crystalline phase: In Fig. S.8A and Fig. S.8B, we show the fixed points in $a_{r}$-$b_{r}$ space corresponding to the non crystalline phase of the system. In both the cases, origin in is the stable fixed point. As in Fig. S.8A, there exists only one fixed point at origin which is stable, therefore any initial condition will end up at this point and there will be no amplification and it corresponds to the non crystalline phase. In Fig. S.8B also, origin is stable fixed point and the other fixed point are far away from it therefore in small amplitude regime it corresponds to the noncrystalline phase.

Subharmonic: In Fig. S.8C and Fig. S.8D, we depict the fixed points in $a_{r}$-$b_{r}$ space corresponding to the subharmonic phase of the system. In this case, origin becomes unstable fixed point along with other stable fixed points. Fixed point diagram for subharmonic phase in Fig. S.8C have four other fixed points (other than zero) in which two (in blue) correspond to the stable and two (in black) correspond to the unstable one. However, the subharmonic phase in Fig. S.8D have eight other fixed points (other than zero) in which four (in blue) corresponds to the stable and four (in black) corresponds to the unstable one.

Anharmonic: Fig. S.8E and Fig. S.8F depict the fixed points in $a_{r}$-$b_{r}$ space corresponding to the anharmonic phase of the system. In the anharmonic phase, the unstable origin is enclosed by a limit cycle. Limit cycle is symmetrically placed around the origin. Fixed point diagram for anharmonic phase in Fig. S.8E have only one unstable fixed point (origin) and a limit cycle enclosing it. Whereas fixed point diagram for anharmonic phase in Fig. S.8F not only have one unstable fixed point (origin) and a limit cycle enclosing it but also there are four other fixed points in which two (in black) are unstable and other two (in blue) are stable.

Biharmonic: Fig. S.8G represents the fixed points diagram in $a_{r}$-$b_{r}$ space corresponding to the biharmonic phase of the system. In the anharmonic phase, the unstable origin is enclosed by two limit cycles. In this case, limit cycles are not symmetrically placed around the origin.

The biharmonic phase can be characterized in two ways. Firstly, in the biharmonic phase one finds the appearance of a peak in the spectrum at the central frequency $\omega$ (Fig. S.6). Secondly, in the quadrature space, as we tune the parameters to move from the anharmonic to biharmonic phase, one finds that the corresponding limit cycle breaks the reflection symmetry about the origin. Consequently, there are two limit cycle solutions to the perturbative equation connected by the transformation $(a_{r},a_{i},b_{r},b_{i})\xrightarrow{}(-a_{r},-a_{i},-b_{r},-b_{i})$ (See Fig. S.8G). This also has the implication that mean position of the limit cycle is no longer at the origin of the quadrature space. In order to see the connection between two characterizations, consider one of the displacements, say $x(t)$. To keep the argument clear, we indulge in abuse of notation and ignore complex conjugation terms to write $x(t)=a(t)e^{i\omega t}$. Now, since for a limit cycle $a(t)$ is periodic, one may write it as a Fourier series,

[TABLE]

Thus the displacement takes a form

[TABLE]

This equation says that if the quadrature space limit cycle contains a frequency component $\omega_{i}$, there is a peak in the spectrum of $x$ at $\omega+\omega_{i}$. As mentioned above, in the biharmonic phase limit cycle is not symmetric about the origin and thus at least one component in the quadrature space has a nonzero time average. A nonzero time average in turn implies a zero frequency mode, or $a_{0}\neq 0$. Thus $a$ has an $\omega_{i}=0$ mode which in turn means that $x$ has a peak at $\omega+0=\omega$.

We remark here that the breaking of the reflection or $\mathbb{Z}_{2}$ symmetry by the limit cycle suggests a connection to Landau’s theory of ferromagnetism. Particularly, the distance of the limit cycle center from the origin functions effectively as an order parameter for the transition. Moreover, the direct forcing term that we have ignored in the analysis may function as the analogue of the conjugate field coupling to the order parameter. We leave this analogy here for study in a future work.

IV Quintessential many-body nature of phases: fluctuation and rigidity

An interacting many-body system shows different time crystalline phases depending on coupling parameter, drive strength and fluctuation present in the system. Due to high mode density, graphene peak interacts with multiple $SiNx$ modes and number depends on the parameters of the experiment.

IV.1 Mode density of $SiNx$ resonator

We do a rough estimation of $SiNx$ mode density of the sample we are working on. Using a lock-in amplifier we send a signal to drive the system weakly. Tuning frequency around the experiment range ($2$ $MHz$ to $4$ $MHz$) we estimate average mode density $88$ peaks/$MHz$. Now we can roughly estimate number of modes showing different time crystalline phases. We will consider modes having frequency in between two graphene drums.

For subharmonic response, frequency difference between graphene peaks are $\sim 100$ $kHz$ and full width of half maxima (FWHM) of a graphene peak is $34.42$ $kHz$. So we have $\sim 12$ $SiNx$ modes in the frequency span of $134.42$ $KHz$. For anharmonic and biharmonic case, difference between graphene peaks are $150$ $KHz$ and $125$ $kHz$ respectively. So, we have $\sim 16$ and $\sim 14$ $SiNx$ modes close to graphene drums when anharmonic and biharmonic phases emerged.

IV.2 Nature of added noise

To test rigidity of time crystalline phases we feed in noise across the capacitive arrangement and record the spectrum with increasing r.m.s. amplitude of the noise ($V_{r.m.s.}$). Experimentally noise is added using a function generator. To understand the nature of noise, we look at Fourier transform of added noise. Blue curve in Fig. S.10A is instrumental noise floor without any added noise for a large range of frequency (0 to 10 MHz). Added noise gives a flat spectrum (red curve in Fig. S.10A) for the same range.

For numerical analysis, we choose a random function to add noise to the system. Amplitude of the function is Gaussian for uniform time step. Fourier transform of this function also gives a flat spectrum. So we will use this function for further numerical analysis.

IV.3 Additive and multiplicative term from capacitive arrangement

Our system is kept in a capacitive arrangement and force is applied through metallic contacts of sample and back gate.

For simplicity of calculation, we rewrite equation S.4a for $\alpha=0$, adding linear damping ($\gamma$) and a force term $F(x,t)$:

[TABLE]

Now, we will calculate force $F(x,t)$ acting on the sample due to applied voltage $V(t)$ to capacitor.

From [50] we adopt the equation for stored elastic energy of the system,

[TABLE]

Here, $E,A,\chi_{0}$ and $L$ represent young modulus, cross sectional area, built in strain and diameter of the suspended material respectively.

So, total force due to this term is

[TABLE]

For a parallel plate capacitor, we can write

[TABLE]

where, $\epsilon_{0}$ and $d_{0}$ represent permittivity of air and mean distance between parallel plates respectively.

Stored electrostatic energy can be expressed as

[TABLE]

where, $V(t)$ is applied voltage across capacitor plates.

So, force due to this term is

[TABLE]

We expand $F_{es}$ around $x_{0}$ up to first order [51];

[TABLE]

So total force acting on the sample perpendicular to the surface is

[TABLE]

If voltage varies over time ($V(t)=V_{DC}+V_{AC}+V_{noise}$), we can rewrite equation S.25:

[TABLE]

A time varying voltage across capacitive arrangement have two contribution to the equation of motion. One term modifies the frequency and other term acts as direct drive term.

IV.4 Noise axis calibration

To compare experiment with numerical result, we need to calibrate noise strength axis. Experimentally we have observed that noise floor increases with increasing added noise strength. Similar things happens in numerical analysis if we consider additive noise term (see Fig. S.11A). We have seen earlier that capacitive force has a leading order term which acts as a force term. Lorentzian function is fitted on average spectrum of fifty iteration for different noise strength (see Fig. S.11A). Noise floor offset values are extracted as a fitting parameter and plotted against noise strength (see Fig. S.11B).

Next term of capacitive force is multiplicative. To understand higher order contribution we turn off additive term in numerical analysis and introduce multiplicative noise term. One must remember that in the case of multiplicative noise Milstein algorithm should be used instead of Euler–Maruyama algorithm. Here again lorentzian function is fitted on average spectrum of fifty iteration for different noise strength (see Fig. S.11C). Unlike experiment we observe non-linear frequency shift with increasing noise strength (see Fig. S.11D). Also note that here noise floor is remaining constant with increasing noise strength. So we can conclude that only leading order term of capacitive forcing is playing role in experiment when noise is added.

Now we can define a new parameter, relative noise spectral density, by taking ratio of extracted noise floor and maxima of thermal peak of the mode. As this is now unit-less quantity, this axis can be used as control parameter to observe phase transition in experiment and numerical analysis.

V Robustness to fluctuations

We focus on subharmonic phase to study robustness of the phase against added noise. In experiment the phase remains stable for certain threshold value of relative noise spectral density (Fig. 4A and B in main text). In numerical analysis, we observe similar nature of spectra (Fig. S.13) but do not observe sharp transition in the two-mode model (Fig. 4C in main text).