††thanks: Equal contribution††thanks: Equal contribution

# Observation of a critical prethermal discrete time crystal created by two-
frequency driving

William Beatrez  Department of Chemistry, University of California, Berkeley,
Berkeley, CA 94720, USA.  Christoph Fleckenstein  Department of Physics, KTH
Royal Institute of Technology, SE-106 91 Stockholm, Sweden.  Arjun Pillai
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Erica Sanchez  Department of Chemistry, University of California,
Berkeley, Berkeley, CA 94720, USA.  Amala Akkiraju  Department of Chemistry,
University of California, Berkeley, Berkeley, CA 94720, USA.  Jesus Alcala
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Sophie Conti  Department of Chemistry, University of California,
Berkeley, Berkeley, CA 94720, USA.  Paul Reshetikhin  Department of Chemistry,
University of California, Berkeley, Berkeley, CA 94720, USA.  Emanuel Druga
Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.  Marin Bukov<EMAIL_ADDRESS>Department of Physics, St.
Kliment Ohridski University of Sofia, 5 James Bourchier Blvd, 1164 Sofia,
Bulgaria.  Max Planck Institute for the Physics of Complex Systems, Nöthnitzer
Str. 38, 01187 Dresden, Germany.  Ashok Ajoy<EMAIL_ADDRESS>Department
of Chemistry, University of California, Berkeley, Berkeley, CA 94720, USA.
Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley,
CA 94720, USA.

###### Abstract

We report the observation of long-lived Floquet prethermal discrete time
crystalline (PDTC) order in a three-dimensional position-disordered lattice of
interacting dipolar-coupled ${}^{13}\textrm{C}$ nuclei in diamond at room
temperature. We demonstrate a novel strategy of ”two-frequency” driving,
involving an interleaved application of slow and fast drives that
simultaneously prethermalize the spins with an emergent quasi-conserved
magnetization along the $\hat{\textbf{x}}$-axis, while enabling continuous and
highly resolved observation of their dynamic evolution when periodically
kicked away from $\hat{\textbf{x}}$. The PDTC order manifests itself in a
robust period doubling response of this drive-induced quasi-conserved spin
magnetization interchanging between $\hat{\textbf{x}}$ and
-$\hat{\textbf{x}}$; our experiments allow a unique means to study the
formation and melting of PDTC order. We obtain movies of the time-crystalline
response with a clarity and throughput orders of magnitude greater than
previous experiments. Parametric control over the drive frequencies allows us
to reach PDTC lifetimes up to 396 Floquet cycles which we measure in a single-
shot experiment. Such rapid measurement enables detailed characterization of
the entire PDTC phase diagram, rigidity and lifetime, informing on the role of
prethermalization towards stabilizing the DTC response. The two-frequency
drive approach represents the simplest generalization of DTCs to multi-
frequency drives; it expands the toolkit for realizing and investigating long-
lived non-equilibrium phases of matter stabilized by emergent quasi-
conservation laws.

##  I Introduction

Periodically driven (Floquet) quantum systems enable the realization of novel
phases of matter far from thermal equilibrium. An interesting example are
discrete time crystals (DTC) [1, 2, 3] – infinitely long-lived non-equilibrium
phases of matter characterized by the dynamical breaking of time-translation
symmetry [4, 5, 6, 7, 8, 9, 10], and stabilized by many-body localization,
signatures of which were observed in recent experiments [11, 12, 13, 14, 15,
16, 17]. In the absence of localization and at high drive frequencies, energy
absorption in Floquet systems is suppressed to exceedingly long times – a
phenomenon known as “Floquet prethermalization” [18, 19, 20], wherein the
periodically driven many-body state has a lifetime ($T_{2}’$) that is
significantly enhanced with respect to the natural system free induction decay
time $T_{2}^{\ast}$. The heating rate $\propto\\!(T_{2}^{\prime})^{-1}$ is
exponentially suppressed for sufficiently local interactions and large driving
frequencies $\omega_{f}{\gg}J$, where $J$ is the intrinsic spin coupling
strength [21, 22, 23]. Recently, it has been suggested that Floquet drives can
yield prethermal discrete time crystals (PDTCs) [24], non-equilibrium
metastable states characterized by a robust subharmonic response in the drive
frequency, and with a parametrically long lifetime that harnesses
prethermalization. This has also inspired extensions of time crystal phenomena
to classical systems [25, 26, 27]. Nevertheless, characterizing the full phase
diagram of the emergent prethermal DTC order, rigidity, and elucidating its
thermalization dynamics towards infinite temperature, remains a challenging
task [28, 29]. At the same time, this is critical for advancing our
fundamental understanding of non-equilibrium order, and for leveraging such
collective phenomena in applications, such as in quantum simulation and
sensing.

Fig. 1: Experimental Implementation. (A) System. Dipolar lattice of
${}^{13}\textrm{C}$ nuclei in diamond. Dashed lines display representative
interspin couplings. Optically pumped NV centers hyperpolarize the
${}^{13}\textrm{C}$ nuclei (blue arrows). (B) Concatenated Floquet drive at 7T
consists of interleaved application of a fast-drive (blue) along
$\hat{\textbf{x}}$ and slow kicks (yellow) along $\hat{\textbf{y}}$ separated
by period $T$. Fast drive is composed of a train of $\vartheta({\neq}\pi)$
pulses (green) separated by period $\tau$ that spin-lock the nuclei along
$\hat{\textbf{x}}$. Kick interval $T{=}N\tau$, and the ratio $N$ is tunable
within the range 1-$10^{3}$. System undergoes evolution under interspin
couplings in intervals between the pulses. Spins are interrogated by RF
induction in $t_{\textrm{acq}}{=}64\mu$s windows between the pulses (red zig-
zag lines); Larmor precession is sampled every 1ns. This effectively yields
the ability to track projections along $\hat{\textbf{x}}$ and
$\hat{\textbf{y}}$ quasi-continuously with period $\tau$. (C) Schematic Bloch
sphere depiction of the reduced density matrix of a single spin in the
rotating frame. Following a kick, spins prethermalize towards the quasi-
stationary state along $\hat{\textbf{x}}$. Continuous tracking of the spins in
$\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane during this process is
highlighted. Kicks shown are for the example $\gamma{=}\pi+\varepsilon$, where
a PDTC state manifests as a periodic switching of the spins between
$+\hat{\textbf{x}}$ and $-\hat{\textbf{x}}$ between successive kicks.

In this paper, we report on a novel experimental approach for the high-
throughput characterization of the formation and melting of PDTCs, which
permits measuring their phase diagram with unprecedented high resolution. Our
experiments are performed on an interacting system of driven, hyperpolarized,
${}^{13}\textrm{C}$ nuclear spins in diamond (Fig. 1A), endowed with long
prethermal lifetimes under Floquet driving. We propose a novel use of a ”two-
frequency” Floquet drive 111The term frequency refers to the inverse periods
of the two superimposed drives, rather than their Fourier decompositions.
applied to the ${}^{13}\textrm{C}$ nuclei (see Fig. 1B): the spins are made to
prethermalize to an effective Hamiltonian $\overline{\mathcal{H}}$ (featuring
a quasi-conserved $\hat{\textbf{x}}$-magnetization [31]), under a fast drive
(with period $\tau$), while also being periodically kicked away from the
prethermal state by a slower drive with period $T$ (${=}N\tau$). The deviation
away from the $\hat{\textbf{x}}$ axis can be continuously monitored in periods
between the fast drive, allowing a means to track the full system dynamics for
long periods (¿14s), corresponding to 450 Floquet cycles or ${>}1.35\times
10^{5}$ fast pulses, without collapsing the quantum state. Compared to point-
by-point measurements in previous experiments, this allows a dimension
reduction for mapping the PDTC phase diagram. In concert with the multiple-
minute-long $T_{2}^{\prime}$ lifetimes of the ${}^{13}\textrm{C}$ nuclei, we
are able to unravel the emergent prethermal DTC order with significantly
higher clarity than previous measurements.

The two-frequency drive developed here represents a generalization of
conventional Floquet driving; at once it allows the dynamical engineering of
an emergent quasi-conserved quantity to stabilize non-equilibrium order at
arbitrary temperatures, while simultaneously exciting the DTC order. Here we
analyze the effect of driving on the long-time melting of the time-crystalline
metastable state and show that the associated heating rates match well to a
theoretical model. Overall therefore, this work enhances the Floquet
engineering toolbox for the experimental realization and observation of non-
equilibrium order beyond the paradigmatic DTC.

Fig. 2: Continuously observed prethermalization and PDTC. (A) Floquet
prethermalization of the kicked ${}^{13}\textrm{C}$ spins, shown for
$\vartheta{=}\pi/2,\tau{=}169\mu$s, $N{=}301$, and $\gamma{=}0.08\pi$. Data
(points) depict single-shot measurement of
$\left\langle\mathcal{I}_{x}\right\rangle$, for four kicks (165ms).
Oscillations depict transient approach to prethermal plateaus (flat regions).
Here, we normalize the signal upon reaching the initial prethermal state. Due
to a slight frequency offset in the pulses applied, the spin-lock axis is
slightly tilted from $\hat{\textbf{x}}$ and deviation from it results in the
signal being apparently greater than unity during the transient. Upper axis
denotes time $t$, while lower axis denotes Floquet cycle number $M$. Inset:
Zoom into region between two kicks (dashed lines) separated by 50.7ms. First
kick is preceded by a 1s-long prethermalization along $\hat{\textbf{x}}$ (not
shown). For movie of full dataset, see ref. [32] for link. (B) Single-shot
PDTC measurement for $\vartheta{=}\pi/2,\tau{=}105\mu$s, $N{=}300$, and
$\gamma{=}0.97\pi$. Panel shows the signal
$\left\langle\mathcal{I}_{x}\right\rangle$; signals from alternate kicks are
depicted by red and blue points. Total time here corresponds to $1.35\times
10^{5}$ periods of the fast drive, and ${\sim}$400 slow kicks (see Fig. 3A for
full data). PTDC decay envelope is approximately monoexponential with
$1/\mathrm{e}$ time constant ${\approx}$4.68s (dashed line). For movie of full
dataset, see ref. [33] for link. (C) PDTC phase response. For characteristic
vertical slices in Fig. 3(A), we plot the phase of the spins $\phi$ on the
Bloch sphere (see Fig. 1C) for the entire 14s period. (i) Phase response at
$\gamma{=}0.01\pi$ showing the quasi-equilibrium of the spins along
$\hat{\textbf{x}}$. (ii). At slight deviations, $\gamma{=}0.06\pi$, the spins
start leaving the prethermal axis $\hat{\textbf{x}}$, and undergo heating,
manifesting in the randomization of their phase at long times (here $t\gtrsim
10$s). (iii) PDTC response at $\gamma{=}1.01\pi$ showing stable period doubled
oscillations, represented by the $\pi$-phase shift between successive
$\gamma$-kicks. (iv) Melting of the PDTC order at $\gamma{=}1.08\pi$,
observable by the randomization of the phase at long times.

##  II System

We consider a lattice of ${}^{13}\textrm{C}$ nuclei (Fig. 1A) in diamond,
optically hyperpolarized for $t_{\textrm{pol}}{=}60$s by surrounding Nitrogen
Vacancy (NV) centers [34, 35]. Hyperpolarization yields a ${\approx}680$-fold
enhancement in ${}^{13}\textrm{C}$ magnetization over thermal equilibrium,
yielding a starting density matrix $\rho_{0}{\sim}\epsilon\mathcal{I}_{z}$,
where $\epsilon{=}0.68\%$ (see Methods). Here,
$\mathcal{I}_{\nu}{=}\sum_{j}I_{j\nu}$ with $\nu{\in}\\{x,y,z\\}$ and
$I_{j\nu}$ refer to spin-1/2 Pauli operators associated with nuclear spin $j$
[36]. The natural abundance (1%) nuclei are not spatially ordered, and are
coupled via dipolar interactions,
$\mathcal{H}_{\textrm{dd}}{=}\sum_{k<\ell}b_{k\ell}(3I_{kz}I_{\ell
z}-\vec{I_{k}}\cdot\vec{I_{\ell}})$. One can define a characteristic energy
scale via the median interspin coupling, $J{=}$0.66kHz [37]. They are also
subject to electron-mediated random on-site fields from surrounding NV and P1
paramagnetic defects, $\mathcal{H}_{z}{=}\sum_{j}c_{j}I_{zj}$ [38, 39].
Despite the concomitant position disorder and random on-site fields, the
three-dimensional long-range nature of the interactions precludes many-body
localization. The long-range interactions also make simulating the exact
Floquet dynamics for a large number of spins inaccessible using classical
computers.

Fig. 1B describes the experimental protocol — hyperpolarized
${}^{13}\textrm{C}$ nuclei are tipped along $\hat{\textbf{x}}$, and subject to
concatenated “slow” and “fast” Floquet drives, characterized by periods $T$
and $\tau$ respectively, with $T{=}N\tau$. First, the fast drive, consisting
of a train of $\vartheta({\neq}l\pi)$ pulses (Fig. 1B), engineers the
internuclear Hamiltonian so that spins, initially aligned along
$\hat{\textbf{x}}$ are rendered quasi-stationary [37]. This is accomplished by
arranging
$\mathcal{H}_{\textrm{dd}}{\rightarrow}\overline{\mathcal{H}}+d_{N}\mathcal{I}_{x}$
to leading order in the Magnus expansion, such that
[$\overline{\mathcal{H}}+d_{N}\mathcal{I}_{x},\mathcal{I}_{x}]{=}0$ [31]. This
quasi-conservation causes the spins to prethermalize along $\hat{\textbf{x}}$
due to the nonintegrable character of $\overline{\mathcal{H}}$. By contrast,
in the absence of the fast drive, evolution under $\mathcal{H}_{\textrm{dd}}$
causes system observables to rapidly (in $J^{-1}$) become indistinguishable
from a featureless infinite temperature state.

Interspersed at period $T$, the slower drive kicks the spins along the
$\hat{\textbf{y}}$ (or $\hat{\textbf{z}}$) axis with angle $\gamma$ (see Fig.
1B). The spins are allowed to prethermalize back along $\hat{\textbf{x}}$
between successive kicks (see Fig. 2A). Fig. 1C shows this visually on the
Bloch sphere (in the rotating frame) for a kick of angle
$\gamma{=}\pi{+}\varepsilon$. Therefore, in the prethermal plateau, the system
is governed by an effective Hamiltonian obtained through an inverse frequency
expansion [31]. For $\gamma{=}\\{0,\pi\\}$, although the slow
$\hat{\textbf{y}}$-kicks do not cause any extra heating, they give rise to
non-equilibrium ordered states. At $\gamma{=}\pi$ the two-cycle time-evolution
operator, $U_{F}^{2}{=}\exp[-i2T\overline{\mathcal{H}}]$, is governed by the
$\mathbb{Z}_{2}$ symmetric many-body Hamiltonian $\overline{\mathcal{H}}$,
where the $\mathbb{Z}_{2}$ symmetry is implemented by flipping the
$\hat{\textbf{x}}$-direction of all spins. This drive-induced symmetry of
$U_{F}^{2}$, together with the discrete time-translation invariance, creates a
spatio-temporal eigenstate order in $U_{F}$. Any initial state that breaks
this symmetry is forced to oscillate with period $2T$, forming a PDTC state.
In the experiment, we additionally observe that interactions stabilize a
finite region near $\gamma{=}\pi$, where a stable PDTC period doubling
response arises, with the spins flipping from $+\hat{\textbf{x}}$ to
$-\hat{\textbf{x}}$ between successive kicks (Fig. 1C).

A distinguishing feature of our experiments is the ability to quasi-
continuously track the prethermalization dynamics after each kick. The spins
are non-destructively interrogated by Nuclear Magnetic Resonance (NMR) in
$t_{\textrm{acq}}({>}0.6\tau)$ windows between the fast pulses (Fig. 1B). The
magnitude and phase of the ${}^{13}\textrm{C}$ Larmor precession is sampled
every 1ns, allowing one to reconstruct the instantaneous projections
$\left\langle\mathcal{I}_{x}(t)\right\rangle$, as well as the phase $\phi(t)$
of the spin vector in the $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane
(see Methods). In reality, the quasi-continuous time variable $t$ is
discretized in units of $\tau$. Hyperpolarization enables high signal-to-noise
(SNR${\gtrsim}10^{2}$) signal acquisition per measurement point.

It is worth emphasizing that, although rotating-frame DTCs are also observable
in our system under a suitable single-frequency drive, either no continuous
measurement can be performed or the spins rapidly decay to infinite
temperature near the DTC point $\vartheta=\pi$ [31]. Two-frequency driving
circumvents this problem since $\vartheta$ can be arbitrarily chosen (except
for $\vartheta{=}\\{0,\pi\\}$) and DTCs instead appear conditioned on
$\gamma$. Concatenated driving, therefore, engineers a separation between the
interaction driven spin-dynamics and breaking of a discrete spatio-temporal
symmetry [31]. Additionally, the ratio $N{=}T/\tau$ is tunable, affording
flexibility in exploring dynamical regimes at small and large $N$. In typical
experiments, $NM{>}1.35\times 10^{5}$ fast pulses are applied, corresponding
to $M{>}$450 Floquet cycles. Our approach portends observing the long-time and
intra-period time-crystalline dynamics continuously and without state
reinitialization. This constitutes a vast improvement in throughput with
respect to contemporary experiments, and yields access to movies that
elucidate the formation and melting of the DTC state (see SI [31]). The
resulting measurement speedup is at least $NM{>}10^{5}$; compared to
experiments probing lab-frame ($\hat{\textbf{z}}$) DTCs [28], these gains
could be as much as $NM(T_{1}/t_{\textrm{pol}}){\sim}10^{7}$.

Fig. 3: Prethermal DTC phase diagram. (A) Movie showing emergence of
prethermal DTCs. 285 traces similar to Fig. 2B are plotted stacked for
different values of kick angle $\gamma$ in $[-1.1\pi,1.1\pi]$. Colors
represent signal $\left\langle\mathcal{I}_{x}\right\rangle$ (see colorbar).
Time $N\tau$ and Floquet cycle number $M$ run vertically and are plotted on a
logarithmic scale. Data are taken to 14s ($M{=}450$ cycles). Central feature
(near $\gamma{=}0$) shows stabilization of long-time spin survival as a result
of drive-induced quasiconservation of $\hat{\textbf{x}}$-magnetization due to
Floquet prethermalization. PDTC response (striped regions) are visible near
${\gamma}={\pm}\pi$. Striped signal denotes spins inverting between
$\hat{\textbf{x}}$ and $-\hat{\textbf{x}}$ in a period-doubled fashion. Spins
flip every $N{=}300$ fast pulses, and rigid PDTC response (peaks) persist to
three decades of the fast drive. Remaining regions are characterized by rapid
spin decay due to dipolar interactions, and are indistinguishable from the
infinite-temperature state. Inset: Zoom into PDTC response in two 1s-windows
centered at $t{=}3.15$s, and $t{=}11.55$s respectively plotted on a linear
scale. Rigid DTC behavior is observable for over $14$s, even at
$\gamma{=}0.97\pi$. See ref. [33] for movie version of same dataset. (B)
Fourier transform of the data in (A), plotted with respect to the inverse
period of the slow drive $\omega=2\pi\leavevmode\nobreak\ T^{-1}$ radians.
Colors represent strength of Fourier peak intensity spanning six orders of
magnitude (colorbar). PDTC (period-doubled) order is evident by the extended
white peaks at $(\gamma,\omega)=(\pm\pi,\pi)$, while the trivial (prethermal)
phase arises for $(\gamma,\omega)=(0,0)$. Dashed lines indicate expected
Fourier peak pattern in the absence of interactions. (C-D) Numerical
simulations analogous to the experiments in A-B for $L{=}14$ interacting spins
on a pseudo random-graph with $J\tau=0.07$ and $N=300$ (see [31]). There is
excellent qualitative agreement with experimental observations.

##  III Prethermalization and discrete time crystals

Fig. 2A first clarifies the prethermalization process during the spin kicks,
essential to ultimately generate the PDTCs. Shown is a single-shot trace
plotting $\left\langle\mathcal{I}_{x}\right\rangle$ after every fast pulse. We
display the Floquet cycle number $M$ on the lower $\hat{\textbf{x}}$-axis and
absolute time on the upper $\hat{\textbf{x}}$-axis for an exemplary 165ms
window. Here $\vartheta{=}\pi/2$, and the $\gamma$-kicks are applied every
50.7ms (denoted by the dashed lines). Each $\gamma$-kick is associated with
transient dynamics of the coupled ${}^{13}\textrm{C}$ nuclei as they
prethermalize along $\hat{\textbf{x}}$, producing the flat plateau-like
regions shown. The inset Fig. 2A(i) zooms into one representative transient;
the oscillation period here is set approximately by the number of fast pulses
required to complete a $2\pi$ rotation. The high temporal resolution revealed
by the high SNR, and the ability of the ${}^{13}\textrm{C}$ nuclear system to
sustain a large number of $\gamma$-kicks (see also Fig. 3A) allows us to track
the kicked prethermalization dynamics for long periods.

Employing an exemplary choice of flip-angle $\gamma{=}0.97\pi$ slightly away
from the perfect DTC point, Fig. 2B demonstrates generation of stable DTC
order, exhibiting period doubling in
$\left\langle\mathcal{I}_{x}(t)\right\rangle$ during the application of
$M{>}$450 $\gamma$-kicks (see Fig. 3A for full data). The data were collected
in a single run of the experiment. Red (blue) colors here represent odd (even)
$\gamma$-kicks respectively, and prethermal plateaus separate successive
kicks. These data comprise ${>}1.35\times 10^{5}$ fast pulses, and rigid DTC
behavior is observable for ${>}$14s (see also Fig. 3A). The decay of the
signal is approximately mono-exponential with a $1/\mathrm{e}$ lifetime
$t{\approx}4.68$s (corresponding to $M_{e}{\approx}149$ Floquet cycles at
$\gamma{=}0.97\pi$), making it amongst the longest DTCs observed in the
literature. Moreover, the $Jt{\approx}10^{4}$ value here is considerably
beyond state-of-art for systems exhibiting DTC order [12, 11, 16, 17, 29],
demonstrating an ability to probe long-time dynamics in our system. This long-
time stability can be attributed to the emergent quasi-conservation of
$\hat{\textbf{x}}$-magnetization under the evolution engineered by the two-
frequency drive. In particular, our driving protocol allows the formation and
observation of DTC order, even if temperature in the prethermal plateau is
(close to) infinite [40], but also at room temperature, which may be well
above the critical temperature associated with symmetry breaking [31]. The
long-range nature of the spin-spin interaction suggests critical DTC order
[9]. On the other hand, the lifetime of the DTC order is parametrically
controlled by the frequency (of switching) of the employed drives (see sec. IV
and SI [31]). Thus, the observed DTC order corresponds to a critical
prethermal DTC.

Continuous observation yields an insightful view into the thermalization
dynamics away from the stable points, a challenging task in other experimental
systems. This is demonstrated in data focusing on the instantaneous phase
$\phi(t){=}\tan^{-1}(\left\langle\mathcal{I}_{y}\right\rangle/\left\langle\mathcal{I}_{x}\right\rangle)$
of the spins in the $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane after
every fast pulse. This is captured by the points in Fig. 2C for the full 14s
experiment, and considering four representative constant-$\gamma$ values (see
Fig. 3A for full data). When $\gamma{=}0.01\pi$, the spins are locked at
$\phi{\approx}0$, reflecting prethermalization along $\hat{\textbf{x}}$. A
slight deviation $\gamma{=}0.06\pi$, reveals transient oscillations in $\phi$
with every $\gamma$-kick but no sign inversion. The transients result in the
observed data spreading around $\phi{=}0$. Fig. 2C(ii) therefore permits
visualization of the ”melting” of the prethermal
$\hat{\textbf{x}}$-magnetization order to infinite temperature, where $\phi$
becomes random; this is observable at $t{\gtrsim}10$s. Analogously, the PDTC
order (Fig. 2C(iii)), here at $\gamma{=}1.01\pi$, is visible as characteristic
$\pi$-phase switching between successive kicks. The last panel (Fig. 2C(iv))
denotes $\gamma{=}1.08\pi$, when again PDTC melting can be observed via the
phase randomization at long times.

Collating 285 such data traces while varying angle
$\gamma{\in}([-1.1\pi,1.1\pi])$, it is possible to construct a movie of the
kicked prethermalizing spins (see ancillary .gif files in SI [31]). The result
is plotted in Fig. 3A, where we display
$\left\langle\mathcal{I}_{x}\right\rangle$ (cf. colorbar), with each vertical
slice corresponding to measurement data as in Fig. 2B. Left and right vertical
axes here refer to the Floquet cycle number $M$ and absolute time,
respectively, on a logarithmic scale. The data highlight the almost three
decades in the slow drive kicks, for which we observe the PDTC order. The
transition into and out of the finite PDTC regions is clearly evident, and not
easily accessible in other experiments, allowing the ability to precisely
characterize the heating dynamics (see Fig. 4), and map the entire non-
equilibrium phase diagram. Near $\gamma{=}0$, we observe a suppression of
heating due to the emergent drive-induced quasiconservation of the
$\hat{\textbf{x}}$-magnetization; the low-frequency kicks approximately wrap
to the identity, and the effective Hamiltonian
$\mathcal{H}_{\textrm{eff}}{\approx}\overline{\mathcal{H}}$ is built solely
from the high-frequency $\hat{\textbf{x}}$-drive. Instead, near
$\gamma{=}\pm\pi$, we observe a region corresponding to PDTC order, where
there is a regular switching of the magnetization from $\hat{\textbf{x}}$ to
-$\hat{\textbf{x}}$ with every slow kick (here 300 fast pulses), a signature
of period doubling response. Interspin interactions play a crucial role to
preserve a uniform switching frequency away from $\gamma\\!=\\!\pi$. The
insets (Fig. 3A(i) and (ii)) show zooms into PDTC regions in a 1.05s window at
$t{=}3.15$s and $t{=}11.55$s respectively. We plot them here on a linear scale
with time for clarity. While faint, the stable periodic DTC response is
markedly clear even in Fig. 3A(ii), highlighting the high SNR in the
experiment.

A complementary view of Fig. 3A is presented in Fig. 3B, where we consider the
Fourier transform of $\left\langle\mathcal{I}_{x}(M)\right\rangle$, the mean
signal value between successive $\gamma$-kicks for each value of $\gamma$
(vertical slices in Fig. 3A). Plotted is the corresponding Fourier intensity
on logarithmic scale, where a span to six orders of magnitude is visible. The
PDTC response appears as a sharp period-doubling peak in frequency at
$(\gamma,\omega)=(\pm\pi,\pi)$. Using a 20% magnitude threshold, we estimate
the finite $\gamma$-extent of the rigid PDTC regions to
$\Delta\varepsilon{=}\pm 0.2\pi$ about $\gamma{=}\pm\pi$. Similarly, the long-
lived prethermal phase at $\gamma{=}0$ appears as a sharp peak at
$(\gamma,\omega)=(0,0)$. For comparison, the dark regions correspond to rapid
state decay near $\gamma{=}\pm\pi/2$, where the effective Hamiltonian no
longer features a quasi-conserved $\hat{\textbf{x}}$-magnetization [31].
Finally, we note in the dashed lines in Fig. 3B the expected position of the
Fourier peaks in the absence of interactions, $J{=}0$. In this diagram, the
peak positions trace a rhombic pattern (dashed lines). The difference in the
experimental data is evident, which indicates the role played by interactions
in our system.

Comparing the experimental measurements to corresponding numerical simulations
(see [31] for details), we find an excellent qualitative agreement: In Fig.
3C-D we display the exact simulation results matching experimental conditions
in the data in Fig. 3A-B. Interestingly, we observe that the long-range
interactions, together with the random spin positions, induce a self-averaging
effect in the simulations so that reliable theoretical results can already be
obtained for moderately small system sizes (here $L{=}14$ spins). By contrast,
the experimental platform comprises a cluster of about $L{\sim}10^{4}$
interacting spins [39] which outcompetes the numerically reachable system
sizes by three orders of magnitude. The presence of a large number of
interacting degrees of freedom is crucial for experimentally observing
collective statistical mechanics phenomena such as symmetry breaking and
thermalization. Thus, our results indicate that Floquet-engineered
${}^{13}\textrm{C}$ nuclei can serve as a competitive quantum simulator of
thermalizing spin dynamics.

The observed PDTC behavior is also insensitive to the initial state. While we
lack microscopic control over the initial state, different and highly non-
trivial initial states can be obtained from letting the initial hyperpolarized
density matrix $\rho_{0}$ evolve under the system Hamiltonian
$\mathcal{H}_{\mathrm{dd}}+\mathcal{H}_{z}$ up to times $\sim\\!J^{-1}$. For a
set of such states we find comparable results between experiment and
simulation [31]. Our system thus satisfies all required landmarks of PDTCs: a
parametrically long lived prethermal window featuring spatio-temporal symmetry
breaking, which is rigid over a finite $\gamma$-region, and insensitive to
fine-tuned initial states.

Fig. 4: Experimental characterization of PDTC rigidity and prethermal
lifetimes. (A-B) System heating rates estimated by the inverse of a
$1/\mathrm{e}$-threshold lifetime. We observe the heating rates follow a power
law $\Gamma{\propto}g\epsilon^{\lambda}+\Gamma_{\mathrm{min}}$ [see text],
with an exponent $\lambda{=}2.06$ for prethermal phase (A) and
$\lambda{=}2.21$ for PDTC phase (B), agreeing well with numerical simulations
[31]. Solid line is a fit and dashed lines are error bars at two standard
deviations. (C-D) Prethermal lifetimes for different values of $N$ and $\tau$.
Dashed lines correspond to power-law fits with exponents $-0.44$ (C) and
$-1.90$ (D). Increasing the drive frequency leads to an increase in lifetime
of the prethermal DTC state.

##  IV  Melting of prethermal order

To quantify the stability of the PDTC order away from the stable point at
$\gamma{=}\pi$, we investigate the influence of finite $\varepsilon$ on its
lifetime, which we define as the $1/\mathrm{e}$ decay time of the signal. In
turn, the inverse lifetime defines the associated heating rate. The long-range
character of the dipolar interactions leads to a logarithmically divergent
total energy of the system in three dimensions that, when combined with the
lack of Lieb-Robinson bounds [23], makes the theoretical analysis of the
heating rates for large systems difficult. For these reasons, to quantify the
heating rates of the two-frequency drive, we perform a series of numerical
simulations, and compare the results against experimental observations [31].

Note that, even at $\gamma{=}\pi$ (where the slow drive does not contribute to
heating) the system slowly heats up due to energy absorption resulting from
the fast (spin-locking) $\hat{\textbf{x}}$-drive. This provides a lower bound
for relevant heating rates, $\Gamma_{\mathrm{min}}(\tau)$, which we find
scales as a power-law of the $\hat{\textbf{x}}$-drive period,
$\Gamma_{\mathrm{min}}{\propto}(J\tau)^{2}$ [37]. Finite $\varepsilon$ opens
an additional channel for the system to absorb energy and, eventually, the
heating rates of both drives conspire to yield a combined overall heating rate
$\Gamma(\tau,\varepsilon)$. In Fig. 4 we extract the heating rates associated
to finite $\varepsilon$, where we display $\Gamma-\Gamma_{\mathrm{min}}$ for
$|\left\langle\mathcal{I}_{x}(t)\right\rangle|$ away from the stable points
($\gamma{=}\\{0,\pi\\})$ obtained from Fig. 3A-B, and plotted on a logarithmic
scale with respect to both heating rate and deviation $\varepsilon$. The
curves indicate a parametrically controllable power-law heating with an
exponent ${\approx}2.21$, consistent with Fermi’s Golden rule. Intriguingly,
we observe that the power law heating rates are close to identical for both
stable points $\gamma{=}\\{0,\pi\\}$, where the relevant timescales are
governed by a Lorentzian
$\Gamma^{-1}{\propto}[g\varepsilon^{2}+\Gamma_{\mathrm{min}}]^{-1}$ for some
system-dependent constant $g$; this behavior is also borne out in numerical
simulations [31].

For the observed prethermal DTC lifetime to be controllable, relevant life-
times are expected to increase with increasing drive frequencies. Indeed, the
heating time in units of Floquet cycles $\Gamma^{-1}/N$ depends sensitively on
the frequency of the Floquet drive $\omega_{F}=2\pi/(N\tau)$. In particular,
at $\varepsilon=0$ our model predicts $\Gamma^{-1}/N\sim 1/(NJ\tau^{2})$.
Increasing the frequency of switching by tuning $N$ ($\tau$), is expected to
lead to a linear (quadratic) increase of $\Gamma^{-1}/N$, which we have
confirmed numerically for our system [31]. In Fig. 4C-D we display
experimental measurements of the dependence of $\Gamma^{-1}/N$ as a function
of $N$ and $\tau$, for finite $\varepsilon$. Although a crisp $1/N$ dependence
is washed out in the presence of finite $\varepsilon$ (see also [31]), we
observe a clear increase in the lifetime. For sufficiently large $\tau$ in the
timescale-separated regime $T/\tau\gg 1$, $\Gamma^{-1}(\tau)/N$ as a function
of $\tau$ agrees well with a power law with an exponent close to $-2$, while
at very small $\tau$ the DTC lifetime increase comes to a halt. This is to be
expected for fixed $\varepsilon$, since decreasing the value of $\tau$ reduces
the many-body nature of the effective Hamiltonian (see [31]). Note that DTC
order is a many-body effect that relies on spontaneous spatio-temporal
symmetry breaking in interacting systems. Thus, when the many-body nature gets
gradually reduced, the stability and rigidity of the DTC order are expected to
decline accordingly.

##  V  Two-frequency Floquet engineering

When the condition $N{=}T/\tau{\gg}1$ is not met, time-scale separation
between the slow and the fast drive is violated, and a complex interplay
between the two drives emerges. In this regime, the two drives cannot be
treated independently of one another; instead, they mutually influence each
other leading to novel effective Hamiltonians that sensitively depend on $N$.
Consequently, the corresponding heating diagrams exhibit significant
differences, as compared to the time-scale separated case, but also for
different values of $N$. In Fig. 5 we show the time-evolution of the
$\hat{\textbf{x}}$-magnetization as a function of $\gamma$ for $N=8$ and
$N=9$. Three features are particularly noteworthy: (i) in contrast to the
timescale separated case, the heating dynamics around $\gamma=\pi$ differs
significantly from that around $\gamma=0$. In particular, (ii) far away from
$\gamma=\pi$, heating depends sensitively on the specific value of $N$, where
even a minimal change (in $N$) can induce completely different heating
behaviours (compare Fig. 5A-B). These features can be explained by the
different effective Hamiltonians (see [31]) forming at different $N$,
$\gamma$, respectively. (iii) Even though interference effects induce new
effective Hamiltonians with case-specific heating properties, the formation of
time-crystalline order is stable against these deformations (Fig. 5C-D).
However, the relevant lifetimes and regimes of rigidity depend sensitively on
$N$ [31]: changing $N=8$ into $N=9$ amounts in a (late) lifetime increase from
$270$ to $396$ Floquet cycles (cf. Fig. 5C-D)). These results provide a proof-
of-principle example for the fine interplay between the two drives which
offers a versatile tool to engineer new kinds of effective Hamiltonians with
orchestrated physical properties, such as DTC order and beyond.

Fig. 5: Experimental characterisation of the heating dynamics at small N.
Complex interplay between the fast and slow periodic drives is seen when the
condition $N=T/\tau\gg 1$ is not met. (A) Movie showing time-series data from
103 experiments with different values of flip angle $\gamma$ in $[0,1.1\pi]$
and $N=8$ $\hat{\textbf{x}}$-pulses between $\hat{\textbf{y}}$-pulses. Colors
represent absolute value of signal $|\langle\mathcal{I}_{x}\rangle|$ (see
colorbar). Floquet cycle number $M$ runs vertically on a logarithmic scale.
Data are taken past $10^{4}$ Floquet cycles. (B) Movie similar to (A) but with
$N=9$ $\hat{\textbf{x}}$-pulses between the $\hat{\textbf{y}}$-pulses. In both
movies (A) and (B), $\hat{\textbf{x}}$-pulse flip angle was calibrated to
$\vartheta=\pi/2$. (C-D) Line cut of (A-B), respectively, at
$\gamma=1.011\pi$, shows a clear DTC signal with even (odd) cycles shown in
red (blue). Late lifetimes (measured starting from Floquet cycle 100, after
the initial transient regime) are 270 Floquet cycles for $N=8$ (C) and 396
Floquet cycles for $N=9$ (D). Due to the complexity of signal processing we
plot absolute values in (A-B) and unwrap the phase information only for the
line cuts shown in (C-D).

##  VI Conclusions and Outlook

Summarizing, we have observed critical prethermal discrete time crystals in
dipolar-coupled nuclear spins in a bulk three-dimensional solid. We developed
a novel protocol to excite and observe the PDTC formation and melting using a
concatenated two-frequency Floquet drive. Parametric control over both drive
frequencies allows us to reach PDTC lifetimes up to 396 Floquet cycles,
observable in a single run of the experiment. This experimental advance
unveils properties of the PDTC with a high degree of clarity, including its
rigidity and melting characteristics, and to map the entire prethermal phase
diagram of the system. Our measurements are in excellent agreement with
numerical simulations, and approximate theoretical predictions [31].

Our study greatly expands the Floquet engineering toolkit: it portends multi-
frequency concatenated drives to excite and stabilize quantum matter far from
equilibrium by engineering emergent quasi-conservation laws that offer
protection against immediate high-temperature melting. The novel Floquet
driving protocol we introduced paves the way to experimentally realizing
intriguing non-equilibrium states, potentially including PDTCs with a stable
fractional subharmonic response [41], or PDTCs and topological prethermal
states stabilized by quasiperiodic driving [42]. An extension of this work to
multi-frequency driving is well within the scope of present-day experimental
capability, and can be used to introduce more degrees of freedom to Floquet
engineering. We envision application of these ideas in quantum simulation and
sensing in AMO platforms as well as with hyperpolarized prethermal spins in
solid-state systems.

##  VII Acknowledgments

We thank J. Bardarson, M. Heyl, C. von Keyserlingk, D. Luitz, R. Moessner, J.
Reimer, and D. Suter for valuable discussions. A.A. acknowledges funding from
ONR under contract N00014-20-1-2806. C.F. acknowledges support from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreements No. 679722 and No.
101001902). M.B. was supported by the Marie Sklodowska-Curie grant agreement
No 890711, and the Bulgarian National Science Fund within National Science
Program VIHREN, contract number KP-06-DV-5 (until 25.06.2021). Computational
work reported on in this paper was performed on the Würzburg HPC cluster.

Data availability – Raw data and scripts to process them is available upon
request.

Author contributions – W.B., C.F., M.B., and A.A. conceived the research.
W.B., A.P., E.S., A.Ak., J.A., S.C., P.R, E.D., and A.A. set up the
experimental apparatus, performed measurements and analyzed the data. C.F.
performed the numerical simulations and the perturbative analysis. A.A. and
M.B. supervised the experiment and the theory work.

Competing interests – Authors declare no competing interests.

##  VIII Methods

1\. Sample and hyperpolarization methodology: Here we employ a CVD grown
single crystal of diamond with a defect density of 1ppm NV centers and natural
abundance ${}^{13}\textrm{C}$. The sample is placed with its [100] face
parallel to the hyperpolarization and readout magnetic fields ( 38mT and 7T
respectively). ${}^{13}\textrm{C}$ nuclei are hyperpolarized via the NV
centers via a protocol described in Ref. [37] involving continuous optical
pumping and swept microwave (MW) irradiation. Polarization is transferred from
the NV electrons to proximal ${}^{13}\textrm{C}$ nuclei via a spin-ratchet
mechanism that involves traversals of Landau-Zener anticrossings in the
rotating frame of the swept MWs. The traversal probabilities are nuclear spin-
state selective and are biased towards one nuclear spin orientation; this
results in hyperpolarization. Spin diffusion serves to transfer this
polarization to nuclear spins in the bulk lattice. Typical bulk polarization
levels reached are about 0.6%.

2\. Data collection and processing: The data processing pipeline follows a
similar approach as described in Ref. [37]. The NMR signal is sampled
continuously in $t_{\textrm{acq}}$ windows between the fast spin-lock pulses,
at a sampling rate of 1GS/s via a Tabor Proteus arbitrary waveform
transceiver. Continuous observation exploits the fact that the NMR coil
produces no backaction on the spins. In typical experiments, the pulses are
spaced apart by $\tau{=}105\;\mu$s and the acquisition windows are
$t_{\textrm{acq}}{=}64\;\mu$s. The ${}^{13}\textrm{C}$ Larmor precession (at
${\approx}$75MHz) is heterodyned to 20MHz prior to digitization. For each
$t_{\textrm{acq}}$ acquisition window, we take a Fourier transform and extract
the magnitude and phase of the 20MHz peak. This corresponds to the application
of a digital bandpass filter with a linewidth of
$t_{\textrm{acq}}^{-1}{\approx}31.2$kHz. Fast digitization hence yields
signal-to-noise (SNR) gains, and typical SNR per point (detection window) is
${\gtrsim}10^{2}$. For the 14s long acquisition in the paper, there are
$\sim$135k such data collection windows.

3\. Extraction of amplitude and phase in Fig. 2 and Fig. 3: The magnitude and
phase of the Fourier transform of the heterodyned precession in each
$t_{\textrm{acq}}$ readout window report respectively on the magnitude and
phase of the spin vector in the $\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$
plane in the lab frame. This corresponds to magnitude
$S_{L}{=}[\left\langle\mathcal{I}_{x}^{2}\right\rangle+\left\langle\mathcal{I}_{y}^{2}\right\rangle]^{1/2}$,
and phase
$\phi_{L}{=}\left\langle\mathcal{I}_{y}\right\rangle/\left\langle\mathcal{I}_{x}\right\rangle$,
where subscript $L$ here refers to the lab frame. It is more convenient to to
instead obtain the phase of the spins $\phi_{R}{=}\phi$ in the rotating frame.
To do this, we note that the phase values $\phi_{L}$ obtained in successive
$t_{\textrm{acq}}$ windows just differ by the (trivial) phase accrued during
the $t_{p}$ spin-locking pulse. Subtracting this global phase allows us to
extract $\phi$, which in combination with the magnitude signal then allows us
to extract the survival probability along the $\hat{\textbf{x}}$ direction in
the rotating frame, $\left\langle\mathcal{I}_{x}(t)\right\rangle$ that we
display in Fig. 2 and Fig. 3 of the main text.

4\. Numerical simulations: To numerically simulate the many-body dynamics of
dipolar interacting ${}^{13}\textrm{C}$ nuclei, we design random graphs of
$L{=}14$ and $L{=}16$ interacting spins-$1/2$ and perform exact time evolution
with up to $10^{6}$ $\hat{\textbf{x}}$-kicks (corresponding to $\approx
10^{3}$ $\hat{\textbf{y}}$-kicks) based on an OMP-parallelized Krylov method
using the open-source python package QuSpin [43]. For further details we refer
the reader to the supplementary material [31].

5\. Comparison with previous work: To date DTC order has been studied in
various systems ranging from cold atoms [12, 29] over superconducting qubits
[17] to systems based on NV centers [11, 16]. The DTC order observed in these
works can be separated into two groups: many-body localized DTCs and
prethermal DTCs. While many-body localized DTCs are assumed to be infinitely
long-lived in the absence of decoherence, prethermal DTCs, as in our
experiments, are ultimately limited by the lifetime of the prethermal plateau.
However, in reality also many-body localized DTCs are subject to decoherence
due to technical limitations. Remarkably, even though our DTC order is of
different physical origin, we find comparable lifetimes to state-of-the-art
many-body-localized DTCs. Moreover, our lifetimes exceed those reported for
prethermal DTCs, both, in units of Floquet cycles as well as absolute time.

Apart from these benchmark parameters, our system shares ingredients with the
systems investigated in Refs. [11] and [16]. Like in our work, these works
investigate systems based on NV centers: Reference [16] examines a quasi one-
dimensional system of 9 dipolar coupled 13C nuclear spins in the many-body
localized regime, while in Ref. [11] effective two-level systems of electronic
states in NV centers are used to implement DTC order. In contrast, our (three-
dimensional) system consists of $10^{3}-10^{4}$ dipolar coupled 13C nuclear
spins. In particular, in comparison to Ref. [11] the normalized interaction
strength $\langle J\rangle/\gamma_{n}^{2}$ in our system is increased by a
factor $4.5\times 10^{4}$ [here $\gamma_{e/n}$ denotes the gyromagnetic ratio
of electronic/nuclear spins]; moreover, at comparable normalized driving
strength, we obtain an improvement of normalized spin lifetimes by a factor of
$5.7\times 10^{2}$ (see also [31] for a detailed table of relevant system
parameters). Induced by the two-frequency drive, the DTC order we observe is
prethermal, i.e., our DTC lifetime is parametrically controlled by the drive
frequency, while no such feature is present in Ref. [11] where a single-
frequency drive is used.

## References

  * Sacha and Zakrzewski [2017] K. Sacha and J. Zakrzewski, Time crystals: a review, Reports on Progress in Physics 81, 016401 (2017).
  * Khemani  _et al._ [2019] V. Khemani, R. Moessner, and S. Sondhi, A brief history of time crystals, arXiv preprint arXiv:1910.10745 (2019).
  * Else  _et al._ [2020a] D. V. Else, C. Monroe, C. Nayak, and N. Y. Yao, Discrete time crystals, Annual Review of Condensed Matter Physics 11, 467 (2020a).
  * Sacha [2015] K. Sacha, Modeling spontaneous breaking of time-translation symmetry, Phys. Rev. A 91, 033617 (2015).
  * Khemani  _et al._ [2016] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems, Phys. Rev. Lett. 116, 250401 (2016).
  * Else  _et al._ [2016] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Phys. Rev. Lett. 117, 090402 (2016).
  * von Keyserlingk  _et al._ [2016] C. W. von Keyserlingk, V. Khemani, and S. L. Sondhi, Absolute stability and spatiotemporal long-range order in floquet systems, Phys. Rev. B 94, 085112 (2016).
  * Yao  _et al._ [2017] N. Y. Yao, A. C. Potter, I.-D. Potirniche, and A. Vishwanath, Discrete time crystals: Rigidity, criticality, and realizations, Phys. Rev. Lett. 118, 030401 (2017).
  * Ho  _et al._ [2017] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin, Critical time crystals in dipolar systems, Phys. Rev. Lett. 119, 010602 (2017).
  * Liao  _et al._ [2019] L. Liao, J. Smits, P. van der Straten, and H. T. C. Stoof, Dynamics of a space-time crystal in an atomic bose-einstein condensate, Phys. Rev. A 99, 013625 (2019).
  * Choi  _et al._ [2017] S. Choi, J. Choi, R. Landig, G. Kucsko, H. Zhou, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, V. Khemani, _et al._ , Observation of discrete time-crystalline order in a disordered dipolar many-body system, Nature 543, 221 (2017).
  * Zhang  _et al._ [2017] J. Zhang, P. Hess, A. Kyprianidis, P. Becker, A. Lee, J. Smith, G. Pagano, I.-D. Potirniche, A. C. Potter, A. Vishwanath, _et al._ , Observation of a discrete time crystal, Nature 543, 217 (2017).
  * Pal  _et al._ [2018] S. Pal, N. Nishad, T. S. Mahesh, and G. J. Sreejith, Temporal order in periodically driven spins in star-shaped clusters, Phys. Rev. Lett. 120, 180602 (2018).
  * Smits  _et al._ [2018] J. Smits, L. Liao, H. T. C. Stoof, and P. van der Straten, Observation of a space-time crystal in a superfluid quantum gas, Phys. Rev. Lett. 121, 185301 (2018).
  * Rovny  _et al._ [2018a] J. Rovny, R. L. Blum, and S. E. Barrett, ${}^{31}\mathrm{P}$ nmr study of discrete time-crystalline signatures in an ordered crystal of ammonium dihydrogen phosphate, Phys. Rev. B 97, 184301 (2018a).
  * [16] J. Randall, C. Bradley, F. van der Gronden, A. Galicia, M. Abobeih, M. Markham, D. Twitchen, F. Machado, N. Yao, and T. Taminiau, Many-body-localized discrete time crystal with a programmable spin-based quantum simulator, Science 374, 1474.
  * Mi  _et al._ [2021] X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, R. Babbush, _et al._ , Time-crystalline eigenstate order on a quantum processor, Nature , 1 (2021).
  * Singh  _et al._ [2019] K. Singh, C. J. Fujiwara, Z. A. Geiger, E. Q. Simmons, M. Lipatov, A. Cao, P. Dotti, S. V. Rajagopal, R. Senaratne, T. Shimasaki, M. Heyl, A. Eckardt, and D. M. Weld, Quantifying and controlling prethermal nonergodicity in interacting floquet matter, Phys. Rev. X 9, 041021 (2019).
  * Rubio-Abadal  _et al._ [2020] A. Rubio-Abadal, M. Ippoliti, S. Hollerith, D. Wei, J. Rui, S. L. Sondhi, V. Khemani, C. Gross, and I. Bloch, Floquet prethermalization in a bose-hubbard system, Phys. Rev. X 10, 021044 (2020).
  * Peng  _et al._ [2021] P. Peng, C. Yin, X. Huang, C. Ramanathan, and P. Cappellaro, Floquet prethermalization in dipolar spin chains, Nature Physics 17, 444 (2021).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. Huveneers, Phys. Rev. Lett. 115, 256803 (2015).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Phys. Rev. Lett. 116, 120401 (2016).
  * Machado  _et al._ [2020] F. Machado, D. V. Else, G. D. Kahanamoku-Meyer, C. Nayak, and N. Y. Yao, Long-range prethermal phases of nonequilibrium matter, Phys. Rev. X 10, 011043 (2020).
  * Else  _et al._ [2017] D. V. Else, B. Bauer, and C. Nayak, Prethermal phases of matter protected by time-translation symmetry, Phys. Rev. X 7, 011026 (2017).
  * Pizzi  _et al._ [2021a] A. Pizzi, A. Nunnenkamp, and J. Knolle, Classical prethermal phases of matter, Phys. Rev. Lett. 127, 140602 (2021a).
  * Ye  _et al._ [2021] B. Ye, F. Machado, and N. Y. Yao, Floquet phases of matter via classical prethermalization, Phys. Rev. Lett. 127, 140603 (2021).
  * Yao  _et al._ [2020] N. Y. Yao, C. Nayak, L. Balents, and M. P. Zaletel, Classical discrete time crystals, Nature Physics 16, 438 (2020).
  * Rovny  _et al._ [2018b] J. Rovny, R. L. Blum, and S. E. Barrett, Observation of discrete-time-crystal signatures in an ordered dipolar many-body system, Phys. Rev. Lett. 120, 180603 (2018b).
  * Kyprianidis  _et al._ [2021] A. Kyprianidis, F. Machado, W. Morong, P. Becker, K. S. Collins, D. V. Else, L. Feng, P. W. Hess, C. Nayak, G. Pagano, N. Y. Yao, and C. Monroe, Observation of a prethermal discrete time crystal, Science 372, 1192 (2021).
  * Note [1] The term frequency refers to the inverse periods of the two superimposed drives, rather than their Fourier decompositions.
  * [31] See Supplemental Material.
  * PDT [2021a] Movie of full dataset from Fig. 2A of main text: https://youtu.be/61ZqLgbCuyo (2021a).
  * PDT [2021b] Movie of full dataset from Fig. 2B and Fig. 3A of main text (first 55 Floquet cycles): https://youtu.be/m5iASnBZ9oo (2021b).
  * Ajoy  _et al._ [2018a] A. Ajoy, K. Liu, R. Nazaryan, X. Lv, P. R. Zangara, B. Safvati, G. Wang, D. Arnold, G. Li, A. Lin, _et al._ , Orientation-independent room temperature optical 13c hyperpolarization in powdered diamond, Sci. Adv. 4, eaar5492 (2018a).
  * Ajoy  _et al._ [2018b] A. Ajoy, R. Nazaryan, K. Liu, X. Lv, B. Safvati, G. Wang, E. Druga, J. Reimer, D. Suter, C. Ramanathan, _et al._ , Enhanced dynamic nuclear polarization via swept microwave frequency combs, Proceedings of the National Academy of Sciences 115, 10576 (2018b).
  * Duer [2004] M. Duer,  _Introduction to Solid-State NMR Spectroscopy_ (John Wiley $\&$ Sons, 2004).
  * Beatrez  _et al._ [2021] W. Beatrez, O. Janes, A. Akkiraju, A. Pillai, A. Oddo, P. Reshetikhin, E. Druga, M. McAllister, M. Elo, B. Gilbert, D. Suter, and A. Ajoy, Floquet prethermalization with lifetime exceeding 90 s in a bulk hyperpolarized solid, Phys. Rev. Lett. 127, 170603 (2021).
  * Reynhardt [2003] E. Reynhardt, Spin lattice relaxation of spin-ï¿½ nuclei in solids containing diluted paramagnetic impurity centers. i. zeeman polarization of nuclear spin system, Concepts in Magnetic Resonance Part A 19A, 20 (2003).
  * Ajoy  _et al._ [2019] A. Ajoy, B. Safvati, R. Nazaryan, J. T. Oon, B. Han, P. Raghavan, R. Nirodi, A. Aguilar, K. Liu, X. Cai, X. Lv, E. Druga, C. Ramanathan, J. A. Reimer, C. A. Meriles, D. Suter, and A. Pines, Hyperpolarized relaxometry based nuclear t1 noise spectroscopy in diamond, Nature Communications 10, 5160 (2019).
  * Luitz  _et al._ [2020] D. J. Luitz, R. Moessner, S. L. Sondhi, and V. Khemani, Prethermalization without temperature, Phys. Rev. X 10, 021046 (2020).
  * Pizzi  _et al._ [2021b] A. Pizzi, J. Knolle, and A. Nunnenkamp, Higher-order and fractional discrete time crystals in clean long-range interacting systems, Nature communications 12, 1 (2021b).
  * Else  _et al._ [2020b] D. V. Else, W. W. Ho, and P. T. Dumitrescu, Long-lived interacting phases of matter protected by multiple time-translation symmetries in quasiperiodically driven systems, Phys. Rev. X 10, 021032 (2020b).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: a Python Package for Dynamics and Exact Diagonalisation of Quantum Many Body Systems part I: spin chains, SciPost Phys. 2, 003 (2017).
  * Khodjasteh and Lidar [2005] K. Khodjasteh and D. A. Lidar, Fault-tolerant quantum dynamical decoupling, Phys. Rev. Lett. 95, 180501 (2005).
  * Witzel and Das Sarma [2007] W. M. Witzel and S. Das Sarma, Concatenated dynamical decoupling in a solid-state spin bath, Phys. Rev. B 76, 241303 (2007).
  * Cai  _et al._ [2012] J. Cai, B. Naydenov, R. Pfeiffer, L. P. McGuinness, K. D. Jahnke, F. Jelezko, M. B. Plenio, and A. Retzker, Robust dynamical decoupling with concatenated continuous driving, New Journal of Physics 14, 113023 (2012).
  * Hayes  _et al._ [2014] D. Hayes, S. T. Flammia, and M. J. Biercuk, Programmable quantum simulation by dynamic hamiltonian engineering, New Journal of Physics 16, 083027 (2014).
  * Vajna  _et al._ [2018] S. Vajna, K. Klobas, T. Prosen, and A. Polkovnikov, Replica resummation of the baker-campbell-hausdorff series, Phys. Rev. Lett. 120, 200607 (2018).
  * Fleckenstein and Bukov [2021a] C. Fleckenstein and M. Bukov, Prethermalization and thermalization in periodically driven many-body systems away from the high-frequency limit, Phys. Rev. B 103, L140302 (2021a).
  * Fleckenstein and Bukov [2021b] C. Fleckenstein and M. Bukov, Thermalization and prethermalization in periodically kicked quantum spin chains, Phys. Rev. B 103, 144307 (2021b).
  * PDT [2021c] Movie of full dataset from Fig. 2B and Fig. 3A of main text (first 25 Floquet cycles): https://youtu.be/_PRlBEPO_54 (2021c).
  * PDT [2021d] Movie of full dataset from Fig. 2B and Fig. 3A of main text (first 155 Floquet cycles): https://youtu.be/YHm4wRz0sfQ (2021d).
  * Magnus [1954] W. Magnus, On the exponential solution of differential equations for a linear operator, Communications on Pure and Applied Mathematics 7, 649 (1954).
  * Haeberlen [1976] U. Haeberlen,  _High Resolution NMR in Solids: Selective Averaging_ (Academic Press Inc., New York, 1976).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J. Oteo, and J. Ros, The magnus expansion and some of its applications, Physics Reports 470, 151 (2009).
  * Ajoy  _et al._ [2020] A. Ajoy, R. Nirodi, A. Sarkar, P. Reshetikhin, E. Druga, A. Akkiraju, M. McAllister, G. Maineri, S. Le, A. Lin, A. M. Souza, C. A. Meriles, B. Gilbert, D. Suter, J. A. Reimer, and A. Pines, Dynamical decoupling in interacting systems: applications to signal-enhanced hyperpolarized readout, arXiv:2008.08323 (2020).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Ikeda and Polkovnikov [2021] T. N. Ikeda and A. Polkovnikov, Fermi’s golden rule for heating in strongly driven floquet systems, Phys. Rev. B 104, 134308 (2021).

Supplementary Information
Observation of a critical prethermal discrete time crystal created by two-
frequency driving
William Beatrez1,∗, Christoph Fleckenstein2,∗, Arjun Pillai1, Erica Sanchez1,
Amala Akkiraju1,

Jesus Alcala,1 Sophie Conti,1 Paul Reshetikhin,1 Emanuel Druga,1, Marin
Bukov3,4, and Ashok Ajoy1,5

1Department of Chemistry, University of California, Berkeley, Berkeley, CA
94720, USA.
2Department of Physics, KTH Royal Institute of Technology, SE-106 91
Stockholm, Sweden.
3Department of Physics, St. Kliment Ohridski University of Sofia, 5 James
Bourchier Blvd, 1164 Sofia, Bulgaria.
4Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Str. 38,
01187 Dresden, Germany.
5Chemical Sciences Division Lawrence Berkeley National Laboratory, Berkeley,
CA 94720, USA.

In this supplementary material we provide additional experimental and
theoretical details for the article “Observation of a long-lived prethermal
discrete time crystal created by two-frequency driving”. In Sec. I, we briefly
introduce the Hamiltonian of the system. We continue with motivating the
employed driving scheme in Sec. II, details on the experimental implementation
(Sec. III) and numerical simulations (Sec. IV). Section V provides additional
experimental data backed up by numerics that studies the thermalization
behavior under $\hat{\textbf{z}}$ kicks, and also proves the equivalence of
$\hat{\textbf{z}}$ and $\hat{\textbf{y}}$ slow kicks. Subsequently, in Sec.
VI, we derive the effective Hamiltonian associated with the fast
$\hat{\textbf{x}}$-, and slow $\hat{\textbf{z}}$-concatenated drives. This is
followed in Sec. VII by a detailed analysis of the formation of time-
crystalline order and the prethermalization dynamics in our system. We
conclude in Sec. VIII with a detailed analysis of the heating timescales and
the lifetime of the prethermal order. Throughout the text, we compare our
theoretical predictions against observed experimental data; in particular, we
compare our results in detail in Sec. V, Sec. VII and Sec. VIII. In Sec. IX we
investigate the implications of two-frequency Floquet engineering. Finally,
Sec. X we discuss the random graph dependence of our numerical simulations.

###### Contents

  1. I Introduction
  2. II System
  3. III Prethermalization and discrete time crystals
  4. IV  Melting of prethermal order
  5. V  Two-frequency Floquet engineering
  6. VI Conclusions and Outlook
  7. VII Acknowledgments
  8. VIII Methods
  9. I  Model
  10. II  Motivation for using a two-frequency drive
  11. III  Details of the experimental implementation
    1. III.1 Pulse calibration and parameter regime
  12. IV Details of the numerical simulations
    1. IV.1 Random graph design
    2. IV.2 Time evolution and initial state
  13. V Thermalizing dynamics: simulation versus experiment
  14. VI  Effective Hamiltonian
    1. VI.1 Analysis of the fast $\hat{\textbf{x}}$-drive
    2. VI.2 Analysis of the slow $\hat{\textbf{z}}$-drive.
  15. VII  Prethermal Floquet time crystalline order
    1. VII.1  Discrete spatio-temporal symmetry breaking and prethermal time crystalline order
    2. VII.2  Thermalizing dynamics in the presence of emergent (quasi-)conservation laws
    3. VII.3 Classification of discrete time crystalline order
  16. VIII Heating timescales and duration of the prethermal plateau
    1. VIII.1 Heating timescales associated with the slow $\hat{\textbf{z}}$-drive
    2. VIII.2 Heating timescales associated with the fast $\hat{\textbf{x}}$-drive
    3. VIII.3 Heating model for the concatenated two-frequency drive
    4. VIII.4 Frequency dependence of heating timescales
  17. IX Two-frequency Floquet engineering
  18. X  Random graph independence

##  I  Model

In this section, we analyze the dynamics of the Floquet drive built out of two
superimposed step drives, that we use in the main text. Since the periods of
repetition, $N$ and $M$, of the two step drives may vary independently, we
refer to this protocol as a “two-frequency” drive [cf. Fig. S1]. Note that the
term frequency here refers to the inverse periods of the two superimposed
drives, rather than the support of their Fourier decompositions. Discrete [44,
45] and continuous [46] versions of such two-frequency drives have been
previously used for dynamical decoupling protocols and Hamiltonian engineering
[47].

We investigate a hyperpolarized lattice of ${}^{13}\textrm{C}$ nuclei in
diamond [34, 35]. In the high-field (secular) approximation, the physics of
the interactions in the system is captured by the dipolar spin Hamiltonian,

$\displaystyle\mathcal{H}_{\mathrm{dd}}=\sum_{j<k}^{L}b_{jk}\left(3I_{jz}I_{kz}-\vec{I}_{j}\cdot\vec{I}_{k}\right),$
(1)

where
$b_{jk}=\frac{\mu_{0}}{4\pi}\hbar\gamma_{n}^{2}(3\cos^{2}\alpha_{jk}-1)\frac{1}{|\vec{r}_{jk}|^{3}}$
with the gyromagnetic ratio $\gamma_{n}=10.7\mathrm{MHz/T}$; the angle of the
interspin vector $\vec{r}_{jk}$ and the external magnetic field $\vec{B}_{0}$,
is
$\alpha_{jk}=\cos^{-1}\left(\frac{\vec{r}_{jk}\cdot\vec{B}_{0}}{|\vec{r}_{jk}||\vec{B}_{0}|}\right)$
. $I_{j\mu}$ describes a spin-1/2 operator in direction $\mu$ at some spatial
position $j$. The number of spins in the system is $L$. In the main text, we
had introduced the coupling strength $J=\left\langle b_{jk}\right\rangle$.

${}^{13}\textrm{C}$ spins located near paramagnetic impurities, such as NV-
centers or other lattice paramagnetic defects (e.g., P1 centers), experience a
hyperfine-mediated magnetic field localized from the defect site. Since this
field falls off as $1/r^{3}$, different ${}^{13}\textrm{C}$ atoms experience a
different shift, based on their location and the distance to the nearest
paramagnetic impurity. This effect can be approximately modeled by introducing
an additional term in the Hamiltonian,

$\displaystyle\mathcal{H}_{z}=\sum_{j=1}^{L}c_{j}I_{jz},$ (2)

where we assume that $c_{j}$ are Gaussian distributed random numbers with a
standard-deviation on the order of the relevant energy scales in Eq. (1) [39].

The total system Hamiltonian in the laboratory frame is composed of
$\mathcal{H}_{\mathrm{dd}}$, $\mathcal{H}_{z}$ and the trivial Zeeman
Hamiltonian. Assuming the pulses applied are on-resonance, however, the system
is conveniently described in the rotating frame by the Hamiltonian,

$\mathcal{H}=\mathcal{H}_{\mathrm{dd}}+\mathcal{H}_{z}.$ (3)

Next, we additionally subject the system to a two-step periodic drive created
by two independent external fields applied along (i) the $\hat{\textbf{x}}$
direction, $U_{x}$, and (ii) the $\hat{\textbf{z}}$-direction, $U_{z}$, cf.
Fig S1 [in the theoretical analysis below, we consider $\hat{\textbf{z}}$
kicks, but in experiments $\hat{\textbf{y}}$-kicks are used as they are easier
to implement; as we show in Sec. VI, $\hat{\textbf{z}}$\- and
$\hat{\textbf{y}}$-kicks create the same effective Hamiltonian to leading
order and thus yield the same effect on the dynamics, cf. Fig. S9]. The
physical system, which generates the time evolution $U_{\mathcal{H}}$, is
exposed to a rapid periodic application of the $\hat{\textbf{x}}$-drive; every
$N$ repetitions of the $\hat{\textbf{x}}$-drive, we apply the
$\hat{\textbf{z}}$-field $U_{z}$ once, which completes one Floquet period (or
cycle). We note that the two constituent drives used here have commensurate
periods. The entire evolution over one driving cycle is then described by the
Floquet unitary $U_{F}$:

$\displaystyle U_{F}$ $\displaystyle=$
$\displaystyle\left(U_{x}U_{\mathcal{H}}\right)^{N}U_{z},$ (4) $\displaystyle
U_{z}$ $\displaystyle=$
$\displaystyle\mathrm{e}^{-i\gamma\mathcal{I}_{z}},\quad
U_{x}=\mathrm{e}^{-i\vartheta\mathcal{I}_{x}},\quad
U_{\mathcal{H}}=\mathrm{e}^{-i\tau\mathcal{H}},$

where we defined the net spin operator $\mathcal{I}_{\nu}=\sum_{j}I_{j\nu}$
with $\nu\in\\{x,y,z\\}$. Note here we have assumed that the pulse action can
be described as instantaneous rotations ($\delta$-pulse approximation). The
angles of the $\hat{\textbf{x}}$ and $\hat{\textbf{z}}$ rotations are denoted
$\vartheta$ and $\gamma$, respectively [Fig S1].

The stroboscopic dynamics of any initial state $|\psi_{0}\rangle$ is thus
determined by successive applications of $U_{F}$

$|\psi(M)\rangle=U_{F}^{M}|\psi_{0}\rangle=\left[\left(U_{x}U_{\mathcal{H}}\right)^{N}U_{z}\right]^{M}|\psi_{0}\rangle,$
(5)

where $M$ denotes the stroboscopic Floquet cycle, and similarly for mixed
initial states.

Fig. S1: Schematic of the two-frequency drive. An effective Hamiltonian is
stabilized using a high-frequency drive with rapidly applied pulses along
$\hat{\textbf{x}}$ direction. Low-frequency $\hat{\textbf{z}}$-kicks are
applied after $N$ fast kicks. This cycle is repeated so as to apply a total of
$M$ slow kicks. Ultimately, the system is driven towards infinite temperature,
and we are able to track this spin thermalization dynamics. Fig. S2: Effect of
a single high-frequency drive (experimental data). The figure shows the net
13C NMR signal obtained using a single-frequency drive consisting of the
application of a train of $\vartheta$-pulses following the initial
$\pi$/2-pulse (adapted from Ref. [37]). Here, $t_{\textrm{acq}}$=32$\mu$s and
pulse spacing $\tau$=100$\mu$s. There are strong signal dips at
$\vartheta{\approx}\\{0,\pi,2\pi\\}$ due to rapid decay from evolution under
the dipolar Hamiltonian. This signal decay is avoided in this work through the
use of the two-frequency drive (see Fig. S1).  Fig. S3: Single-frequency
driving PDTC. Experimental data showing PDTC with traditional single-frequency
driving after initializing the spins with a polarization along
$\hat{\textbf{z}}$. Each data point represents the integrated signal
$\left\langle\mathcal{I}_{z}\right\rangle$, readout by multiple-pulse spin-
locking after the application of $M$ pulses with flip angle $\gamma=0.97\pi$
and period
$\tau={31.5\textrm{ms},3.15\textrm{ms},1.05\textrm{ms},0.315\textrm{ms}}$
(A-D). For $\tau=31.5ms$ (A), the DTC response persists for less than 30
Floquet cycles before reaching 1/e signal decay. The optimal PDTC lifetime is
found for $\tau=1.05$ms (C) where the system sees a quick decay during the
first 10 Floquet cycles followed by a long plateau. When driving is too fast,
the Floquet Hamiltonian becomes effectively single-particle [cf. Eq. (VI.1)]
and a split beat pattern emerges, quickly destroying the DTC response (see
(D), $\tau=0.315$ms). These data motivate the use of a two-frequency in our
system which enables the observation of long-lived PDTC order.

##  II  Motivation for using a two-frequency drive

Let us now elucidate upon the experimental and theoretical motivation for
using the two-frequency drive in this work. We begin by noting that the
${}^{13}\textrm{C}$ nuclei are dipolar coupled with the Hamiltonian
$\mathcal{H}_{\textrm{dd}}$ above, and when prepared in an initial state
$\rho_{0}\propto\mathcal{I}_{x}$, undergo dipolar evolution under it. This
results in the rapid free induction decay, with a time constant
$T_{2}^{\ast}\sim 1.5$ms.

Consider now the application of a single frequency drive, consisting of a
train of $\vartheta$-pulses spin locked with the ${}^{13}\textrm{C}$ nuclei.
Figure S2 shows the integrated signal obtained for different values of
$\vartheta$. While the resulting spin decay lifetimes $T_{2}^{\prime}$ can be
far in excess of $T_{2}^{\ast}$, the fastest signal decay occurs for
$\vartheta=\\{0,\pi,2\pi\\}$. This stems from the fact that the dipolar
interaction is invariant under rotations by $\pi$; this yields an average
Hamiltonian under the Floquet drive,
$\bar{\mathcal{H}}^{(0)}=\mathcal{H}_{\textrm{dd}}$. For other angles of
$\vartheta$, it is easy to show instead that the dipolar interaction is
effectively engineered to a form that commutes with the initial state (see
Sec. VI). This fast decay at $\vartheta=\pi$ would normally prevent
observation of a PDTC in the rotating frame with a single frequency drive
alone. In our work, this problem is essentially circumvented by the use of the
two frequency drive.

We note that it is still possible to observe the DTC-like behavior in the
laboratory frame with a single frequency drive. This is shown experimentally
in Fig. S3. However, obtaining data similar to Fig. 2A of the main text then
has to be constructed in a point-by-point fashion, and the experiment has to
be re-initialized between successive measurement kicks. This makes observation
of the formation and melting of the PDTC order at high resolution extremely
challenging (cf. Fig. 3 of main text). In summary, our two-frequency drive
approach allows the ability to excite the PDTC order while it simultaneously
permits the tracking of the spin thermalization continuously, yielding a
unique means to observe the driven interacting spin dynamics away from
equilibrium.

From the theoretical perspective, multi-frequency drives are especially
interesting as they allow to Floquet engineer novel kinds of effective
Hamiltonians and non-equilibrium ordered states that might be inaccessible
with single-frequency drives. To illustrate this idea, let us assume a two-
step (i.e., single-frequency) drive composed of the Hamiltonians $H_{1}$ and
$H_{2}$ repeatedly applied with amplitudes $J$ and $\vartheta$ such that time
evolution over one Floquet period is given by

$\displaystyle U_{F}=\exp\left(-iJH_{1}\right)\exp\left(-i\vartheta
H_{2}\right).$ (6)

Provided $||JH_{1}||\sim||\vartheta H_{2}||\ll 1$, the Baker-Campbell-
Hausdorff expansion predicts a Floquet Hamiltonian of the form
$H_{F}=(JH_{1}+\vartheta H_{2})/(J+\vartheta)+\mathcal{O}(T^{\prime})$ with
$T^{\prime}\sim J,\vartheta$. Although, $H_{F}$ might provide a new effective
Hamiltonian that mixes the properties of $H_{1}$ and $H_{2}$, there is only a
limited amount of flexibility: the lowest order effective Hamiltonian will
always be given by the average Hamiltonian of $JH_{1}$ and $\vartheta H_{2}$.
In contrast, a much more interesting regime is found when $||JH_{1}||\ll 1$
but $||\vartheta H_{2}||\sim 1$. In this case, the standard Baker-Campbell-
Hausdorff expansion leads to an infinite series of terms $\propto J$ (i.e.,
terms $\propto J\vartheta,J\vartheta^{2},J\vartheta^{3}$, etc.), and more
sophisticated techniques are required to find closed form expressions [48].
Eventually, the series resummation leads to new types of (interaction) terms
which can depend on $\vartheta$ in a non-linear way (see also Sec. VI.1). This
ansatz allows us to engineer new effective Hamiltonians with case-specific
properties beyond the average Hamiltonian: for instance, in this work, we
engineer a many-body effective Hamiltonian with an emergent $U(1)$ and
$\mathbb{Z}_{2}$ symmetry, which is subject to a second drive to implement
non-equilibrium order – in our case a prethermal discrete time-crystalline
order – in the system.

Fig. S4: Calibration of flip angle. Experimental data points represent
integrated signal readout via spin-locking following initial excitation pulse
of finite width. Experiments were performed for different pulse widths and
plot represents a Rabi oscillation, with flip angle $\pi$ being calibrated to
the zero-crossing and the rest of the flip angles scaling linearly with pulse
width.

##  III  Details of the experimental implementation

###  III.1 Pulse calibration and parameter regime

The flip angles $(\vartheta,\gamma)$ in the experiments are calibrated using a
${}^{13}\textrm{C}$ nuclear Rabi oscillation. The first tipping pulse (to be
calibrated) is followed by a $\vartheta{\approx}\pi/2$ spin-lock train. We
plot the integrated signal in Fig. S4. The SNR is very good because of the
long resulting rotating frame lifetimes $T_{2}^{\prime}$. Furthermore, in
order to make the experiments less susceptible to precise calibration of the
first $\pi/2$ pulse (that tips the spins onto the
$\hat{\textbf{x}}{\text{-}}\hat{\textbf{y}}$ plane), in Fig. 1B [main text] we
apply the first phase $\gamma$-kick to the spins after a period of 1s, so that
the spins have already prethermalized along $\hat{\textbf{x}}$ before the
$\gamma$-kicks are applied. This ensures that the initial state of the spins
in the experiment is exactly $\rho_{0}\propto\epsilon\mathcal{I}_{x}$, with
the experimentally measured value $\epsilon=0.68\%$. Fig. S5 shows
experimental data for the system’s natural free induction decay time and
maximum possible lifetime extension under Floquet control. Further relevant
system and drive parameters are collected in Tab. S1.

Fig. S5: Free induction decay and pulsed spin-lock data. Measurements shown for free induction decay (green line and points) and for pulsed spin-locking Floquet control (red line and points). Free induction decay lifetime $T_{2}^{*}=1.5$ms is extended with high-frequency pulsing by over 5 orders of magnitude, $T_{2}^{\prime}=90.9$s. Adapated from Ref. [37].  Parameter | Value
---|---
Inter-spin spacing | $1\mathrm{nm}$
Initial state $\hat{\textbf{x}}$-polarization $\epsilon$ | $0.68\%$
Median interaction strength $\langle J\rangle$ | $0.66\mathrm{kHz}$
Drive strength $\Omega$ | $7.6\mathrm{kHz}$
Total spin lifetime $T_{2}^{\prime}$ | $90.9\mathrm{s}$
Normalized interaction strength $\langle J\rangle/\gamma_{n}^{2}$ | $5.8\times 10^{-12}\mathrm{HzT^{2}}$
Normalized drive strength $\Omega/\gamma_{n}$ | $7.1\times 10^{-}4\mathrm{T}$
Normalized spin lifetime $\gamma_{n}T_{2}^{\prime}$ | $9.7\times 10^{-8}\mathrm{T^{-1}}$
Sampling rate | $(0.105\mathrm{ms})^{-1}$
Number of measured points per kick | 299
Total observation time | $14\mathrm{s}$
Signal to noise ratio per point | $>10^{2}$
Table S1: Summary of experimentally relevant system and drive parameters.

##  IV Details of the numerical simulations

Fig. S6: Schematic for the generation of a pseudo-random graph used in
numerical simulations: Starting off with a spin $\vec{I}_{1}$ in some random
position, we propose a new position for the next spin $\vec{I}_{2}$ [see
text]. Since $r_{\mathrm{min}}<r_{12}<r_{\mathrm{max}}$ we accept
$\vec{I}_{2}$ and add it to the graph. Next, we propose $\vec{I}_{3}$, which
also satisfies all requirements since all inter-spin vectors
$r_{13},r_{23}>r_{\mathrm{min}}$, while $r_{13}<r_{\mathrm{max}}$. Hence, we
also accept $\vec{I}_{3}$. The next proposal, $\vec{I}_{4}$, only has inter-
spin vectors which are all larger than $r_{\mathrm{max}}$ and, thus,
constitutes an invalid choice; hence, it is discarded. If we aim to simulate
$L=4$ spins (or more) we would have to continue proposing new spins until a
valid fourth spin position is found and added to the graph.

###  IV.1 Random graph design

We perform exact numerical simulations of the system with $L$
${}^{13}\textrm{C}$ nuclear spins, placed on a pseudo-random graph. In
creating the lattice graph, we respect two conditions: (i) we ensure a certain
average spin density at the system size $L$; (ii) we avoid isolated spin
positions (‘outliers’), that would not contribute to the many-body dynamics as
their coupling constant becomes negligible in comparison to their “bulk” spins
counterparts due to the spatial decay of the dipolar interactions.

To iteratively generate the spin positions on the graph, we apply the
following rule: first, we randomly generate a new spin position, and we check
if it has a minimum distance $r_{\mathrm{min}}$ to all other spins. In
addition, we also require that at least one of the mutual interspin vectors
satisfies $|\vec{r}_{jk}|<r_{\mathrm{max}}$ [cf. Fig. S6]. If these conditions
are met, we accept the new spin position and add it to the list of spins on
the graph. Otherwise we discard the proposed spin position and start over with
a new proposal. In the main text we use $r_{\mathrm{min}}=0.7$,
$r_{\mathrm{max}}=0.8$, both in units of
$\sqrt[3]{\mu_{0}\hbar\gamma_{n}^{2}}$.

###  IV.2 Time evolution and initial state

In the numerical simulations using this pseudo-random graph of spins, unless
stated explicitly otherwise, we initialize the system in the
$\hat{\textbf{x}}$-polarized pure product state
$|\psi_{0}\rangle\\!=\\!\bigotimes_{j=1}^{L}\frac{1}{\sqrt{2}}\left(|\\!\uparrow_{j}\rangle\\!+\\!|\downarrow_{j}\rangle\right)$,
and perform numerically exact time evolution according to the protocol given
in Eq. (4). While most of our results are based on this initial state, we also
show that the effects we observe are qualitatively initial-state independent
[cf. Sec. VII.2 and Fig. S10].

Since the interaction is long-range and the positions of the nuclei are not
spatially ordered, the characteristic energy scale of the model is not
immediately obvious, and needs to be extracted from the dynamics. To this end,
we define a characteristic energy $J=1/\tau_{d}$, where $\tau_{d}$ is the
timescale on which single-particle observables decay by a factor of
$1/\mathrm{e}$, as they approach their equilibrium value when evolved under
$U_{\mathcal{H}}$. We note in passing that $J$ may vary depending on the
initial state, i.e., different values might be obtained for $J$ when probed
with other initial states than $|\psi_{0}\rangle$. Therefore, the value of $J$
only serves as a rough estimate of relevant energy scales.

To further increase the ergodicity of the drive and diminish finite-size
effects in the dynamics, we add a small uniformly distributed random “noise”
$\delta\tau\in[-0.05\tau,0.05\tau]$ to the duration $\tau$ for which we apply
$\mathcal{H}$, i.e. $\tau\rightarrow\tilde{\tau}=\tau+\delta\tau$ so that
$\tilde{\tau}$ is slightly different in each cycle of $(U_{x}U_{\mathcal{H}})$
[49, 50]. After each application of $U_{x}$, $U_{z}$, respectively, we compute
the expectation value of single particle observables

$\langle x\rangle=\frac{2}{L}\langle\mathcal{I}_{x}\rangle,\quad\langle
y\rangle=\frac{2}{L}\langle\mathcal{I}_{y}\rangle,\quad\mathrm{and}\quad\langle
z\rangle=\frac{2}{L}\langle\mathcal{I}_{z}\rangle.$

##  V Thermalizing dynamics: simulation versus experiment

Fig. S7: Experimental data showing a linear scale depiction of Fig. 3A of the
main text. 285 traces similar to Fig. 2B of main paper are plotted stacked for
different values of kick angle $\gamma$ in $[-\pi,\pi]$. Colors represent
signal $\left\langle\mathcal{I}_{x}\right\rangle$ (see colorbar). Time $N\tau$
and Floquet cycle number $M$ run vertically and are plotted here on a linear
scale [cf. left and right $\hat{\textbf{x}}$-axes]. Insets show PDTC behavior
(evidenced by oscillations between positive and negative signals every 300
fast pulses) in 0.16$\pi$\- radian-wide windows near $-\pi$ (i) and $\pi$
(ii). The PDTC response is strong at early times (i) and persists to long
times (ii) despite signal decay caused by heating. See ref. [51, 33, 52] for
movies of data for first 25, 55 and 155 Floquet cycles, respectively.  Fig.
S8: Comparison of numerical simulation and experiment –(a)-(b) Time evolution
of the $\hat{\textbf{x}}$-polarized initial state under the drive generated by
Eq. (4) at $\gamma=\varepsilon$ (a) and $\gamma=\pi+\varepsilon$ (b) for
different values of $\varepsilon$: we plot the expectation value of the
$\hat{\textbf{x}}$-magnetization, $\langle
x\rangle=\frac{2}{L}\langle\mathcal{I}_{x}\rangle$, vs. time given in number
of kicks. (c) $|\langle x\rangle|$ as a function of the
$\hat{\textbf{z}}$-kick angle $\gamma$ and number of kicks. To generate the
random graph of spins we used $r_{\mathrm{min}}=0.7$, $r_{\mathrm{max}}=0.8$
(in units of $\sqrt[3]{\mu_{0}\hbar\gamma_{n}^{2}}$) with $L=14$. The single
particle energies $c_{j}$ are normally distributed random numbers with mean
$\langle c_{j}\rangle=\overline{b_{jk}}$ and standard deviation
$\sigma_{c}=10\times\overline{b_{jk}}$, where $\overline{b_{jk}}$ is the
median of all coupling constants of the graph. (d)-(f) experimental data
corresponding to (a) and (b) and (c) respectively. Parameters are as in Fig. 3
of the main text. In all panels, we plot logarithmically many points in the
number of kicks. $N=300$ and $\tau J=0.07$ for all figures.  Fig. S9:
Equivalence of $\hat{\textbf{z}}$ and $\hat{\textbf{y}}$ kicks. Experimental
results analogous to the dataset shown in Fig. S8, except here the low-
frequency driving consists of composite pulses along $\hat{\textbf{z}}$-axis
rather than simple $\hat{\textbf{y}}$-pulses (see Fig. S1 for driving
protocol). NMR signal $\left\langle\mathcal{I}_{x}\right\rangle$ is measured
between fast pulses and composite $\hat{\textbf{z}}$-rotation is applied
between every $N=176$ fast pulses. Single-shot experiment yields one line, and
experiment is repeated for different values of $\gamma$ in $[0,\pi]$. Just as
with the experiment with $\hat{\textbf{y}}$-pulses, the dynamics show long
lifetimes for small flip angle deviations $\varepsilon$, $\gamma=\varepsilon$
(a) and $\gamma=\pi+\varepsilon$ (b).

We now elucidate in greater detail the thermalization dynamics of the spins in
a comparison between theory and experiment. First, we begin by displaying Fig.
S7 that shows the data corresponding to the experiments in Fig. 3B of the main
text, but with the time-axis plotted in a linear (as opposed to a logarithmic)
scale. We plot the survival probability of the spins
$\left\langle\mathcal{I}_{x}\right\rangle$ under the concatenated action of
the slow and fast drives. The stable points around $\gamma{=}0$ and $\pm\pi$
are clearly visible, along with the large region in the center where the spins
decay rapidly. The insets (reproduced from the main text), show the PDTC
behavior in more detail for two ${\approx}1$s-long windows centered at
$t{=}3.15$s and $t{=}11.55$s.

Figure S8 displays the dynamics of the system, obtained numerically ((a)-(c))
and in the experiment ((d)-(f)). In the numerical simulations, the initial
state $|\psi_{0}\rangle$ is evolved in time up to $5\times 10^{4}$ kicks for
$N=300$ [throughout the supplementary material, we also show results for a few
more values of $N$]. Note that here we treat $\hat{\textbf{x}}$-kicks and
$\hat{\textbf{z}}$-kicks on equal footing: each application of $U_{x}$,
$U_{z}$, respectively, adds $+1$ to the total number of kicks. We show the
simulated time evolution curves for $\gamma=0+\varepsilon$ (a) and
$\gamma=\pi-\varepsilon$ (b), where different curves correspond to different
values of $\varepsilon$ (see colorbar): with an increasing value of
$\varepsilon$ the heating time gets gradually reduced until no stable
magnetization is visible anymore when $\varepsilon\to\pi/2$.

In the experiment, the initial state after the $\pi/2$-pulse about the
$y$-axis [see Fig. 1, main text] is the mixed state
$\rho_{0}\propto\epsilon\mathcal{I}_{x}$. The time evolution of the
$\hat{\textbf{x}}$-magnetization shown Fig. S8 (d) and (e) is measured in the
experiment. Despite the different initial states, the curves are thoroughly
comparable on a qualitative basis, with the notable difference that the
experimental curves yield much cleaner results in the small $\varepsilon$
regime. This can be traced back to large differences in the system size in
experiment ($L\sim\mathcal{O}(10^{4})$) and theory
($L\sim\mathcal{O}(10^{1})$).

In Fig. S8 (c) (Fig. S8 (f)) we display theoretical (experimental) results of
the absolute value of the $\hat{\textbf{x}}$-magnetization, $|\langle
x\rangle|$, for a wide range of $\gamma$-values. In both panels, we can
distinguish between three overall regimes: around $\gamma=0$ ($\gamma=2\pi$)
heating is suppressed as $U_{z}$ almost wraps up to the identity so that the
slow-frequency $\hat{\textbf{z}}$-kicks become ineffective, and the effective
Hamiltonian (discussed in Sec. VI) is engineered solely from the high-
frequency $\hat{\textbf{x}}$-drive. Around $\gamma=\pi$, we observe a similar
behaviour; there, the formation of the PDTC [cf. Sec. VII.1] leads to slow
thermalization dynamics. In between these regimes the $\hat{\textbf{z}}$-kicks
lead to rapid thermalization to infinite temperature resulting in a decay of
the signal within very few cycles (cf. Fig S7).

We note finally that quantitatively identical results are obtained using a
$\hat{\textbf{z}}$-kick to the spins instead of a $\hat{\textbf{y}}$-kick. In
the experiment, however, this requires each kick to be constructed out of a
composite rotations along the $\hat{\textbf{x}}$ and $\hat{\textbf{y}}$ axes.
Data corresponding to such $\hat{\textbf{z}}$-kicks are shown in Fig. S9
around the two stable points, $\gamma=\epsilon$ and $\gamma=\pi-\epsilon$.

Overall, the excellent qualitative agreement between theory and experiment is
remarkable considering the three orders of magnitude difference in the system
size between the two. From a theory perspective, this justifies the results
obtained using exact simulation of the dynamics in the modelled system, and
indicates that these numerical techniques can be used to make reliable
predictions about experimental systems. Nonetheless, the experiment covers
system sizes that are infeasible to reach on any classical computer in the
foreseeable future. Therefore, the experiment remains indispensable for
probing statistical mechanics concepts related to collective phenomena, such
as thermalization dynamics or symmetry breaking, which are only well-defined
in the thermodynamic limit.

##  VI  Effective Hamiltonian

Intuitively, the physics of the two-frequency drive is easy to comprehend when
the repetition rates of the two different drives are well separated from one
another. In this case, the repetition number of the $\hat{\textbf{x}}$-drive
is $N\gg 1$ [cf. Fig. S1]. In this regime, we can consider the
$\hat{\textbf{x}}$-drive to be fast, and the $\hat{\textbf{z}}$-drive – slow,
w.r.t. the energy scale $J$ of the Hamiltonian.

For the subsequent analysis, it will prove useful to rewrite Eq. (4) in the
so-called toggling frame:

$\left(U_{x}U_{\mathcal{H}}\right)^{N}U_{z}=\left(U_{x}U_{\mathcal{H}}U_{x}^{-1}\right)\left(U_{x}^{2}U_{\mathcal{H}}U_{x}^{-2}\right)\dots\times
U_{x}^{N}U_{z}.$ (7)

Denoting the toggling frame Hamiltonians by
$\mathcal{H}_{n}=U_{x}^{n}\mathcal{H}U_{x}^{-n}$ leads to the compact
expression

$U_{F}=\left(U_{x}U_{\mathcal{H}}\right)^{N}U_{z}=\left(\prod_{n=1}^{N}\exp(-i\tau\mathcal{H}_{n})\right)U_{x}^{N}U_{z}.$
(8)

The toggling-frame expansion is based on the Magnus expansion [53, 54, 55] in
$N^{-1}$, and provides a tool to compute an approximation for the generator of
the $N$-cycle unitary $U^{N}=\left(U_{x}U_{\mathcal{H}}\right)^{N}$.

Apriori, the toggling-frame expansion aims to directly approximate the
generator of $U^{N}$ after $N$ cycles, which is different from the
$\hat{\textbf{x}}$-drive Floquet unitary $U=U_{x}U_{\mathcal{H}}$ associated
with the repetition of the $\hat{\textbf{x}}$-drive only. With this, we draw
the readers’ attention to the subtle distinction between the Magnus expansion
(which does not require a periodic drive), and the Floquet-Magnus expansion
(which is the specialization of the former to Floquet systems). The toggling-
frame expansion is based on the Magnus expansion.

Upon adding the second drive, the proper stroboscopic Floquet period is
increased. Because of this, rather coincidentally, the toggling frame
expansion w.r.t. $N$ will actually approximate the proper Floquet unitary
$U_{F}=\left(U_{x}U_{\mathcal{H}}\right)^{N}U_{z}$, up to the overall single-
particle initial kick $U_{x}^{N}U_{z}$, cf. Eq. (8) and the subsequent
analysis in Sec. VI.1. Then, one can consider applying the toggling frame
expansion once more, this time w.r.t. the $\hat{\textbf{z}}$-drive repetition
rate $M$, cf. Eq. (5). This directly provides us an approximate series
expansion in powers of $M^{-1}$ for the generator of $[U_{F}]^{M}$, [Sec.
VI.2].

###  VI.1 Analysis of the fast $\hat{\textbf{x}}$-drive

Toggling-Frame Expansion – Let us assume the $\hat{\textbf{x}}$-kicks to be
high-frequency compared to the energy scales of the system, i.e., $\tau J\ll
1$. Then, applying the Baker-Campbell-Hausdorff (BCH) formula to leading-
order, the toggling frame unitary can be written as

$\displaystyle\prod_{n=1}^{N}\exp(-i\tau\mathcal{H}_{n})=\exp\left[-iN\tau\bar{\mathcal{H}}^{(0)}+\mathcal{O}\left((\tau
J)^{2}\right)\right],$ (9)

where the effective Hamiltonian $\bar{\mathcal{H}}^{(0)}$ to order
$\mathcal{O}\left((\tau J)^{2}\right)$, is given by [56]

$\displaystyle\bar{\mathcal{H}}^{(0)}=\frac{1}{N}\sum_{n=1}^{N}\mathcal{H}_{n}=\bar{\mathcal{H}}^{(0)}_{\mathrm{dd}}+\bar{\mathcal{H}}^{(0)}_{z},$
(10)

where

$\displaystyle\bar{\mathcal{H}}^{(0)}_{\mathrm{dd}}$ $\displaystyle=$
$\displaystyle\sum_{j<k}b_{jk}\bigg{(}\frac{3}{2}\bigg{[}\mathcal{H}_{\mathrm{ff}}+\mathcal{G}_{\mathrm{c}}(N,\vartheta)\mathcal{H}_{\mathrm{dq}}-\tilde{\mathcal{H}}_{\mathrm{ff}}\mathcal{G}_{\mathrm{s}}(N,\vartheta)\bigg{]}-\vec{I}_{j}\vec{I}_{k}\bigg{)},$
$\displaystyle\bar{\mathcal{H}}^{(0)}_{z}$ $\displaystyle=$
$\displaystyle\sum_{j}c_{j}\left(\mathcal{G}_{\mathrm{c}}(N,\vartheta/2)I_{jz}-\mathcal{G}_{\mathrm{s}}(N,\vartheta/2)I_{jy}\right),$
(11)

with

$\displaystyle\mathcal{H}_{\mathrm{ff/dq}}$ $\displaystyle=$ $\displaystyle
I_{jz}I_{kz}\pm I_{jy}I_{ky},$
$\displaystyle\tilde{\mathcal{H}}_{\mathrm{ff}}$ $\displaystyle=$
$\displaystyle I_{jz}I_{ky}+I_{jy}I_{kz},$ (12)

and

$\displaystyle\mathcal{G}_{\mathrm{c}}(N,\vartheta)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\frac{\sin(N\vartheta)}{\sin(\vartheta)}\cos\left((N+1)\vartheta\right),$
(13) $\displaystyle\mathcal{G}_{\mathrm{s}}(N,\vartheta)$ $\displaystyle=$
$\displaystyle\frac{1}{N}\frac{\sin(N\vartheta)}{\sin(\vartheta)}\sin\left((N+1)\vartheta\right).$

Notice that $\mathcal{G}_{\mathrm{c/s}}(N,\vartheta)$ decays as $1/N$ for
large $N$. Thus, as long as $\vartheta\not\approx l\pi$ ($l\in\mathbb{N}$), we
find $\mathcal{G}_{\mathrm{c/s}}(N,\vartheta)\rightarrow 0$ in the $N\gg 1$
regime, which simplifies to,

$\displaystyle\bar{\mathcal{H}}^{(0)}\simeq\overline{\mathcal{H}}=\sum_{j<k}b_{jk}\bigg{(}\frac{3}{2}\mathcal{H}_{\mathrm{ff}}-\vec{I}_{j}\cdot\vec{I}_{k}\bigg{)}.$
(14)

Importantly, we find $\left[\mathcal{I}_{x},\overline{\mathcal{H}}\right]=0$,
which implies that the drive preserves any polarization in $\hat{\textbf{x}}$
direction to leading order in $\tau J$. One can readily convince oneself that
the smaller higher-order terms break this emergent conservation law. However,
at least to order $\mathcal{O}((J\tau)^{2})$, these terms are additionally
suppressed by a factor of $1/N$ [56]. Thus, we expect the leading order
effective Hamiltonian to capture well all the interesting physics that plays
out in the prethermal state.

Note that the above analysis breaks down as $\vartheta\rightarrow l\pi$.
Indeed, in this limit we lose the quasi-conservation of the
$\hat{\textbf{x}}$-magnetization already in the leading-order Hamiltonian (cf.
Fig. S2). This implies that initially polarized states will decay rapidly to
their equilibrium value close to zero magnetization. In the present
experiment, the state of the system would then appear indistinguishable from
an infinite-temperature state, since we can only measure single-particle
operators, whose expectation values match the infinite-temperature value in
the quench dynamics governed by $\overline{\mathcal{H}}$. Fortunately, this
behavior can be prevented by setting $\vartheta$ sufficiently far away from
the $0,\pi,2\pi$ by employing a two-frequency drive. In this work, we
therefore work at $\vartheta=\pi/2$.

Replica Expansion – An alternative approach to the dynamics generated by the
fast $\hat{\textbf{x}}$-drive is given by the Floquet Hamiltonian
$\mathcal{H}_{F,x}$. Unlike the effective Hamiltonian, $\mathcal{H}_{F,x}$
only governs the dynamics over a single period of $\hat{\textbf{x}}$ driving,
$U_{x}U_{\mathcal{H}}$. Formally, $\mathcal{H}_{F,x}$ is defined by,

$U_{x}U_{\mathcal{H}}=\mathrm{e}^{-i(\vartheta+\tau)\mathcal{H}_{F,x}}=\mathrm{e}^{-i\vartheta\mathcal{I}_{x}}\mathrm{e}^{-i\tau\mathcal{H}}.$
(15)

Typically, $\mathcal{H}_{F,x}$ can be obtained via a BCH expansion. However,
having $\tau J\ll\vartheta$ complicates the standard procedure. In fact, here,
the BCH expansion produces infinitely many terms $\propto\tau\vartheta$. A
closed form expression requires a resummation of these terms. A neat way how
this can be achieved was recently put forward in Ref. [48] by using a replica
resummation. This trick utilizes a replica expression of the logarithm

$\displaystyle(\vartheta+\tau)\mathcal{H}_{F,x}=i\log(U_{x}U_{\mathcal{H}})=\lim_{\rho\rightarrow
0}\frac{i}{\rho}\left((U_{x}U_{\mathcal{H}})^{\rho}-1\right).$ (16)

In the next step, the right hand side of Eq. (16) can be expanded in powers of
$\tau$. Eventually, to linear order in $\tau$ this yields

$\displaystyle(\vartheta+\tau)\mathcal{H}_{F,x}$ $\displaystyle=$
$\displaystyle\vartheta\mathcal{I}_{x}+\tau\bigg{[}\sum_{j<k}b_{jk}\bigg{(}\frac{1}{2}(1\\!-\\!3\vartheta\cot(\vartheta))I_{jy}I_{ky}$
(17) $\displaystyle+$
$\displaystyle\frac{1}{2}(1\\!+\\!3\vartheta\cot(\vartheta))I_{jz}I_{kz}\\!-\\!I_{jx}I_{kx}\\!-\\!\frac{3}{2}\vartheta(I_{jz}I_{ky}\\!+\\!I_{jy}I_{kz})\bigg{)}$
$\displaystyle-$
$\displaystyle\sum_{j}c_{j}\left(\frac{\vartheta}{2}I_{jy}\\!+\\!\frac{\vartheta}{2}\cot(\vartheta/2)I_{jz}\right)\bigg{]}+\mathcal{O}(\tau^{2}).$

At $\vartheta=\pi/2$, the above expression simplifies to

$\displaystyle\left(\frac{\pi}{2}+\tau\right)\mathcal{H}_{F,x}$
$\displaystyle=$
$\displaystyle\frac{\pi}{2}\mathcal{I}_{x}+\tau\overline{\mathcal{H}}-\tau\bigg{(}\frac{3}{4}\pi\sum_{j<k}b_{jk}\tilde{\mathcal{H}}_{\mathrm{ff}}+\frac{\pi}{4}\sum_{j}c_{j}(I_{jy}-I_{jz})\bigg{)}$

Equation (VI.1) shines a more microscopic light on the dynamics generated by
the fast drive. With this approach, we even obtain the expansion up to first
order in $\tau$. Notice that these corrections are all linear in $\tau$,
however they can be highly non-linear in $\vartheta$ as they emerge from
resummations of commutators. The appearance of such terms can be traced back
to the particular structure of the drive where $J\tau\ll 1$ while
$\vartheta\sim\mathcal{O}(1)$. Moreover, Eq. (VI.1) allows us to study the
limit $\tau\rightarrow 0$. There, the system effectively becomes a single
particle problem where no notion of (pre-)thermalization exists. Likewise,
PDTCs become ill-defined since they require many-body interactions [cf. Sec.
VII.1]. This renders the $\tau\rightarrow 0$ regime unfavorable for us as it
decreases the regime of rigidity [see Fig. S3 (D)]. On the other hand, the
large-$\tau$ regime is equally bad, since the onset of higher-order terms
spoils the quasi-conservation of $\hat{\textbf{x}}$-magnetization which
protects the PDTC order, cf. Sec. VII. Hence, eventually we are left with a
window of suitable $\tau$ values so that the many-body nature of the effective
Hamiltonian is guaranteed while simultaneously the
$\hat{\textbf{x}}$-magnetization remains quasi-conserved.

###  VI.2 Analysis of the slow $\hat{\textbf{z}}$-drive.

To sum up, in the $N\gg 1$ regime, we can use the toggling-frame expansion to
approximate Eq. (8) to order $\mathcal{O}((J\tau)^{2})$ as

$\displaystyle
U_{F}\approx\mathrm{e}^{-iT\overline{\mathcal{H}}}U_{x}^{N}U_{z}=U_{0}U_{x}^{N}U_{z},\quad
U_{0}=\mathrm{e}^{-iT\overline{\mathcal{H}}}$ (19)

with $T=N\tau$. Describing the effects of the slow $\hat{\textbf{z}}$-drive
requires a careful analysis of the induced modification to the effective
Hamiltonian, which is a prerequisite to understand the thermalizing dynamics
of the system in the experiment. Therefore, let us now investigate the
repeated application of the slow $\hat{\textbf{z}}$-kicks.

The unitary which describes the time evolution up to $M$ Floquet cycles is
given by

$\displaystyle[U_{F}]^{M}=\left(U_{0}U_{x}^{N}U_{z}\right)^{M}.$ (20)

For large values $N\gg 1$, we have $\left[U_{0},U_{x}\right]=0$, where the
unitary $U_{0}$ is defined in Eq. (19). Moreover, $U_{x}^{N}$ denotes a
rotation by $N\vartheta=N\pi/2$. Thus, we have $U_{x}^{N}\equiv
U_{x}^{\left(N\;\mathrm{mod}\;8\right)}$. Hence, in the limit $N\gg 1$, the
action of the high-frequency $\hat{\textbf{x}}$-drive can be summarized into a
single unitary

$\displaystyle U_{0}U_{x}^{N}\\!$ $\displaystyle\equiv$
$\displaystyle\\!U_{0x}=\exp\bigg{[}-iT\overline{\mathcal{H}}-\frac{i\pi}{2}\left(N\;\mathrm{mod}\;8\right)\mathcal{I}_{x}\\!\bigg{]}.$
(21)

To derive the leading-order effective Hamiltonian of the combined
$\hat{\textbf{x}}$\- and $\hat{\textbf{z}}$-drive, we apply the toggling frame
expansion once again, this time w.r.t. $M$. To do so, we multiply Eq. (20)
from the left by $\mathds{1}=U_{z}^{M}U_{z}^{-M}$, and re-group the terms:

$\displaystyle\left(U_{0}U_{x}^{N}U_{z}\right)^{M}$ $\displaystyle=$
$\displaystyle U_{z}^{M}\left(U_{z}^{-M}U_{0x}U_{z}^{M}\right)\times$ (22)
$\displaystyle\times\left(U_{z}^{-(M-1)}U_{0x}U_{z}^{M-1}\right)\times\dots\times\left(U_{z}^{-1}U_{0x}U_{z}\right)$
$\displaystyle=$ $\displaystyle
U_{z}^{M}\left(\prod_{m=M}^{1}U_{z}^{-m}U_{0x}U_{z}^{m}\right),$

where the $M$-toggling frame Hamiltonians read as

$\displaystyle\mathcal{H}_{m}(\gamma)$ $\displaystyle=$ $\displaystyle
U_{z}^{-m}\left(\overline{\mathcal{H}}\\!+\\!\frac{\pi}{2T}(N\leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ 8)\mathcal{I}_{x}\right)U_{z}^{m}$
$\displaystyle=$
$\displaystyle\sum_{j<k}b_{jk}\bigg{(}\frac{3}{2}I_{jz}I_{kz}\\!+\\!\frac{3}{2}\bigg{[}\cos^{2}(m\gamma)I_{jy}I_{ky}\\!+\\!\sin(m\gamma)^{2}I_{jx}I_{kx}$
$\displaystyle-\sin(m\gamma)\cos(m\gamma)\left(I_{jx}I_{ky}+I_{jy}I_{kx}\right)\bigg{]}-\vec{I}_{j}\vec{I}_{k}\bigg{)}$
$\displaystyle+\frac{\pi\leavevmode\nobreak\ (N\leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\
8)}{2T}\left(\cos(m\gamma)\mathcal{I}_{jx}-\sin(m\gamma)\mathcal{I}_{jy}\right).$

Note that we cannot simply follow similar steps as in the analysis of the
high-frequency $\hat{\textbf{x}}$-drive: since $TJ>1$, the lowest order Magnus
expansion of the toggling frame Hamiltonians is a poor approximation to the
total effective Hamiltonian that governs the time evolution. However, there is
one notable exception: at $\gamma=l\pi$ all higher orders of the Magnus
expansion w.r.t. $M$ are zero and the lowest order becomes exact (up to
corrections from the fast $\hat{\textbf{x}}$-drive). This can be seen right
away from Eq. (VI.2) since
$\mathcal{H}_{m}(l\pi)=\overline{\mathcal{H}}+(-1)^{ml}\pi(N\leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ 8)/(2T)\mathcal{I}_{x}$. Then, the effective
unitary describing evolution over $M$ periods becomes (up to corrections from
the fast drive),

$\displaystyle[U_{F}]^{M}\big{|}_{\gamma=l\pi}$ $\displaystyle\approx$
$\displaystyle
U_{z}^{M}\exp\big{[}-iMT\overline{\mathcal{H}}+d_{N}(M,l)\mathcal{I}_{x}\big{]},$
$\displaystyle d_{N}(M,l)$ $\displaystyle=$
$\displaystyle\frac{\pi(-1)^{l}(\\!(\\!-\\!1)^{lM}\\!-\\!1)(N\leavevmode\nobreak\
\mathrm{mod}\leavevmode\nobreak\ 8)}{2((-1)^{l}-1)}.$ (24)

Away from $\gamma=l\pi$, the Magnus expansion cannot be used to obtain
generically reliable results.

##  VII  Prethermal Floquet time crystalline order

For $l$ even and $\gamma=l\pi$, the $\hat{\textbf{z}}$-kicks are inactive as
they give each spin a complete rotation. For these values of $\gamma$, the
effective Hamiltonian is time-independent (up to corrections
$\mathcal{O}(1/N)$ emerging from the fast-drive, which we neglect). We leave
this case aside for the time being, and will return to study its heating rates
in Sec. VIII.

###  VII.1  Discrete spatio-temporal symmetry breaking and prethermal time
crystalline order

Much of the intriguing behavior of the system is for odd $l$, where it
exhibits prethermal time-crystalline order [3] [in this case, the effective
Hamiltonian depends on the parity of $M$: for odd $M$, we can obtain a
contribution from the second term of Eq. (VI.2); for even $M$, the terms in
the second line in Eq. (VI.2) vanish]. In particular, at kick angles
$\gamma\approx\pi$, we argue that the driven-system satisfies all
prerequisites for prethermal discrete time-crystals:

(i) Spatio-temporal Symmetry Breaking: After two Floquet cycles ($M\\!=\\!2$),
the time-evolution operator Eq. (VI.2) takes a particularly simple form

$\displaystyle[U_{F}]^{2}\big{|}_{\gamma=\pi}=\mathrm{e}^{-i2T\overline{\mathcal{H}}}.$
(25)

Notice that $[\overline{\mathcal{H}},\hat{P}_{z}]=0$, where
$\hat{P}_{z}=\mathrm{e}^{-i\pi\mathcal{I}_{z}}$ and thus the two-cycle
effective Hamiltonian $\overline{\mathcal{H}}$ obeys a discrete
$\mathbb{Z}_{2}$ symmetry (defined by $\hat{P}_{z}$), arising as a consequence
of the Floquet drive. This $\mathbb{Z}_{2}$ symmetry operation is induced by
flipping the $\hat{\textbf{x}}$-direction of all spins. In fact, it can be
generically demonstrated that this emergent symmetry conspires with the
discrete time-translation symmetry to yield a spatio-temporal ordering of
eigenstates in the Floquet operator itself [24]: at $\gamma\\!=\\!\pi$,
$U_{F}$ can be approximated exponentially-well by

$\displaystyle
U_{F}\simeq\mathcal{U}\hat{P}_{z}\mathrm{e}^{-iT\hat{D}}\mathcal{U}^{\dagger},$
(26)

where $\hat{D}$ is an effective Hamiltonian with the property
$[\hat{P}_{z},\hat{D}]=0$ and $\mathcal{U}$ is a time-independent (though
many-body) rotation close to the identity. Then, the effective drive generated
by $\hat{P}_{z}\mathrm{e}^{-iT\hat{D}}$ obeys spatio-temporal eigenstate
order. While it might be difficult to obtain $\mathcal{U}$ and $\hat{D}$
generically, there is an exception for our long-range interacting model at
$N\mathrm{mod}\leavevmode\nobreak\ 8\\!=\\!0$, where $\mathcal{U}=1$ and
$\hat{D}\equiv\overline{\mathcal{H}}$. For this case, let us denote the joint
eigenstates by $|n,p\rangle$: then
$\overline{\mathcal{H}}|n,p\rangle=E_{n}|n,p\rangle$, and
$\hat{P}_{z}=p|n,p\rangle$ ($p\\!\in\\!\\{\pm 1$}). The Floquet unitary can be
eigen-decomposed as

$\displaystyle U_{F}$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{n,p\in\\{\pm
1\\}}\\!\\!\\!\exp\left[-iT\left(E_{n}+\frac{\pi(1-p)}{2T}\right)\right]|n,p\rangle\langle
n,p|,$ (27)

where each eigenstate $|n,p\rangle$ has a partner $|n,-p\rangle$ whose quasi-
energy is shifted by exactly $\pi/T$. In particular, this implies that any
symmetry-broken initial state will necessarily also break the discrete time-
translation invariance of the Floquet drive, and instead oscillate with a
period of $2T$, exhibiting period-doubling [3].

(ii) Insensitivity to the initial state appears as consequence of spatio-
temporal eigenstate order. In fact, as we will demonstrate below in Sec.
VII.2, any initial state which is not an eigenstate of the Floquet unitary and
possesses a finite $\hat{\textbf{x}}$-magnetization, even mixed states, is
subject to stable long-lived period-doubled oscillations.

(iii) A parametrically long-lived prethermal time-window, as analyzed in
detail in Sec. VIII.

(iv) Robustness to perturbations: when $\gamma=\pi+\varepsilon$ is detuned
from the sweet-spot ($\gamma=\pi$), the system remains symmetry-broken due to
the presence of interactions. We have verified this numerically, using a
Fourier-transform of the time-evolution curves of observables in units of
stroboscopic cycles $M$

$\displaystyle\langle O\rangle(\omega)=\sum_{j=0}^{M-1}\mathrm{e}^{-i\omega
jT}\langle O\rangle(jT),$ (28)

where $\omega=2\pi k/(MT)$ with $k=0,\dots,(M-1)$ and $\langle O\rangle(jT)$
is the expectation value of the observable of interest at stroboscopic times.
The results are displayed in Fig. 3 B and D of the main text for experimental
and theoretical data. Around $\gamma=\pi$ we clearly observe a single rigid
peak at $\omega=\pi/T$ of finite extent ($\approx\pi/5$), corresponding to a
period-doubled oscillation in both panels. In addition, before the system
heats up to infinite temperature, the prethermal time-crystalline order is
also robust to small random perturbations in the time duration $\tau$, as we
demonstrate numerically in Sec. VIII.

It follows therefore that our ${}^{13}\textrm{C}$ nuclear spin system, subject
to the two-frequency drive at $\gamma\approx\pi$, features prethermal discrete
time crystalline (PDTC) order.

###  VII.2  Thermalizing dynamics in the presence of emergent
(quasi-)conservation laws

Consider again the two-cycle time-evolution operator of Eq. (25). We leave the
discussion of unconstrained heating to infinite temperature for Sec. VIII.
Since the effective Hamiltonian $\overline{\mathcal{H}}$ is ergodic, according
to the Eigenstate Thermalization Hypothesis (ETH) [57], we expect an initial
state to thermalize with respect to $\overline{\mathcal{H}}$. After the non-
universal initial transient dynamics, the state of the system is described by
a Gibbs state, which is determined by the corresponding Lagrange multipliers
(the obvious one being the inverse temperature $\beta$). Thermal states are
time-translation invariant, which implies Lagrange multipliers that are
constant in time. Similarly, Floquet prethermal states are described by slowly
and continuously changing Lagrange multipliers [49, 50].

At first sight, this appears at odds with the period-doubling oscillations of
observables, induced by the broken spatio-temporal symmetry near $\gamma=\pi$.
To resolve this conundrum, note first that the thermal properties of the state
are defined w.r.t. the period-doubled unitary, Eq. (25): indeed, on the
timescale $2M$ (discrete) time-translation invariance is inherent. Yet, this
does not a priori imply time-invariance of the state on smaller timescales,
i.e. over a single Floquet cycle $M$, where we observe a strong alternating
pattern in time. As we explain below, the intriguing mechanism that enables a
time-dependence in the (approximate prethermal) Gibbs state can be understood
to arise from a fine interplay between two conspiring emergent symmetries of
$\overline{\mathcal{H}}$, both induced by the two-frequency drive: (i) the
spatio-temporal Ising $\mathbb{Z}_{2}$ symmetry in combination with (ii) the
$\hat{\textbf{x}}$-magnetization quasi-conservation.

Fig. S10:  Numerical Simulation: (a)-(b) Time evolution of different initial
states $|\psi\rangle_{i}^{t_{d}}$ prepared by evolving the fully polarized
state for a time $t_{d}$ with $\mathcal{H}$:
$|\psi\rangle_{i}^{t_{d}}=\mathrm{e}^{-it_{d}\mathcal{H}}\bigotimes_{j=1}^{L}\frac{1}{\sqrt{2}}\left(|\\!\uparrow_{j}\rangle-|\downarrow_{j}\rangle\right)$.
In (b) we display the absolute value of the $\hat{\textbf{x}}$-magnetization
using the data from (a), normalized to its initial magnitude. The simulation
parameters are $L=16$, $N=255$, $\gamma=\pi+0.02$, $J\tau=0.2$, and the
remaining parameters are the same as in Fig. S8. Experiment: (c) Corresponding
experimental time-evolution curves for different initial states obtained using
the same procedure: the initially prepared density matrix (see main text) is
evolved with the bare dipolar Hamiltonian $\mathcal{H}$ up to time $t_{d}$
before the two-frequency drive is started. The inset provides a zoom in of the
marked region to visualize individual curves. The parameters of the experiment
are the same as in Fig. 3 of the main text, and $\gamma=0.96\pi$.

If we initialize the system in a $\mathbb{Z}_{2}$-symmetry-broken state with
non-zero magnetization, we expect the system to prethermalize only w.r.t. the
corresponding symmetry-broken sector [23]. In the grand-canonical formulation
of ETH, the subsystem density matrix $\hat{\rho}_{A}$ in the prethermal
plateau resembles a thermal state with respect to this symmetry-broken sector
of the subsystem (effective) Hamiltonian $\overline{\mathcal{H}}_{A}$ [49,
58]. Notice that, in this case, $\hat{\rho}_{A}$ is not time-translation
invariant under evolution with $U_{F}$: instead,
$U_{F}\hat{\rho}_{A}U_{F}^{\dagger}=\hat{P}_{z}\hat{\rho}_{A}\hat{P}_{z}^{\dagger}=\overline{\hat{\rho}_{A}}$,
where $\overline{\hat{\rho}_{A}}$ corresponds to the (prethermal) subsystem
density matrix in the complementary $\mathbb{Z}_{2}$-symmetry sector. Indeed,
the fact that $\hat{\rho}_{A}$ is thermal only w.r.t. one of the two
$\mathbb{Z}_{2}$ symmetry sectors is equivalent to having non-zero
magnetization expectations of $\hat{\rho}_{A}$ and
$\overline{\hat{\rho}_{A}}$.

More precisely, according to ETH, any given initial state $|\psi_{0}\rangle$
exposed to the drive of Eq. (27) is expected to locally resemble a Gibbs
state. In the thermodynamic limit, a subsystem $A$ will then be described by

$\displaystyle\hat{\rho}_{A}=\frac{1}{\mathcal{Z}}\exp\left(-\beta\overline{\mathcal{H}}_{A}-\mu\mathcal{I}_{x,A}\right),$
(29)

where
$\mathcal{Z}=\mathrm{Tr}[\exp(-\beta\overline{\mathcal{H}}_{A}-\mu\mathcal{I}_{x,A})]$
with the inverse temperature $\beta$. Here $\overline{\mathcal{H}}_{A}$ and
$\mathcal{I}_{x,A}$ are the Hamiltonian and $\hat{\textbf{x}}$-magnetization
restricted to subsystem $A$, respectively. The magnetization potential $\mu$
is determined self-consistently from the initial state:

$\displaystyle\mathrm{Tr}[\overline{\mathcal{H}}\hat{\rho}(t=0)]$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\mathrm{Tr}\big{[}\overline{\mathcal{H}}_{A}\mathrm{e}^{-\beta\overline{\mathcal{H}}_{A}-\mu\mathcal{I}_{x,A}}\big{]},$
(30) $\displaystyle\mathrm{Tr}[\mathcal{I}_{x}\hat{\rho}(t=0)]$
$\displaystyle=$
$\displaystyle\frac{1}{\mathcal{Z}}\mathrm{Tr}\big{[}\mathcal{I}_{x,A}\mathrm{e}^{-\beta\overline{\mathcal{H}}_{A}-\mu\mathcal{I}_{x,A}}\big{]}.$
(31)

So long as $\mu\neq 0$, the state described by Eq. (29) is not invariant under
time evolution with $U_{F}$. Instead, time evolution over a single Floquet
cycle maps $\mu\rightarrow-\mu$ as a result of the spatio-temporal
$\mathbb{Z}_{2}$ symmetry.

Now, consider initializing the system at $t=0$ in a state where,
$\mathrm{Tr}[\mathcal{I}_{x}\hat{\rho}(t=0)]=0$. For such an initial
condition, let us solve Eq. (31) for $\mu$. Because of the
$\hat{\textbf{x}}$-magnetization symmetry, $\overline{\mathcal{H}}_{A}$ and
$\mathcal{I}_{x,A}$ possess a common set of eigenvectors, so that
$\overline{\mathcal{H}}_{A}|s,m\rangle=E_{s,m}|s,m\rangle$, and
$\mathcal{I}_{x,A}|s,m\rangle=m|s,m\rangle$ with $m$ fixing the magnetization
sector. Then, Eq. (31) (assuming vanishing LHS) can be rewritten as

$\displaystyle 0=-\frac{\partial}{\partial\mu}\ln\bigg{(}\sum_{s,m}\exp[-\beta
E_{s,n}-\mu m]\bigg{)}.$ (32)

Since $\overline{\mathcal{H}}_{A}$ has the additional $\mathbb{Z}_{2}$ Ising
symmetry $\hat{P}_{z}$ [cf. Sec. VII], we know that $E_{s,m}=E_{s,-m}$ and
thus Eq. (32) boils down to

$\displaystyle\sum_{s,m>0}2m\mathrm{e}^{-\beta E_{s,m}}\sinh(\mu m)=0,$ (33)

which has a unique solution, namely $\mu\equiv 0$. To sum up, states that obey
$\mathrm{Tr}[\mathcal{I}_{x}\hat{\rho}(t=0)]=0$, are invariant under time-
evolution with $U_{F}$. In particular all eigenstates of $U_{F}$ satisfy this
condition, since they possess a spatio-temporal ordering.

As an example, consider the following initial state:

$\displaystyle|C_{\pm}\rangle=\frac{1}{\sqrt{2}}\bigg{(}\bigotimes_{j=1}^{L}|\\!\rightarrow_{j}\rangle\pm\bigotimes_{j=1}^{L}|\\!\leftarrow_{j}\rangle\bigg{)},$
(34)

where $|\\!\rightarrow_{j}\rangle$ ($|\\!\leftarrow_{j}\rangle$) represents an
eigenstate of $I_{jx}$ with eigenvalue $+1/2$ ($-1/2$). One can convince
onself that $|C_{\pm}\rangle$ obeys the required ordering and is indeed an
eigenstate of $U_{F}$ with $\mathrm{Tr}[\mathcal{I}_{x}\hat{\rho}(t=0)]=0$.
Hence, when starting from the state $|C_{\pm}\rangle$ (or likewise any other
cat-like linear combination of two $\hat{\textbf{x}}$-magnetization
eigenstates), no subharmonic oscillations can appear: the state is time-
translation invariant under evolution with $U_{F}$. In contrast, any symmetry-
broken initial state which obeys $\mu\neq 0$ necessarily breaks the discrete
time-translation symmetry.

To validate this conclusion numerically, we initialize the system in various
different states by time-evolving the $\hat{\textbf{x}}$-polarized state for
some transient time $t_{d}$ with $\mathcal{H}$, i.e.
$|\psi\rangle_{i}^{t_{d}}=\mathrm{e}^{-it_{d}\mathcal{H}}\bigotimes_{j=1}^{L}|\leftarrow_{j}\rangle$.
In that way we obtain experimentally accessible, yet highly non-trivial,
initial states with non-zero $\hat{\textbf{x}}$-magnetization, provided
$t_{d}J\lesssim 1$. The results are displayed in Fig. S10 (a) and (b). We plot
the absolute value of the signal $|\left\langle x\right\rangle|$ for clarity.
All chosen initial states lead to long-lived subharmonic oscillations of the
signal; yet, we find differences in their lifetime: while states with large
initial $\hat{\textbf{x}}$-magnetization tend to increase their lifetime with
increasing $t_{d}$, at small initial $\hat{\textbf{x}}$-magnetization, the
lifetime behavior becomes less systematic. For values very close to zero
initial $\hat{\textbf{x}}$-magnetization we observe no PDTC order (data not
shown).

Experimentally, we pursue the same strategy. We create different initial
states by evolving the rotated starting density matrix $\rho_{0}$ (see main
text) with the bare Hamiltonian $\mathcal{H}$ for different times $t_{d}$ to
obtain a set of different and highly non-trivial initial states, which are
subsequently subject to two-frequency driving. Each state corresponds to a
different run of the experiment. We display the results in Fig. S10 (c). We
show only absolute values, yet all tested initial states create persistent (in
time) and rigid (in $\varepsilon$) subharmonic oscillations with similar
heating times. This is a notable difference to the data obtained using
numerical simulations, where different initial states can induce considerably
different heating times. We hypothesize that the difference in the heating
behavior of various initial states is related to the relatively small system
sizes that can be reached in the numerical simulations.

###  VII.3 Classification of discrete time crystalline order

As outlined in detail in the previous section, time crystalline order is
expected to emerge for all states which initially possess a finite
$\hat{\textbf{x}}$-magnetization expectation value: this is because, within
the prethermal plateau, thermalization happens with respect to
$\overline{\mathcal{H}}_{A}$, which preserves the initial
$\hat{\textbf{x}}$-magnetization. Similar to the discussion of Ref. [40], this
emergent $\hat{\textbf{x}}$-magnetization conservation additionally constrains
thermalization, which implies that no thermal phase transition in
$\overline{\mathcal{H}}_{A}$ can preclude the PDTC order to form. This
mechanism of stabilizing prethermal time-crystalline order is different from
the conventional approach where the formation of PDTC order is solely
dependent on the energy density and the associated inverse temperature $\beta$
of the initial state with no additional protection. In our Floquet system,
such a phenomenon can occur when the repetition number $N$ is decreased and
$\hat{\textbf{x}}$-non-conserving terms $\propto 1/N$ start gaining
importance, cf. Eq. (13).

Importantly, the extent to which $\overline{\mathcal{H}}$ describes the
dynamics of the system is parametrically controlled by the frequency of the
fast $\hat{\textbf{x}}$-drive: increasing the drive frequency (i) ensures a
suppression of $\hat{\textbf{x}}$-magnetization conservation breaking terms
appearing in higher order corrections of the inverse frequency expansion of
the two-cycle effective Hamiltonian, and (ii) it parametrically increases the
lifetime of the prethermal plateau for which the DTC order can be observed
(see Sec. VIII below). However, since our system exhibits long-range
interactions with a critical exponent, we obtain a power-law suppression of
heating rates as a function of frequency (see Sec. VIII). The DTC order in our
system thus combines features of prethermal time-crystalline order [23, 29]
and critical time-crystalline order [11, 9].

##  VIII Heating timescales and duration of the prethermal plateau

The lifetime of a PDTC, like that of any prethermal order, is predominantly
determined by the amount of energy absorbed from the nonequilibrium drive. To
quantify this energy absorption, we define the heating time empirically, as
follows. Given the time evolution curve of some observable $O$, the heating
time is the time required to reach $1/\mathrm{e}$ of the initial value of that
observable. In principle, one would have to define the heating time with
respect to the prethermal value of a given observable. However, since we only
investigate the quasi-conserved $\hat{\textbf{x}}$-magnetization, the initial
and prethermal expectation values coincide.

The concatenated two-frequency drive offers two obvious mechanisms for the
system to heat up to infinite temperature: the slow and fast drives
independently cause energy absorption, each with a potentially different rate
corresponding to its own timescale. Numerically and experimentally, we observe
an interplay of both timescales. However, theoretically we can disentangle
them and analyze each effect separately. To this end, in what follows, we
discuss the influence of heating caused by the slow $\hat{\textbf{z}}$-drive,
and then the fast $\hat{\textbf{x}}$-drive, before we conclude with the
general case for the two-frequency drive.

###  VIII.1 Heating timescales associated with the slow
$\hat{\textbf{z}}$-drive

To analyze heating caused by the slow drive, we assume the (approximate)
effective Hamiltonian emerging from the fast drive, Eq. (21), to be exact.
This allows us to numerically study a periodic drive of the form of Eq. (19).

We further increase the ergodicity of the relatively small systems accessible
in our simulations and diminish the related finite-size effects by adding a
small uniformly distributed noise $\delta$ to the driving period [49, 50]

$\displaystyle
U_{F}=\mathrm{e}^{-iN(\tau+\delta\tau)\overline{\mathcal{H}}}U_{z},$ (35)

where $\delta$ takes a different value in each Floquet cycle.

We then evolve the fully $\hat{\textbf{x}}$-polarized initial state under Eq.
(35) up to $5\times 10^{4}$ kicks. The results are displayed in Fig. S11. We
find a clear power-law scaling of the heating times $\Gamma_{z}^{-1}$ which is
close to a Fermi-golden-rule scaling, $\Gamma^{-1}_{z}\sim\varepsilon^{-2}$.
This can be understood as the $\hat{\textbf{z}}$-drive is effectively
operating in the low-frequency regime for $N\gg 1$ [the short-range
interacting exponential suppression typically occurs in the high-frequency
regime]. Note that, at $\varepsilon\\!=\\!0$, although the energies of many-
body states can differ by multiples of the drive frequency, heating is
suppressed by the complete absence of matrix elements between these states.
Finite $\varepsilon\\!>\\!0$ induces finite matrix elements which in turn
results in intense energy absorption according to Fermi’s Golden rule. Note
that we observe approximately the same golden rule $\varepsilon$-scaling for
different values of $N$ (Fig. S11a-b).

Fig. S11: Numerical simulation – $\varepsilon$ dependence of the heating
timescales $\Gamma_{z}^{-1}$ associated with the slow
$\hat{\textbf{z}}$-drive: (a) $N=300$, (b) $N=400$ around
$\gamma=\pi+\varepsilon$. The system size is $L=14$, and $J\tau=0.2$. Further
simulation parameters are the same as in Fig. S8.  Fig. S12: Numerical
simulation – Time-evolution for random graphs of different system sizes at
$\gamma=-0.039$ for $\delta=0$ (a) and $\delta=0.05$ (b). $N=300$ and
$J\tau=0.2$. Remaining parameters are as in Fig. S8.

The power-law suppressed heating in our system, causes the finite stable
regions around $\gamma=0,\pi$ to shrink much faster, as compared to the
exponential suppression of heating at high-frequencies. However, in practice
the timescales are still parametrically controlled by $\varepsilon$ and the
prethermal physics is governed by the effective Hamiltonian derived in Sec.
VI. In fact, the experimentally observed lifetimes for our PDTC readily exceed
state-of-art lifetimes of PDTCs in the high-frequency regime [28, 29].

We note in passing that finite noise $\delta$ in the driving protocol has no
influence on the scaling of heating times, but it is capable of removing long-
time finite size synchronization effects [49, 50]. This can be seen in Fig.
S11 where finite $\delta$ causes more points to align on the straight line
without changing its slope, whereas points that are already in-line remain
unchanged. Yet another important aspect of finite $\delta$ is shown in Fig.
S12: while the curves of the noise-free evolution in (a) – especially their
long-time heating behaviour – are affected by finite-size effects, the
addition of a small random $\delta$ almost collapses all curves so that the
heating times become basically insensitive to the system size $L$. By
contrast, the experiment offers the advantage of working with a large enough
system size which eliminates finite-size effects for practical purposes, and
hence the extra $\delta$-noise is neither necessary, nor used.

###  VIII.2 Heating timescales associated with the fast
$\hat{\textbf{x}}$-drive

The physics of the thermalization dynamics associated with the fast
$\hat{\textbf{x}}$-drive only, is not immediately obvious due to the interplay
between three effects: (i) the power of the dipolar interaction term is
critical in three dimensions, which implies a logarithmic divergence of the
total energy in the thermodynamic limit. (ii) In the limit $\tau\rightarrow
0$, the effective Hamiltonian of the system approaches an integrable model,
for which the notion of thermalization itself is not well-defined.
Integrability breaking, though, occurs at order $\mathcal{O}(\tau)$ in the
effective Hamiltonian, although its not pronounced at the accessible system
sizes in the numerical simulations. (iii) Even when $J\tau\ll 1$, we still
have $\vartheta\sim\mathcal{O}(1)$ so that the high-frequency regime
($\vartheta\ll 1$ and $J\tau\ll 1$) is practically inaccessible.

For these reasons, it is difficult to make predictions for rigorous bounds on
the heating rates, given the above properties of the model [23]. Nonetheless,
it is possible to study the energy absorption numerically. To this end, we
investigate the drive generated by the repeated application of
$U_{x}U_{\mathcal{H}}$ at $\gamma=0$ with fixed $\vartheta=\pi/2$ for
different values of $J\tau$. From the time-evolution curves we extract the
heating times $\Gamma^{-1}_{\mathrm{min}}(J\tau)$ and display them as a
function of $\tau$ in Fig. S13 for different random graph realizations. We
find that – similar to the slow $\hat{\textbf{z}}$-drive – the fast
$\hat{\textbf{x}}$-drive shows a power-law scaling,
$\Gamma^{-1}_{\mathrm{min}}\propto(J\tau)^{\kappa}$ with a Golden-rule
exponent $\kappa\approx-2$. This result is in agreement with recent
experimental measurements reported in Ref. [37].

Fig. S13: Numerical simulation – Heating times
$\Gamma^{-1}_{\mathrm{min}}(J\tau)$ associated with the fast
$\hat{\textbf{x}}$-drive, extracted from corresponding time-evolution curves
(data not shown), as a function of $J\tau$. Different colors (red, green and
blue) correspond to different random graph realizations. The system size is
$L=14$. The remaining simulation parameters are the same as in Fig. S8.

###  VIII.3 Heating model for the concatenated two-frequency drive

Now that we have analyzed the fast and slow drives independently, we can move
on with the discussion of the heating rates of the more complex concatenated
two-frequency drive.

Fig. S14:  Comparison of scaling behaviour in theory and experiment: (a)-(b)
Numerical simulation – Numerically extracted heating times averaged over 20
graph realizations,
$\overline{\Gamma}^{-1}=N_{\mathrm{seed}}^{-1}\sum_{s=1}^{N_{\mathrm{seed}}}\Gamma_{s}^{-1}$
with $N_{\mathrm{seed}}=20$, as function of $\varepsilon$ at
$\gamma=\pi+\varepsilon$ (a) and $\gamma=0+\varepsilon$ (b). We plot
$\overline{\Gamma}-\overline{\Gamma}_{\mathrm{min}}$ on the $y$ axis to check
for corresponding scaling. $\overline{\Gamma}_{\mathrm{min}}^{-1}$ is defined
as the heating time at $\varepsilon=0$:
$\overline{\Gamma}_{\mathrm{min}}^{-1}\equiv\overline{\Gamma}^{-1}(\varepsilon=0)$.
We fit the data to
$\overline{\Gamma}=g/N\varepsilon^{\lambda}+\overline{\Gamma}_{\mathrm{min}}$.
(fit values of (a): $g=0.56$, $\overline{\Gamma}_{\mathrm{min}}=3.67\times
10^{-5}$, fit values of (b): $g=0.52$,
$\overline{\Gamma}_{\mathrm{min}}=3.40\times 10^{-5}$). The simulation
parameters are $L=16$, $N=255$, $J\tau=0.2$. Remaining parameters are as in
Fig. S8. (c)-(d) Experiment – Corresponding heating rates extracted from
experimental data for $\gamma=\pi-\varepsilon$ (c) and $\gamma=\varepsilon$
(d). In (d) the three smallest data points are not included in the fit since
these curves did not reach the $1/\mathrm{e}$-threshold. We fit the data to
$\Gamma=g/N\varepsilon^{\lambda}+\Gamma_{\mathrm{min}}$ (fit parameters of
(c): $g=6.15$, $\Gamma_{\mathrm{min}}=1.75\times 10^{-5}$, fit parameters of
(d): $g=4.17$, $\Gamma_{\mathrm{min}}=7.34\times 10^{-6}$). The remaining
parameters are the same as in Fig. 3 of the main text.

To understand the combined heating timescales we deploy a simple model:
assuming that there are no correlation effects between the different drives
(which applies in the time-scale separated regime $N\gg 1$), a simple model
that approximates the dynamics of the magnetization decay can be formulated
using a discretized-in-time rate equation that captures the influence of the
$MN$-th kick on the magnetization:

$\langle x\rangle((MN\\!+\\!1)\tau)-\langle
x\rangle(MN\tau)=-\bigg{[}\frac{\Gamma_{z}(\varepsilon)}{N}+\Gamma_{\mathrm{min}}(\\!J\tau\\!)\bigg{]}\langle
x\rangle(MN\tau).$ (36)

Here we denote the heating rates of the slow $\hat{\textbf{z}}$-drive and the
fast $\hat{\textbf{x}}$-drive by $\Gamma_{z}(\varepsilon)$ and
$\Gamma_{\mathrm{min}}(J\tau)$, respectively. Recalling that $N$ and $M$ are
the repetition numbers of the $\hat{\textbf{x}}$\- and $\hat{\textbf{z}}$\-
drives, Eq. (36) has the simple solution

$\langle
x\rangle(MN\tau)=\exp\left[-\left(\frac{\Gamma_{z}(\varepsilon)}{N}+\Gamma_{\mathrm{min}}(J\tau)\right)MN\tau\right].$

Notice that the exact functional form of the true time-dependence of the
magnetization may not resemble a simple exponential. However, if we are only
interested in heating timescales and taking into account the experimentally
observed mono-exponential behavior, Eq. (36) provides a valid approximation,
since the time required to reach a value of $1/\mathrm{e}$ is exactly given by
$\Gamma^{-1}=[\Gamma_{z}(\varepsilon)/N+\Gamma_{\mathrm{min}}(J\tau)]^{-1}$.
At $\varepsilon\\!=\\!0$ the slow drive is ineffective so that
$\Gamma^{-1}|_{\varepsilon=0}=\Gamma_{\mathrm{min}}^{-1}(J\tau)$, and we
recover the heating rate of the fast drive. For a fixed fast-drive frequency,
we can thus identify the minimum heating rate as
$\Gamma_{\mathrm{min}}(J\tau)$. Together with the conclusions from Sec. VIII.1
we, therefore, expect a scaling approximately given by

$\displaystyle\Gamma^{-1}=\frac{1}{\frac{g}{N}\varepsilon^{2}+\Gamma_{\mathrm{min}}(J\tau)},$
(37)

with some numerical constant $g$. Rearranging this expression as
${\Gamma}=g/N\varepsilon^{\lambda}+{\Gamma}_{\mathrm{min}}$ we fit a straight
line on a log-log plot to both numerical and experimental data, with fitting
parameters $g$, $\lambda$, and $\Gamma_{\mathrm{min}}$.

In Fig. S14 we present fits of the heating times extracted from both numerical
simulations (Fig. S14 (a)-(b)) and experimental data (Fig. S14 (c)-(d)). The