## Quantum thermalization and Floquet engineering in a spin ensemble with a clock transition

Mi Lei, Rikuto Fukumori, Chun-Ju Wu, Edwin Barnes, 3 4, 1, 2,

1, 2 1, 2 1, 2 3

Sophia Economou, Joonhee Choi, ∗ and Andrei Faraon † 1 Kavli Nanoscience Institute and Thomas J. Watson, Sr., Laboratory of Applied Physics, California Institute of Technology, Pasadena, California 91125, USA 2 Institute for Quantum Information and Matter, California Institute of Technology, Pasadena, California 91125, USA 3 Department of Physics, Virginia Technology, Blacksburg, VA 24061, USA 4

Department of Electrical Engineering, Stanford University, Stanford, California 94305, USA

Studying and controlling quantum many-body interactions is fundamentally important for quantum science and related emerging technologies. Optically addressable solid-state spins offer a promising platform for exploring various quantum many-body phenomena due to their scalability to a large Hilbert space. However, it is often challenging to probe many-body dynamics in solid-state spin systems due to large on-site disorder and undesired coupling to the environment. Here, we investigate an optically addressable solid-state spin system comprising a strongly interacting ensemble of millions of ytterbium-171 ions in a crystal. Notably, this platform features a clock transition that gives rise to pure long-range spin-exchange interactions, termed the dipolar XY model. Leveraging this unique feature, we investigate quantum thermalization by varying the relative ratio of interaction strength to disorder, dynamically engineering the XY model into other many-body Hamiltonian models, and realizing a time-crystalline phase of matter through periodic driving. Our findings indicate that an ensemble of rare-earth ions serves as a versatile testbed for many-body physics and offers valuable insights for advancing quantum technologies.

Many-body quantum phenomena arise when multiple particles interact, playing a crucial role in contemporary physics and engineering 1,2 . Among the systems exhibiting many-body physics, optically addressable solid-state spins have become a particularly rich platform to study due to their scalability to very large numbers of spins 3 and their potential applications in quantum simulations 4 , sensing 5,6 , and information processing . 7

To achieve versatility in quantum applications, it is essential to manipulate the time evolution of a many-body system in a programmable manner . However, achieving 8 local individual control of spins in solid-state systems is often challenging due to their nanometric proximity. To address this issue, a global control sequence-composed of either pulsed or continuous driving-can be applied periodically to the entire system, effectively manipulating many-body dynamics 9,10 . This periodic driving method, known as Floquet Hamiltonian engineering, enables the efficient realization of various many-body Hamiltonians with different interaction types and strengths, transformed from the original system Hamiltonian 11-14 . Notably, periodically driven quantum systems can also exhibit a variety of exotic non-equilibrium phases of matter, such as Discrete Time Crystals (DTCs) 15-19 , stabilized by many-body interactions.

as promising candidates for a many-body testbed due to their highly coherent optical and spin transitions at cryogenic temperatures, scalability, and ease of integration into photonic devices 20-22 . Moreover, the diversity in REI species, hosts, and concentrations, combined with naturally arising dipole-dipole interactions 23 or engineered cavity-mediated interactions 21 , introduces novel experimental control parameters and offers theoretical insights for exploring large-scale quantum many-body systems and their dynamics. However, despite numerous spectroscopic studies on large ensembles of REIs that have focused on identifying materials with high coherence times 24,25 for quantum technologies such as quantum transducers 26 and memories 27 , less effort has been devoted to the microscopic understanding, control, and engineering of these systems.

In this context, a strongly interacting ensemble of spins with robust global control protocols can be employed to explore many-body physics and out-of-equilibrium dynamics. Rare-earth ions (REIs) doped in solids emerge

∗ joonhee.choi@stanford.edu

† faraon@caltech.edu

In this work, we report on the characterization and control of quantum many-body dynamics in a dense ensemble of approximately 10 6 ytterbium-171 ions in a nanophotonic cavity with a yttrium orthovanadate host crystal ( 171 Yb 3+ :YVO ). 4 These high-density REIs are randomly positioned with an average distance of ≈ 9 nm (or equivalently, ≈ 86 ppm in concentration) in an effective three-dimensional volume defined by the cavity mode (Fig. 1a). Each Yb ion provides both microwave and optical transitions for spin state manipulation and stateselective readout, respectively. The nanophotonic optical cavity enables high-fidelity initialization and fast readout of Yb ions, while coherent microwave control is achieved via a coplanar waveguide. Specifically, the ground and optically excited states of Yb comprise {| 0 〉 , | 1 〉 , | Aux 〉} and {| 0 〉 e , | 1 〉 e , | Aux 〉 e } , respectively, resulting from the hybridization of its electron and nuclear spins (Fig. 1b) 35 .

FIG. 1. Many-body platform based on rare-earth ions.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-98e1ca21-6ebc-461b-886b-be0290b20bb9.jpg)

a, An ensemble of ≈ 10 6 rare-earth ion spins (green arrows) is coupled to an optical cavity (red mirrors) and a coplanar waveguide (gold stripes). The optical input and output (red arrows) are sent from the same port for spin initialization and readout. b, The energy levels of each spin, with four ground states ( 0 , | 〉 | 1 , 〉 and | Aux ) and four optically excited states 〉 ( 0 | 〉 e , | 1 〉 e , and | Aux 〉 e ). Both | Aux 〉 and | Aux 〉 e are doubly degenerate states, whereas | 0 , 〉 | 0 〉 e , | 1 , 〉 and | 1 〉 e are clock states. We define | 0 〉 and | 1 〉 as a qubit and read out their state via the optical transition around 984 nm. c, An interacting spin ensemble where spins i and j interact via pairwise interactions with strength J ij . d, Control of the average interaction strength J by varying the population distribution between the qubit manifold {| 0 〉 , | 1 〉} , and the auxiliary states | Aux 〉 via optical pumping (Supplementary Information). e, Our spin system permits pure spin-exchange interactions with no Ising interactions. f, Benchmarking our many-body platform to other solid-state electronic spin systems. The bar chart shows the ratio of the average interaction strength to disorder, J/W , for different systems: 1: REIs with nonzero first-order-Zeeman shift (Non-ZEFOZ REIs) 28,29 ; 2: Defects in hexagonal boron nitride (hBN) 30 ; 3: Pink diamond 31,32 ; 4: Black diamond 33 ; 5: P1 centers 34 ; 6: ZEFOZ REIs (our work). Inset: comparison of the absolute values of J and W .

In this study, the spin transitions are defined within the ground state manifold using {| 0 〉 , | 1 〉} , serving as an effective 'qubit' with a microwave transition frequency of 675 MHz. For spin state readout, we utilize the cavityenhanced, resonant optical transition between | 1 〉 and | 0 〉 e at a wavelength of 984 nm, enabling the optical detection of state-selective photoluminescence signals (Supplementary Information).

In essence, the goal of studying many-body physics in an experiment is to observe coherent quantum phenomena over extended durations governed by the target unitary dynamics of the system, while minimizing undesired incoherent coupling to external environments. In this regard, our REI platform offers distinct advantages over other solid-state spin systems for exploring the dynamics of a closed many-body quantum system, as outlined below.

First, the chosen spin states, {| 0 〉 , | 1 〉} , are first-order insensitive to external electromagnetic fluctuations at zero magnetic fields, known as the 'clock' transition 35 . Consequently, the spins are less susceptible to both decoherence and inhomogeneity induced by the external environment. Second, the electronic spin g factor, which determines the strength of the dipole moment, is approximately three times higher than that of a single electron 36 . This results in stronger dipole-dipole interactions between REIs (Fig. 1c). Third, we can control the overall interaction strength of the system, J , defined using the average nearest-neighbor distance (Supplementary Information), by adjusting the effective density of REIs within the qubit manifold {| 0 〉 , | 1 〉} (Fig. 1d). This manipulation is achieved by varying the population distribution between the qubit manifold {| 0 〉 , | 1 〉} and the nonparticipating auxiliary ground states | Aux 〉 through optical initialization (Supplementary Information). Lastly, unlike conventional dipole-dipole interactions, which include both spin-exchange and energy-shifting Ising interactions, our spin system realizes a pure spin-exchange interaction without the Ising component (Fig. 1e). This Hamiltonian is known as the dipolar XY model 37 , which has recently been experimentally investigated in various platforms to explore fundamental many-body phenomena and their applications, such as continuous symmetry breaking 38 and spin squeezing 39 . The naturally occurring XY model in three-dimensional dipolar systems offers a novel and complementary configuration for studying similar physics in the solid-state setting.

The combination of strong spin-spin interactions, weak coupling to the external bath, cavity-enhanced fast readout, and effective spin density control positions our system as a promising platform for studying many-body spin dynamics. Notably, our spin system features the highest ratio of interaction strength to transition inhomogeneity among solid-state electronic spin systems, underscoring its significance in probing many-body dynamics where collective interaction effects are minimally hindered by on-site disorder (Fig. 1f).

FIG. 2. Characterization of decoherence dynamics. We characterize spin system decoherence arising from interaction and disorder using spin-echo, Ramsey, and /epsilon1 -CPMG sequences. We compare three different cases: an interacting spin ensemble with large J (blue) and small J (orange), and an isolated single ion as a reference (black). a, In the spin-echo measurement, spin-spin interactions drive decoherence with decay rates dependent of J , while disorder is decoupled to the zeroth order. b, In the Ramsey measurement, both interaction and disorder contribute to decoherence. In a and b , coherence is measured as a function of the free evolution time, τ , and normalized to their respective maximum coherence. The last π/ 2 pulse with a variable rotation axis θ is used to extract coherence from photoluminescence signals (Supplementary Information). c , In the /epsilon1 -CPMG measurement, a π + /epsilon1 pulse is applied k times at a fixed τ . Coherence is measured as a function of the rotation angle offset, /epsilon1 , and then normalized to the maximum coherence observed in the spin-echo measurement. Note that /epsilon1 can control the decoherence contributions caused by disorder and interaction. The interacting spin ensemble reveals unconventional behavior where coherence is maximized when /epsilon1 ≈ ± π/ 2. The error bars are obtained from fits to the experimental data, and simulations are shown as solid lines (Supplementary Information). τ = 300 ns and k = 8 are chosen for both the large and small J cases, while τ = 5 8 . µ s and k = 400 are chosen for the single ion. In c , we rescale the simulation results by an /epsilon1 -independent prefactor &lt; 1 to match the experimental data.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-053469ca-fb0d-496c-9bcb-3c4acb3737ae.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-c1ae2940-5d6a-4168-9416-cd985b7c029b.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-70138695-2e0f-44cc-a535-e9fcbe2b4e3a.jpg)

## CHARACTERIZATION OF A STRONGLY INTERACTING SPIN SYSTEM

The very first step before we utilize our REI spin system as a many-body testbed is to quantitatively characterize the interaction and disorder strengths in the system. To this end, we describe the dynamics of our interacting spin system using the following Hamiltonian (in a rotating frame), ˆ , H defined within the qubit manifold {| 0 〉 , | 1 〉} :

ˆ S i ± = ˆ S i x ± iS ˆ i y are the creation and annihilation operators for spin excitation of ion i . We define the strength of on-site disorder, W , as the full width at half maximum of the probability distribution function of ∆ . i

$$\hat { H } = \hat { H } _ { \text{dis} } + \hat { H } _ { \text{int} } \quad \quad ( 1 ) \quad \underset { / \mathfrak { w } } { o u }$$

where ˆ H dis = ∑ N i ∆ ˆ i S i z is the on-site disorder Hamiltonian with spin detuning ∆ i for ion i , and ˆ H int = ∑ N ij, i&gt;j J ij ( ˆ S S i x ˆ j x + ˆ S S i y ˆ ) = j y ∑ N ij, i&gt;j J ij 2 ( ˆ S i + ˆ S j -+ ˆ S i -ˆ S j + ) is the long-range, dipolar spin-exchange Hamiltonian with position- and orientation-dependent pairwise interaction strength J ij between two ions i and j (see Supplementary Information for details). Here, ˆ S i µ is the spin1/2 operator of ion i along the µ -axis ( µ = x, y, z ), and

In strongly interacting spin systems, the decoherence of a spin ensemble is influenced not only by the random, inhomogeneous on-site fields but also by interactions with the rest of the spin system, which act as an intrinsic bath 40 . To isolate the decoherence effect arising from spin-spin interactions, we employ the celebrated spin-echo sequence 41 , which effectively cancels out the on-site disorder in the Hamiltonian at short times (Fig. 2a). Here, we observe rapid decay of ensemble coherence within ≈ 1 µs , which is much shorter than that of an isolated single spin independently measured from a reference sample. We confirm that the decay time is highly dependent on the effective spin density, where the decoherence rate increases as J increases (orange/blue markers, Fig. 2a; see Supplementary Information for additional experimental data). This implies that intrinsic spin-spin interactions dominate the decoherence mecha-

nisms.

To extract the effective spin density and the corresponding average interaction strength, J , we conduct numerical simulations based on a closed many-body system (Supplementary Information). We find that the simulations show good agreement with the experiment, revealing spin concentrations of ≈ 46 ppm and ≈ 25 ppm for the cases with large J ≈ 2 π × 0 35 MHz and small . J ≈ 2 π × 0 19 . MHz, respectively (orange/blue lines, Fig. 2a). These spin densities within the qubit manifold are reasonable given the total spin densities of all ground states of ≈ 86 ppm, as independently measured by mass spectrometry 21 .

Having characterized the interaction strengths, we now focus on identifying the on-site disorder strength, W , using the disorder-sensitive Ramsey sequence (Fig. 2b). Experimental data show that the Ramsey signal decays faster than the spin-echo signal in both the small and large J regimes, due to the additional contribution of disorder-induced decoherence. By comparing the experimental data to the corresponding Ramsey sequence simulation, we estimate an on-site disorder strength of W ≈ 2 π × 0 65 MHz, which is independent of interaction . strengths (orange/blue lines, Fig. 2b).

To further corroborate the significance of spin-spin interactions in our system, we employ the so-called /epsilon1 -CPMG sequence 42 , a variant of the conventional CPMG sequence (Fig. 2c). Specifically, the /epsilon1 -CPMG sequence purposely uses an imperfect echo pulse with a rotation angle of π + /epsilon1 with nonzero /epsilon1 . For either an isolated spin or a non-interacting spin ensemble affected only by disorder, such imperfect spin rotation with a sizable /epsilon1 results in non-ideal dynamical decoupling, leading to rapid decoherence compared to /epsilon1 = 0 (black markers, Fig. 2c). In contrast, we observe drastically different behavior in our high-density spin system, where nonzero /epsilon1 pulses better preserve ensemble coherence (orange/blue markers, Fig. 2c).

This feature originates from the different sensitivities of the /epsilon1 -CPMG sequence to on-site disorder and manybody interactions (Supplementary Information). Specifically, when /epsilon1 = 0, spin-spin interactions still lead to rapid decoherence of the spin ensemble because the π pulses do not alter their spin-exchange interaction Hamiltonian, ˆ H int (whereas the effects from the on-site disorder Hamiltonian, ˆ H dis , can be decoupled). However, when /epsilon1 = ± π/ 2, effectively corresponding to a π/ 2 pulse, ˆ H int transforms into a Hamiltonian that is the sum of the Ising interaction along the y -axis and the Heisenberg interaction (Supplementary Information). The initially prepared spins along the y -axis then become an eigenstate under this Hamiltonian. Consequently, using π/ 2 pulses instead of π pulses in the /epsilon1 -CPMG sequence significantly extends the coherence of the interacting spin ensemble (Supplementary Information). Using the Hamiltonian parameters extracted from the Ramsey and echo sequences, we reproduce the coherence dependence on /epsilon1 that matches experimental data (orange/blue

## MICROSCOPIC UNDERSTANDING OF DECOHERENCE MECHANISMS

When probing the spin-echo dynamics over longer timescales, we observe that the decoherence profile starts deviating from a simple exponential decay, displaying much slower relaxation at later times (orange markers, Fig. 3b). Remarkably, a numerical simulation with calibrated disorder and interaction strengths (Model III) shows excellent agreement with the experimental data even at longer timescales (orange line, Fig. 3b). To comprehend these late-time observations, we further consider two simpler theoretical models, I and II, which consider two-body and many-body systems with no disorder, respectively (Fig. 3a; see Supplementary Information for details). We find that both fail to capture the late-time slowdown in decoherence (purple/red lines, Fig. 3b).

Notably, the cross-over from fast to slow decoherence occurs after a time approximately given by the inverse

FIG. 3. Microscopic understanding of decoherence mechanisms. a, Schematic of Models I, II, and III for spinecho sequence simulations. Model I considers only the pair of spins with the largest interaction strength and excludes disorder. Models II and III involve an ensemble of N spins with many-body interactions; Model II excludes disorder, while Model III includes disorder with strength W = 2 π × 0 65 MHz . (calibrated from Fig. 2b). b, Comparison of the experimental spin-echo data (markers) against the three models. The case with interaction strength J ≈ 2 π × 0 19 MHz is considered for . comparison. In each realization of the Monte Carlo numerical simulations, N = 9 spins are randomly positioned based on the lattice structure and the given spin density. The interaction strengths between each pair of spins are calculated based on their positions (Supplementary Information). The error bars are obtained from fits to the experimental data, and simulations are shown as solid lines.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-b4fbb7e2-9c78-46f8-8ce3-2b51e0c807f9.jpg)

FIG. 4. Controlling system evolution via dynamic Hamiltonian engineering. a, Hamiltonian engineering protocol based on control pulse sequences. The first π/ 2 pulse around the x -axis initializes the system into a globally polarized state along the y -axis. The spin states evolve as | Ψ( T ) 〉 = ˆ U | Ψ(0) 〉 over the total interrogation time T , where ˆ U = e -iH ˆ eff T is the timeevolution operator governed by the sequence-dependent effective Hamiltonian, ˆ H eff , obtained through averaged Hamiltonian theory (Supplementary Information). The spin polarization along the y -axis is extracted by fitting the photoluminescence signal as a function of the rotation axis θ of the last π/ 2 pulse (Supplementary Information). b, Ramsey, spin-echo, WAHUHA-echo, and spin-locking sequences for dynamic Hamiltonian engineering. A Rabi frequency of Ω y ≈ 2 π × 10 MHz is used in the spin-locking sequence. The base sequence period of the WAHUHA-echo sequence is 132 ns, limited by finite pulse durations (Supplementary Information). c, We decompose the effective Hamiltonian of each sequence into a sum of different Hamiltonians (see Eq. (2)). The relative weights between the different Hamiltonian terms are represented by the lengths of the bars in the plot. See the main text for details. d, Comparison of decoherence profiles as a function of T under different control sequences. Coherence is normalized to the maximum coherence for each sequence, with error bars indicating the standard deviation of the experimental data. Solid lines represent simulations for the spin-echo and Ramsey sequences, while single-exponential fits are applied to the spin-locking and WAHUHA-echo sequences.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-944f4673-0dde-49f7-8bfb-10a82c4e4811.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-69a0bc43-d2ae-472f-9381-f4a4495604e9.jpg)

disorder strength, τ ∼ 1 /W (Fig. 3b). We attribute this to the breakdown of perfect decoupling of on-site disorder when the free evolution time becomes comparable to the inverse disorder strength, implying the critical role of on-site disorder in late-time many-body dynamics. Specifically, high-order disorder effects, arising from the non-zero commutator between the disorder and interaction Hamiltonians ([ ˆ H dis , H ˆ int ] = 0), become pronounced at late times, imposing energetic penalties on spin-exchange processes and thereby slowing down the decoherence rate (Supplementary Information).

## CONTROL OF SYSTEM EVOLUTION VIA DYNAMIC HAMILTONIAN ENGINEERING

/negationslash

Our microscopic analysis allows us to dissect the underlying decoherence mechanisms of the interacting spin system influenced by both interactions and disorder. Moreover, the analysis strongly supports an exceptionally high degree of coherence in a many-body regime, as evidenced by its excellent agreement with simulations based on a closed many-body system.

Having established that the decoherence dynamics of an individual spin are governed by interactions within the system, we use the decoherence profile as a proxy to explore dynamically engineered many-body Hamiltonians with control pulse sequences (Fig. 4a-c). Intuitively, each spin can be viewed as a quantum sensor sensitive to interactions with surrounding spins, experiencing interactioninduced dephasing, which is often termed quantum thermalization 43 .

Specifically, we prepare an initial state in a globally polarized state along the y -axis, let the spin system evolve under a control pulse sequence over an interrogation time T , and apply a π/ 2 pulse with a variable rotation axis θ at the end of the sequence to measure the mean coherence of all individual spins. Through this experiment, we can probe the decoherence dynamics of a spin ensemble as a function of T , defined as the average spin polarization along the y -axis, P T ( ) = 1 N 〈 Ψ( T ) | ˆ S tot y | Ψ( T ) . 〉 Here, ˆ S tot y = ∑ N i ˆ S i y is the total spin operator along

the y -axis with N being the total number of spins, and | Ψ( T ) 〉 = e -iH ˆ eff T | Ψ(0) 〉 is the quantum state at time T obtained through time evolution under the dynamicallyengineered effective Hamiltonian ˆ H eff .

Crucially, with dynamic Hamiltonian engineering of the original dipolar XY Hamiltonian, we can realize a wide class of different effective Hamiltonians, ˆ H eff , parameterized as

$$\underset { \mu = x, y, z } { \text{arnererze as} } & \text{androd} \\ \hat { H } _ { \text{eff} } = w _ { \text{Heis} } \hat { H } _ { \text{Heis} } + \sum _ { \mu = x, y, z } \left ( w _ { \text{on-site} } ^ { \mu } \hat { H } _ { \text{on-site} } ^ { \mu } + w _ { \text{Ising} } ^ { \mu } \hat { H } _ { \text{Ising} } ^ { \mu } \right ) \cdot \text{on many} \\ \text{interact} \text{i} \\ \text{broken s} ] \\ \text{ian con} ;$$

Here, w Heis , w on-site , and w Ising are relative weights between the different Hamiltonians: the Heisenberg Hamiltonian ˆ H Heis = ∑ N ij, i&gt;j J ij /vector S i · /vector S j , the on-site Hamiltonian ˆ H µ on-site = ∑ N i h S i µ ˆ i µ , and the Ising Hamiltonian ˆ H µ Ising = ∑ N ij, i&gt;j J ij ˆ S S i µ ˆ j µ , where /vector S i = ( ˆ S , S i x ˆ i y , S ˆ ) is i z a vectorized spin operator, and /vector h i = ( h , h i x i y , h i z ) is an effective on-site field for ion i . Note that the Heisenberg Hamiltonian ˆ H Heis is isotropic, whereas the other two Hamiltonians exhibit directionality along the µ -axis.

As depicted in Figs. 4b and 4c, we implement and compare four different effective Hamiltonians using Ramsey, spin-echo, Waugh-Huber-Haberlen (WAHUHA)-echo 44 , and spin-locking sequences. The on-site Hamiltonian is determined by the inhomogeneous disordered field along the z -axis ( /vector h i = (0 0 ∆ )) for the Ramsey and spin-echo , , i sequences, and the homogeneous control field along the y -axis with strength Ω y ( /vector h i = (0 Ω , y , 0)) for the spinlocking sequence. Meanwhile, the WAHUHA-echo sequence is designed to realize the pure Heisenberg Hamiltonian (to leading order), thereby protecting initially polarized spins from interaction-induced decoherence. This protection arises because a globally polarized spin state is an eigenstate of the Heisenberg Hamiltonian 10 . Additionally, the spin-locking sequence prevents spin dephasing by employing a strong pinning field along the y -axis with strength Ω y /greatermuch J after the initial state preparation, effectively arresting spin dynamics.

In line with these theoretical expectations, both the WAHUHA-echo and spin-locking sequences demonstrate significantly prolonged coherence times compared to the Ramsey and spin-echo sequences, confirming the effective engineering of the underlying many-body Hamiltonian (Fig. 4d). However, we observe that the coherence of the WAHUHA-echo and spin-locking sequences still decays over time, with 1 /e time constants of ≈ 20 µs and ≈ 73 µs , respectively. We attribute these decays to imperfections in spin ensemble control caused by finite pulse duration and rotation angle errors, as well as fast time-dependent fluctuations of on-site fields originating from the external spin bath (Supplementary Information). These issues could be addressed by implementing a more robust control sequence and improving control pulse fidelity 10 .

.

## SIGNATURES OF DISCRETE TIME-CRYSTALLINE PHASE

Studying non-equilibrium dynamics under periodic driving unveils the microscopic mechanisms of driven many-body phenomena, such as DTCs 15-19 , revealing an interplay among disorder, dimensionality, interactions, and robustness to both intrinsic and extrinsic perturbations. The canonical model for observing DTCs relies on many-body localization (MBL) and disordered Ising interactions, which lead to the emergence of symmetrybroken spin-glass ordering 45 . Since our system Hamiltonian consists of pure spin-exchange interactions, we employ a spin-locking sequence to effectively realize an Ising interaction along the y -axis (Fig. 5a). In particular, the spin-locking sequence enables long-lived many-body coherence, providing an ideal setting to study the emergent dynamics at late times under periodic driving (Supplementary Information). While we do not expect our threedimensional spin system to exhibit MBL, we investigate the possibility of inducing time-crystalline behavior in a disordered system, akin to the critical DTCs where the relaxation of spin ordering occurs exponentially slowly 46 .

Specifically, we apply a Floquet pulse sequence for probing DTCs as follows: we first prepare the globally polarized spin state along the y -axis, immediately drive the polarized spins along the same direction for a variable spin-locking duration τ , rotate every spin by an angle π + /epsilon1 around the x -axis, and repeat this 'spin-lockand-rotate' operation k times (Fig. 5a). Subsequently, using the last π/ 2 pulse with two opposite phases of x and -x , we measure the normalized spin polarization, P k ( ), at stroboscopic times t = kτ , where k is an integer. The repeated base sequence imposes a discrete timetranslation symmetry with a period of τ , during which spins are allowed to interact (ignoring the finite spinrotation duration), providing the desired setting to test whether our many-body system can break this discrete time-translation symmetry.

/negationslash

Experimentally, we find that when /epsilon1 ≈ 0, corresponding to near-perfect π rotations, with a short interaction period τ ≈ 0, the spin polarization, P k ( ), oscillates up and down over time, exhibiting a 2 τ -periodic oscillation (Fig. 5b). The decay envelope of the oscillation is caused by imperfections in the π pulse, as well as the dephasing of the spin-locking signal itself. The discrete Fourier transform (DFT) of the oscillating P k ( ) reveals a subharmonic peak at a frequency of ν = 1 2 in units of the base period τ . However, this subharmonic peak arises trivially from the finetuned condition of /epsilon1 = 0; as soon as we introduce a non-zero systematic rotation angle error, /epsilon1 = 0, as a perturbation to the system, the subharmonic response is disrupted, and instead, the system exhibits a beat note where the corresponding DFT spectrum shows /epsilon1 -dependent frequency splitting (Fig. 5c). Crucially, however, the subharmonic oscillation of the polarization can be restored by allowing for a longer interaction time of τ ≈ 0 145 . × 2 π J ≈ 425 ns at a nonzero /epsilon1 , which is in-

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-6e202aa0-ad31-4b64-938f-2aac061b3115.jpg)

(

)

FIG. 5. Experimental signatures of discrete timecrystalline states. a, Floquet pulse sequence for probing DTC states. All spins are initially polarized along the y -axis by the first π/ 2 pulse around the x -axis, followed by the periodic repetition of a base sequence with periodicity τ , applied k times. The base sequence contains spin-locking along the y -axis for a duration τ and a global spin rotation by an angle of π + /epsilon1 around the x -axis (ignoring the finite spin-rotation duration). The photoluminescence signals, C + and C -, are collected under the last π/ 2 pulse with opposite phases, x and -x , respectively. We estimate the normalized spin polarization, P k ( ) = | C + -C -| | C + + C -| , at stroboscopic times t = kτ . b-d, The polarization, P k ( ), as a function of Floquet cycle k , and the corresponding Fourier spectrum, | S ν ( ) | 2 , are shown for various values of τ and /epsilon1 : b , τ ≈ 0, /epsilon1 ≈ 0; c , τ ≈ 0, /epsilon1 ≈ 0 03 . π ; d , τ ≈ 425 ns, /epsilon1 ≈ 0 03 . π . The subharmonic oscillation is observed in d , despite the nonzero angle offset /epsilon1 , attributed to the stabilization of DTC states by spin-spin interactions. e, DTC phase diagram constructed using the subharmonic peak intensity at ν = 0 5, i.e., . | S ν ( = 0 5) . | 2 . We observe characteristic linear phase boundaries (dashed lines), determined by identifying the critical perturbation strength (markers) where | S ν ( = 0 5) . | 2 &lt; 0 4 (Supplementary Information). .

dicative of a DTC phase exhibiting robustness against perturbations (Fig. 5d).

We proceed to investigate the stability of this subharmonic behavior more systematically for various values of interaction time τ and perturbation strength /epsilon1 by constructing a phase diagram using the subharmonic peak intensity at ν = 0 5 (Fig. 5e). . The resulting DTC phase diagram shows a characteristic linear phase boundary (red dashed lines, Fig. 5e), consistent with observations from other many-body platforms 18,19 . We further substantiate the robustness of the observed DTC phases by confirming the persistent subharmonic oscillations when varying the initial spin states through global rotations away from the y -axis (Supplementary Information).

## CONCLUSION AND OUTLOOK

Our experimental demonstrations showcase that a REI system provides a versatile and flexible testbed for manybody physics. We envision that REI platforms hold greater potential to serve as large-scale analog quantum simulators in the solid state, offering unique features compared to other solid-state spin systems. First, REI platforms offer flexible engineering options, including a variety of fabrication-friendly choices for ion species, host crystals, and doping concentrations spanning a wide dynamic range from ppb to a few percent. Second, codoping different REI species allows for the simultaneous engineering of two distinct groups of many-body systems, facilitating the study of heterogeneous spin-spin interactions 47-49 , analogous to the dual-species experiments using neutral atom arrays 50 . Third, proximal nuclear spins within the crystal can serve as an additional quantum register. The nuclear spin-spin interactions can be mediated through engineered hyperfine interactions with a Yb ion, offering additional functionality such as a quantum memory 36 . Lastly, all of these control knobs, combined with on-chip integration using nanotechnology 21,22 , enable a scalable REI system, opening up a range of practical applications from quantum simulation to networking and sensing.

## ACKNOWLEDGMENTS

We thank S. L. N. Hermans and T. Xie for reading the paper and providing useful feedback, and A. Ruskuc and A. Beckert for useful discussions. Funding: This work was supported by U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Co-design Center for Quantum Advantage (contract number DE-SC0012704), Institute of Quantum Information and Matter, an NSF Physics Frontiers Center (PHY-1733907) with support from the Moore foundation EPI program, NSF UCLA, NSF QUIC-TAQS (Award Number: 2137984). The initial device and experimental setup development was

supported by the Office of Naval Research award no. N00014-19-1-2182 and N00014-22-1-2422. The device nanofabrication was performed in the Kavli Nanoscience Institute at the California Institute of Technology. R.F. acknowledges the support from the Quad fellowship. E.B. acknowledges support from the National Science Foundation (grant no. 1847078). J.C. acknowledges support from the Terman Faculty Fellowship at Stanford. Author Contributions: M.L., J.C., and A.F. conceived the idea and experiment. M.L. built the ex-

- [1] J. Eisert, M. Friesdorf, and C. Gogolin, Quantum manybody systems out of equilibrium, Nature Physics 11 , 124 (2015).
- [2] H. Tasaki, Physics and mathematics of quantum manybody systems , Vol. 66 (Springer, 2020).
- [3] D. D. Awschalom, R. Hanson, J. Wrachtrup, and B. B. Zhou, Quantum technologies with optically interfaced solid-state spins, Nature Photonics 12 , 516 (2018).
- [4] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Reviews of Modern Physics 86 , 153 (2014).
- [5] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Reviews of Modern Physics 89 , 035002 (2017).
- [6] J. Rovny, S. Gopalakrishnan, A. C. B. Jayich, P. Maletinsky, E. Demler, and N. P. de Leon, New opportunities in condensed matter physics for nanoscale quantum sensors, arXiv preprint arXiv:2403.13710 (2024).
- [7] D. P. DiVincenzo, The physical implementation of quantum computation, Fortschritte der Physik: Progress of Physics 48 , 771 (2000).
- [8] B. Fauseweh, Quantum many-body simulations on digital quantum computers: State-of-the-art and future challenges, Nature Communications 15 , 2123 (2024).
- [9] J. H. Shirley, Solution of the schr¨dinger equation with a o hamiltonian periodic in time, Physical Review 138 , B979 (1965).
- [10] J. Choi, H. Zhou, H. S. Knowles, R. Landig, S. Choi, and M. D. Lukin, Robust dynamic hamiltonian engineering of many-body spin systems, Physical Review X 10 , 031002 (2020).
- [11] N. Goldman and J. Dalibard, Periodically driven quantum systems: effective hamiltonians and engineered gauge fields, Physical Review X 4 , 031027 (2014).
- [12] P. Bordia, H. L¨ uschen, U. Schneider, M. Knap, and I. Bloch, Periodically driving a many-body localized quantum system, Nature Physics 13 , 460 (2017).
- [13] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, et al. , Controlling quantum many-body dynamics in driven rydberg atom arrays, Science 371 , 1355 (2021).
- [14] S. Geier, N. Thaicharoen, C. Hainaut, T. Franz, A. Salzinger, A. Tebben, D. Grimshandl, G. Z¨rn, and u M. Weidem¨ uller, Floquet hamiltonian engineering of an isolated many-body spin system, Science 374 , 1149 (2021).
- [15] D. V. Else, B. Bauer, and C. Nayak, Floquet time crystals, Physical Review Letters 117 , 090402 (2016).
- [16] V. Khemani, A. Lazarides, R. Moessner, and S. L. Sondhi, Phase structure of driven quantum systems,

perimental setup, performed the measurements, and analyzed the data. R.F. performed the numerical simulations. C.J.W. took the data for a reference sample. M.L., R.F., E.B., S.E.E., J.C., and A.F. interpreted the results. M.L., J.C., and A.F. wrote the manuscript with inputs from all authors. All work was supervised by J.C. and A.F. Competing interests: The authors declare no competing interests. Data and materials availability: All data are available in the manuscript or the supplementary information.

Physical Review Letters 116 , 250401 (2016).

- [17] N. Y. Yao and C. Nayak, Time crystals in periodically driven systems, Physics Today 71 , 40 (2018).
- [18] J. Zhang, P. W. 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).
- [19] 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).
- [20] C. W. Thiel, T. B¨ttger, and R. Cone, Rare-earth-doped o materials for applications in quantum information storage and signal processing, Journal of Luminescence 131 , 353 (2011).
- [21] M. Lei, R. Fukumori, J. Rochman, B. Zhu, M. Endres, J. Choi, and A. Faraon, Many-body cavity quantum electrodynamics with driven inhomogeneous emitters, Nature 617 , 271 (2023).
- [22] S. Ourari, L. Dusanowski, S. P. Horvath, M. T. Uysal, /suppress C. M. Phenicie, P. Stevenson, M. Raha, S. Chen, R. J. Cava, N. P. de Leon, et al. , Indistinguishable telecom band photons from a single Er ion in the solid state, Nature 620 , 977 (2023).
- [23] B. Merkel, P. Cova Fari˜a, and A. Reiserer, Dynamical n decoupling of spin ensembles with strong anisotropic interactions, Physical Review Letters 127 , 030501 (2021).
- [24] T. B¨ttger, C. Thiel, Y. Sun, and R. Cone, Optical decoo herence and spectral diffusion at 1.5 µ m in Er 3+ :Y SiO 2 5 versus magnetic field, temperature, and Er 3+ concentration, Physical Review B 73 , 075101 (2006).
- [25] C. Thiel, W. Babbitt, and R. Cone, Optical decoherence studies of yttrium oxyorthosilicate Y 2 SiO 5 codoped with Er 3+ and Eu 3+ for optical signal processing and quantum information applications at 1.5 microns, Physical Review B 85 , 174302 (2012).
- [26] L. A. Williamson, Y.-H. Chen, and J. J. Longdell, Magneto-optic modulator with unit quantum efficiency, Physical Review Letters 113 , 203601 (2014).
- [27] M. Businger, L. Nicolas, T. S. Mejia, A. Ferrier, P. Goldner, and M. Afzelius, Non-classical correlations over 1250 modes between telecom photons and 979-nm photons stored in 171 Yb 3+ :Y SiO , Nature Communications 2 5 13 , 6438 (2022).
- [28] S. Gupta, X. Wu, H. Zhang, J. Yang, and T. Zhong, Robust millisecond coherence times of erbium electron spins, Physical Review Applied 19 , 044029 (2023).
- [29] T. Xie, J. Rochman, J. G. Bartholomew, A. Ruskuc, J. M. Kindem, I. Craiciu, C. W. Thiel, R. L. Cone,

- and A. Faraon, Characterization of Er 3+ :YVO 4 for microwave to optical transduction, Physical Review B 104 , 054111 (2021).
- [30] R. Gong, G. He, X. Gao, P. Ju, Z. Liu, B. Ye, E. A. Henriksen, T. Li, and C. Zu, Coherent dynamics of strongly interacting electronic spin defects in hexagonal boron nitride, Nature Communications 14 , 3299 (2023).
- [31] G. He, B. Ye, R. Gong, Z. Liu, K. W. Murch, N. Y. Yao, and C. Zu, Quasi-floquet prethermalization in a disordered dipolar spin ensemble in diamond, Physical Review Letters 131 , 130401 (2023).
- [32] G. He, B. Ye, R. Gong, C. Yao, Z. Liu, K. W. Murch, N. Y. Yao, and C. Zu, Experimental realization of discrete time quasi-crystals, arXiv preprint arXiv:2403.17842 (2024).
- [33] J. Choi, S. Choi, G. Kucsko, P. C. Maurer, B. J. Shields, H. Sumiya, S. Onoda, J. Isoya, E. Demler, F. Jelezko, et al. , Depolarization dynamics in a strongly interacting solid-state spin ensemble, Physical Review Letters 118 , 093601 (2017).
- [34] E. J. Davis, B. Ye, F. Machado, S. A. Meynell, W. Wu, T. Mittiga, W. Schenken, M. Joos, B. Kobrin, Y. Lyu, et al. , Probing many-body dynamics in a twodimensional dipolar spin ensemble, Nature Physics 19 , 836 (2023).
- [35] J. M. Kindem, J. G. Bartholomew, P. J. Woodburn, T. Zhong, I. Craiciu, R. L. Cone, C. W. Thiel, and A. Faraon, Characterization of 171 Yb 3+ :YVO 4 for photonic quantum technologies, Physical Review B 98 , 024404 (2018).
- [36] A. Ruskuc, C.-J. Wu, J. Rochman, J. Choi, and A. Faraon, Nuclear spin-wave quantum register for a solid-state qubit, Nature 602 , 408 (2022).
- [37] K. R. Hazzard, M. van den Worm, M. Foss-Feig, S. R. Manmana, E. G. Dalla Torre, T. Pfau, M. Kastner, and A. M. Rey, Quantum correlations and entanglement in far-from-equilibrium spin systems, Physical Review A 90 , 063622 (2014).
- [38] C. Chen, G. Bornet, M. Bintz, G. Emperauger, L. Leclerc, V. S. Liu, P. Scholl, D. Barredo, J. Hauschild, S. Chatterjee, et al. , Continuous symmetry breaking in a two-dimensional rydberg array, Nature 616 , 691 (2023).
- [39] G. Bornet, G. Emperauger, C. Chen, B. Ye, M. Block, M. Bintz, J. A. Boyd, D. Barredo, T. Comparin, F. Mezzacapo, et al. , Scalable spin squeezing in a dipolar rydberg atom array, Nature 621 , 728 (2023).
- [40] G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, et al. , Critical thermalization of a disordered dipolar spin system in diamond, Physical Review Letters 121 , 023601 (2018).
- [41] E. L. Hahn, Spin echoes, Physical Review 80 , 580 (1950).
- [42] W. K. Schenken, S. A. Meynell, F. Machado, B. Ye, C. A. McLellan, M. Joos, V. Dobrovitski, N. Y. Yao, and A. C. B. Jayich, Long-lived coherences in strongly interacting spin ensembles, arXiv preprint arXiv:2309.15444 (2023).
- [43] L. S. Martin, H. Zhou, N. T. Leitao, N. Maskara, O. Makarova, H. Gao, Q.-Z. Zhu, M. Park, M. Tyler, H. Park, et al. , Controlling local thermalization dynamics in a floquet-engineered dipolar ensemble, Physical Review Letters 130 , 210403 (2023).
- [44] J. S. Waugh, L. M. Huber, and U. Haeberlen, Approach to high-resolution nmr in solids, Physical Review Letters
- 20 , 180 (1968).
- [45] 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 601 , 531 (2022).
- [46] W. W. Ho, S. Choi, M. D. Lukin, and D. A. Abanin, Critical time crystals in dipolar systems, Physical review letters 119 , 010602 (2017).
- [47] M. Lephoto, O. Ntwaeaborwa, S. S. Pitale, H. Swart, J. Botha, and B. M. Mothudi, Synthesis and characterization of BaAl O : 2 4 Eu 2+ co-doped with different rare earth ions, Physica B: Condensed Matter 407 , 1603 (2012).
- [48] A. El Fakir, M. Sekkati, G. Schmerber, A. Belayachi, Z. Edfouf, M. Regragui, F. Cherkaoui El Moursli, Z. Sekkat, A. Dinia, A. Slaoui, et al. , Influence of rare earth (Nd and Tb) co-doping on ZnO thin films properties, Physica Status Solidi c 14 , 1700169 (2017).
- [49] C. C. Baker, A. Burdett, E. J. Friebele, D. L. Rhonehouse, W. Kim, and J. Sanghera, Rare earth co-doping for increased efficiency of resonantly pumped er-fiber lasers, Optical Materials Express 9 , 1041 (2019).
- [50] K. Singh, S. Anand, A. Pocklington, J. T. Kemp, and H. Bernien, Dual-element, two-dimensional atom array with continuous-mode operation, Physical Review X 12 , 011040 (2022)

## Supplementary Information

## Contents

| 1   | Materials and Methods                             | Materials and Methods                                                                              |
|-----|---------------------------------------------------|----------------------------------------------------------------------------------------------------|
|     | 1.1                                               | Experimental setup and device . . . . . . . . . . . . . . . . . . . . . . . . . . 2                |
|     | 1.2                                               | Experimental sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2               |
|     | 1.3                                               | Coherence measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6                  |
|     | 1.4                                               | Microscopic numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . 6                 |
| 2   | Theoretical modeling and analysis 11              | Theoretical modeling and analysis 11                                                               |
|     | 2.1                                               | Derivation of the dipolar XY Hamiltonian . . . . . . . . . . . . . . . . . . . . 11                |
|     | 2.2                                               | Tunable dipolar XXZ Hamiltonians . . . . . . . . . . . . . . . . . . . . . . . 12                  |
|     | 2.3                                               | Analysis of influences from the | Aux 〉 states . . . . . . . . . . . . . . . . . . . 13           |
| 3   | Detailed analysis for /epsilon1 -CPMG sequence 14 | Detailed analysis for /epsilon1 -CPMG sequence 14                                                  |
|     | 3.1                                               | /epsilon1 = - π/ 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15    |
|     | 3.2                                               | /epsilon1 = 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16     |
|     | 3.3                                               | Comparsion between /epsilon1 = - π/ 2 and /epsilon1 = 0 cases . . . . . . . . . . . . . . . . . 16 |
| 4   | Models I and II for the spin-echo measurement 17  | Models I and II for the spin-echo measurement 17                                                   |
|     | 4.1                                               | Two-spin case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18             |
|     | 4.2                                               | Three-spin case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20             |
| 5   | Additional experimental data                      | 22                                                                                                 |
|     | 5.1                                               | WAHUHA-echo measurement . . . . . . . . . . . . . . . . . . . . . . . . . . 22                     |
|     | 5.2                                               | Spin-locking measurement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24                |
|     | 5.3                                               | Robustness of DTC phases to initial states . . . . . . . . . . . . . . . . . . . . 25              |

## 1 Materials and Methods

## 1.1 Experimental setup and device

The experimental setup is shown in Fig. S1. For state initialization and optical readout of the ground-state spin, we address two optical transitions (labeled A and F) with a separation of 6 118 . GHz (Fig. S2a). A 984 nm laser is locked to the frequency of the A transition. A frequency sideband at the F transition is generated by an electro-optic modulator (EOM), which is driven by radio frequency (RF) pulses with a carrier frequency of 6 118 . GHz. These RF pulses are shaped by an RF switch, controlled via a transistor-transistor logic (TTL) signal from an arbitrary waveform generator (AWG), and amplified to drive the EOM, ensuring that the optical carrier power at the A transition frequency is minimized. The optical pulses addressing the A and F transitions are shaped by two sequential acousto-optic modulators (AOMs) with up to 200 MHzchirping range (we compromised the power to achieve a higher bandwidth, covering the optical inhomogeneous linewidth of 150 MHz [51]). Both AOMs are driven directly by the amplified RF signal output from the AWG. The optical light is sent to the nanophotonic device in the dilution refrigerator, where Yb ions are incorporated, and the reflected light is collected by a superconducting nanowire single photon detector (SNSPD). A gating AOM is used before the SNSPD to selectively attenuate the intense reflected input pulses (the RF drive for this AOM is not shown in Fig. S1a).

For microwave electronics, we need to coherently drive two spin transitions: f g ≈ 0 675 . GHz for the ground state manifold and f e ≈ 3 37 . GHz for the excited state manifold (Fig. S2a). For f g ( f e ), RF pulses with a carrier frequency of 200 ( 150 ) MHz are mixed with the local oscillator (LO) frequency of 875 ( 3220 ) MHz for frequency up-conversion to the target frequency. The unwanted sideband and LO components are then filtered out, and the target frequency signal is amplified. The RF signals for f g and f e are combined using a diplexer and then sent to the device through a coplanar waveguide (Fig. S1b).

Further details about the device can be found in [51].

## 1.2 Experimental sequences

The general experimental sequences are shown in Fig. S2, including state initialization, spin dynamics control, and optical readout. Fig. S2a shows the energy diagram of Yb, where both the ground and optically excited states exhibit fine structures labeled as {| 0 〉 , | 1 〉 , | Aux 〉} and {| 0 〉 e , | 1 〉 e , | Aux 〉 } e , respectively. As illustrated in Fig. S2b, the sequences begin by driving the

Figure S1: Experimental setup and device. a, The red line indicates the beam path of a 984 nm laser. The laser is locked to the frequency of the A transition. A frequency sideband at the F transition is generated by an EOM. Optical pulses are generated using AOMs. Both the EOM and AOMs are driven by gated RF sources. The light passes through a circulator to the device, and the reflected light is directed to a SNSPD for time-resolved photon counting. The blue lines indicate microwave signal delivery. RF pulses for the ground state spin transition ( f g ≈ 0 675 . GHz) and excited state spin transition ( f e ≈ 3 37 . GHz) are generated using frequency upconversion from the local oscillator signal mixed with the output of the AWG. The signals after the mixers pass through band-pass filters and amplifiers, and are then combined using a diplexer before being sent to the device. b, Scanning electron microscope image of the chip in the dilution fridge, surrounded by a coplanar waveguide. c, Scanning electron microscope image of the nanophotonic device.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-5d951daf-dfe5-40d4-8fa2-3402cc9bacff.jpg)

F transition for 100 µ s, during which the carrier frequency is swept over a chirping range of 200 MHz, with a repetition number of 10, to transfer population from the auxiliary state | Aux 〉 to the qubit manifold. Then, a sequence involving driving on transition A for 10 µ s, followed by a π pulse resonant with the excited spin transition at frequency f e , is repeated 10 times to polarize the population to | 0 〉 . After initialization, many-body spin dynamics within the qubit manifold are interrogated using spin control sequences, such as a spin-echo sequence applied to the f g transition. At the end of the sequence, a 1 µ s optical pulse is applied to the A transition

Figure S2: Experimental sequences. a, The simplified energy level diagram for a 171 Yb 3+ ion in a YVO 4 crystal. Both the ground and optically excited states exhibit fine structures labeled as {| 0 〉 , | 1 〉 , | Aux 〉} and {| 0 〉 e , | 1 〉 e , | Aux 〉 } e , respectively. Microwave spin transitions occur at frequencies f g ≈ 0 675 . GHz and f e ≈ 3 37 . GHz for the ground and excited states, respectively. The optical transitions A and F occur near a wavelength of 984 nm, with a separation of 6 118 . GHz. Note that transition A is coupled to the cavity mode used for fast spin state readout, while transition F is not coupled to the cavity mode and is driven for optical initialization. b, Experimental sequences include state initialization, spin dynamics control, and optical readout. c, Control of the average spin-spin interaction strength J within the qubit manifold {| 0 〉 , | 1 〉} (orange shaded area). Adjusting the driving amplitude on the F transition during initialization effectively controls the transfer of population from the auxiliary state | Aux 〉 to the qubit manifold {| 0 〉 , | 1 〉} (blue arrows). Subsequently, a combination of optical pulses driving the A transition and microwave pulses driving the f g transition is applied to polarize the spin state to | 0 〉 (red arrows). A higher (smaller) population in the qubit manifold corresponds to a larger (smaller) average interaction strength J .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-f1be405e-4b8e-4984-a74b-acfa3ec2b313.jpg)

Figure S3: Decoherence dominated by spin-spin interactions in spin-echo measurements. a, Normalized coherence decay of the spin ensemble as a function of evolution time in the spin-echo sequence. The early-time decay within 1 µs accelerates with increasing J , achieved by increasing the driving amplitude on the F transition during state initialization. Error bars represent the standard deviation of the experimental measurements (see Sec. 1.3 for details), while the solid lines are simple exponential fits used to extract the 1 /e decoherence times. b, The fitted decoherence times monotonically decrease with increasing drive amplitude on the F transition, implying that a larger population within the qubit manifold leads to stronger spin-spin interactions and thus faster decoherence.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-fc4e9037-8f4c-430f-a027-4cea0e308b35.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-203547fb-4d23-42e9-ab5b-2a89ffb4f2f5.jpg)

## for reading out the population in | 1 〉 .

As shown in Fig. S2c, we can adjust the average interaction strength J between spins by tuning the population within the qubit manifold. This adjustment is achieved by varying the driving amplitude on the F transition during initialization, which controls the amount of population transferred from the auxiliary state | Aux 〉 to the qubit manifold {| 0 〉 , | 1 〉} . Specifically, we experimentally confirm that a higher (smaller) population in the qubit manifold leads to a larger (smaller) average interaction strength J , as evidenced by the spin density-dependent decoherence rates observed in the spin-echo sequence (Fig. S3a). The early-time decay rate within the first 1 µ s is governed by the average spin-spin interaction strength J , which exhibits a monotonic scaling with the driving amplitude on the F transition (Fig. S3b). This arises from the accelerated pump-out rate from the | Aux 〉 state due to strong driving on transition F.

## 1.3 Coherence measurement

For all the experiments presented in the main text, the initial state is prepared as a globally polarized state along the y -axis. The ensemble coherence of the many-body spin system can be understood as the mean coherence of individual spins, quantified by examining the residual spin polarization along the y -axis at the end of spin dynamics interrogation (Fig. S4a). In other words, measuring coherence is equivalent to measuring the polarization along the y -axis. Note that the π/ 2 analyzer pulse converts the y -axis polarization into z -axis polarization, which can be read out using photoluminescence signals C . By sweeping the phase angle, θ , of the analyzer pulse (at a fixed interrogation time), we obtain sinusoidal oscillations of photoluminescence as a function of θ and fit them with C θ ( ) = C amp cos θ + C offset (Fig. S4b). The coherence, ˜ C , is then calculated as ˜ = C C max -C min C max + C min = C amp C offset , and the corresponding error bar is extracted from the self-covariance of the fit (Fig. S4c).

Figure S4: Coherence characterization in spin-echo measurements. a, The spin-echo sequence with a final π/ 2 analyzer pulse with a variable phase angle θ . The phase angle θ defines the rotation axis of the π/ 2 pulse relative to the initial π/ 2 pulse. While maintaining a fixed free evolution time τ , θ is varied from 0 to 2 π to quantify the residual coherence of the spin ensemble. b, Photoluminescence signals from the spin system exhibit a sinusoidal oscillation as θ is swept, from which coherence is extracted as the contrast of the oscillation. c, Coherence decay profiles as a function of evolution time τ for the cases of large J ≈ 2 π × 0 35 . MHz and small J ≈ 2 π × 0 19 . MHz.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-7675918c-81c0-43fe-b2e6-3927e7493045.jpg)

## 1.4 Microscopic numerical simulation

Here, we present the details of the numerical simulations used throughout this work. To set up the spin system in simulations, we begin by generating the lattice structure of YVO 4 (Fig. S5a).

Wecenter the lattice on a Y ion and replace it with a Yb ion, then the lattice size is appropriately extended to accommodate the doping concentration of Yb ions for N ions. These N ions are randomly substituted into the Y sites (including the central one) with random spin frequency detuning ∆ i . The values of ∆ i are randomly sampled from a Lorentzian distribution with a full width at half maximum W . This defines the on-site disorder Hamiltonian, ˆ H dis , given by

$$\hat { H } _ { \text{dis} } = \sum _ { i = 1 } ^ { N } \Delta _ { i } \hat { S } _ { z } ^ { i }. & & ( S 1 ) \\ \text{nin-s} \text{interactions the native interaction strength } I \cdot \text{between}$$

Regarding the modeling of spin-spin interactions, the pairwise interaction strength J ij between each pair of Yb ions is calculated based on their dipolar XY interaction Hamiltonian, ˆ H int :

$$\hat { H } _ { \text{int} } = \sum _ { i j, \ i > j } ^ { N } J _ { i j } ( \hat { S } _ { x } ^ { i } \hat { S } _ { x } ^ { j } + \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } ), \\ \mu _ { 0 } \mu _ { R } ^ { 2 } g _ { \text{II} } ^ { 2 } \,. \quad \text{e}$$

$$J _ { i j } = - \frac { \stackrel { \psi \cdot \cdot \cdot \cdot } { \mu _ { 0 } } \mu _ { B } ^ { 2 } g _ { \| } ^ { 2 } } { 8 \pi r _ { i j } ^ { 3 } } ( 3 z _ { i j } ^ { 2 } - 1 ).$$

Here, µ 0 is the magnetic permeability of free space, g ‖ = 6 08 . is the g-factor along the crystal c -axis [52] (defining the c -axis to be the z -axis), r ij is the distance between ions i and j , and z ij ≡ | /vector r ij · /vector e z | /r ij is the z -directional cosine between the ions ( - ≤ 1 z ij ≤ 1 ), where /vector e z is the unit vector along the z -axis. Since J ij is dependent on both distance and orientation between each spin pair, we define the average spin-spin interaction J as follows:

$$J \equiv \frac { \mu _ { 0 } \mu _ { B } ^ { 2 } g _ { \| } ^ { 2 } } { 4 \pi } \frac { 1 } { \left \langle r _ { i j } \right \rangle ^ { 3 } } = 2 \pi \times 4 8 0 \, \text{MHz} \cdot \text{nm} ^ { 3 } \times \frac { 1 } { \left \langle r _ { i j } \right \rangle ^ { 3 } }, \\ \psi \sim \dots \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \color \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col. \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col.: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \col: \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colot \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colots \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \colon \, \text{.}$$

where 〈 r ij 〉 is the average nearest-neighbor distance in units of nm, and we have used the condition z ij = ± 1 to maximize the interaction strength in different orientations. For a given crystal structure, the average interaction strength is linearly proportional to the spin density n s , i.e., J ∝ n s . For YVO 4 crystal, the lattice parameters are a = 7 119 . ˚ A , c = 6 290 . ˚ . A Considering that there are 4 ions per unit cell, the relation between the average distance, 〈 r ij 〉 in units of nm, and the spin concentration, n s in units of ppm, can be expressed as 〈 r ij 〉 = 0 629 . (4 n s ) 1 / 3 × 10 -2 nm. Since the g-factor of a Yb ion is three times larger than that of an electron, we expect the interaction strength in our system to be approximately 9 times larger compared to the case of an NV center, where J = 2 π × 52 MHz · nm 3 × 1 〈 r ij 〉 3 [53].

Given a random realization of a disordered spin system composed of N Yb ions, the manybody spin states can be fully characterized by solving for noiseless dynamics under the system

a

Figure S5: Numerical simulation results of the spin-echo decoherence dynamics. a, The crystal structure of YVO . 4 b, System size scaling of the spin-echo decoherence profiles for different numbers of simulated spins from N = 2 to N = 10 ions. The decay profiles start to converge approximately when N &gt; 8 ions. c, Effects of on-site disorder on spin-echo decoherence dynamics, simulated for different disorder strengths W . After the early-time transient ( t &lt; 1 µ s) where the decay rate is independent of on-site disorder, spin decoherence slows down, indicating a crossover from an interaction-dominated to a disorder-dominated regime. d, Effects of initial spin polarization, η pol , on spin-echo decoherence dynamics. e, Normalized coherence for the decay time traces shown in d . Upon rescaling, the decay profiles overlap for different polarization values of η pol , indicating that the decoherence dynamics are independent of initial spin polarization. 8

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-99d5b771-1f25-426d-a86d-af158af5f7f8.jpg)

Hamiltonian ˆ H = ˆ H dis + ˆ H int , incorporating control pulse sequences such as spin-echo or Ramsey sequences. To mitigate numerical artifacts stemming from boundary effects in finitesize simulations, only the final state of the central Yb ion is considered for readout. This process is repeated across multiple Monte Carlo runs, where the positions and on-site detunings of the ions are randomized each time, to capture ensemble-averaged dynamics.

To determine the number of ions needed to simulate the system accurately, a convergence test on system size N is conducted (Fig. S5b). We find that the dynamics start to converge for system sizes larger than 8 ions, as adding ions farther away does not significantly contribute to the dynamics of the central readout spin. Henceforth, we fix the system size to N = 9 ions in Fig. S5c-e, as well as the simulation results presented in the main text.

Based on the many-body simulation of a disordered spin ensemble, we extract the average interaction strength, J , and on-site disorder, W , by comparing the simulation results to the experimentally observed dynamics. First, the early-time decay rate within 1 µs does not depend on on-site disorder strength W ; it is only governed by the spin-spin interaction J (Fig. S3c). This allows us to determine the average interaction strength, where the only unknown parameter is n s , the concentration of Yb ions. By matching the decay rate at short times ( t &lt; 1 µs ) from the simulation to that of the experiment, we can extract the ion concentration within the qubit manifolds for different cases, yielding 25 ppm for the small J case and 46 ppm for the large J case (see Fig. 2a in the main text). The extracted effective spin densities within the qubit manifold are estimated to be lower than an independently measured concentration of 86 ppm using secondary ion mass spectrometry, which reflects the total density of Yb ions [54]. This difference is attributed to the remaining untransferred population residing in the doubly degenerate | Aux 〉 states outside the qubit manifold. We then extract W = 2 π × 0 65 . MHz, independent of the interaction strength J , by finding the best agreement in the experiment and simulation under the Ramsey sequence (see Fig. 2b in the main text). As shown in Fig. S5c, the presence of on-site disorder slows down the decoherence dynamics at late times in the spin-echo sequence. We successfully match this extended-time behavior using the disorder strength W , extracted from Ramsey data (Fig. S6).

These calibrated interaction and disorder strengths reproduced all the experimental data across all cases presented in the main text figures. This strongly suggests that the decoherence dynamics of the spin ensemble are driven by the intrinsic spin-spin interactions and on-site disorder described by the closed many-body Hamiltonian, without additional dephasing due to extrinsic coupling to external environments.

Figure S6: Comparison of long-time decoherence dynamics between experiment and simulation under the spin-echo dynamics. The numerical simulations with calibrated parameters (solid lines) show excellent agreement with the experiments (markers) over extended long times for both the small J (red) and large J (blue) cases.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-8f8afc94-6c97-4495-bb79-e10f9f33c9e4.jpg)

We also investigate whether the initial spin distribution within the qubit manifold influences the resulting decoherence processes. Specifically, as described earlier, we employ both optical transition A and microwave transition f e to initialize the Yb ion's population into the | 0 〉 state. We can quantify the fidelity of this process as initial polarized rate η pol ≡ p 0 p 0 + p 1 , where p 0 and p 1 correspond to the populations in | 0 〉 and | 1 〉 after initialization, respectively. The relation between η pol and the measured coherence ˜ C (Sec. 1.3) is ˜ C = 2 η pol -1 . According to the experimental data shown in Fig. S4c near the τ = 0 point, we estimate that η pol ≈ 0 75 . . From numerical simulations investigating the dependence of decoherence dynamics on initial polarized rates, we find that the decoherence profiles do not depend on the value of η pol , except for a reduction in the contrast of the coherence signal due to imperfect initialization (Fig. S5d). We confirm that all cases overlap when each imperfect case is rescaled by a factor of 1 / (2 η pol -1) or equivalently by the maximum contrast at τ = 0 (Fig. S5e). The independence of the initial polarized rate arises because, despite individual spins being randomly oriented along either the + y or -y directions, the dipolar XY Hamiltonian (Eq. S2) induces dephasing dynamics through terms like ˆ S S i x ˆ j x , which flip ions i and j , regardless of whether they are initially aligned parallel or antiparallel along the y -axis.

## 2 Theoretical modeling and analysis

## 2.1 Derivation of the dipolar XY Hamiltonian

Here, we provide a theoretical analysis for constructing the dipolar XY Hamiltonian in our system (Eq. S2). We start with the magnetic dipole-dipole interaction and derive an effective spin-spin interaction within the qubit manifold. The ground-state qubit states {| 0 〉 | , 1 〉} at zero magnetic field can be expressed as:

$$| 0 \rangle = \frac { 1 } { \sqrt { 2 } } ( | \uparrow \Downarrow \rangle - | \downarrow \uparrow \rangle )$$

$$| 1 \rangle = \frac { 1 } { \sqrt { 2 } } ( | \uparrow \Downarrow \rangle + | \downarrow \uparrow \rangle )$$

where {|↑〉 , |↓〉} are electron spins and {|⇑〉 , |⇓〉} are nuclear spins associated with a Yb ion [52]. The magnetic dipole-dipole interaction Hamiltonian between the spins of ions i and j is described as follows:

$$\hat { H } _ { d d } ^ { i j } = - \frac { \mu _ { 0 } } { 4 \pi r _ { i j } ^ { 3 } } ( 3 ( \vec { \mu } _ { i } \cdot \vec { e } _ { r _ { i j } } ) ( \vec { \mu } _ { j } \cdot \vec { e } _ { r _ { i j } } ) - \vec { \mu } _ { i } \cdot \vec { \mu } _ { j } )$$

where /vector e r ij ≡ /vector r ij /r ij is the normalized unit distance vector from ion i to ion j , and /vector µ j is the magnetic dipole moment operator for ion j , given as

$$\text{ment operator for ion } \jmath, \text{given as} \\ \vec { \mu } _ { j } = \mu _ { B } \mathfrak { g } \cdot \mathbf S _ { 0 } ^ { j } = \mu _ { B } \begin{pmatrix} g _ { \perp } \hat { S } _ { 0, x } ^ { j } \\ g _ { \perp } \hat { S } _ { 0, y } ^ { j } \\ g _ { \| } \hat { S } _ { 0, z } ^ { j } \end{pmatrix} & & ( S 8 ) \\ \text{hr magneton, } \mathfrak { g } = \text{diag} ( g _ { \perp }, g _ { \perp }, g _ { \| } ) \text{ is the anisotropic g-factor tensor with}$$

Since we are interested in deriving the effective interaction Hamiltonian within the qubit manifold, we first calculate the matrix elements of the magnetic dipole moment with respect to the qubit {| 0 〉 , | 1 〉} basis states:

where µ B is the Bohr magneton, g = diag ( g ⊥ , g ⊥ , g ‖ ) is the anisotropic g-factor tensor with g ⊥ = 0 85 . and g ‖ = 6 08 . , and S j 0 = ( ˆ S j 0 ,x , S ˆ j 0 ,y , S ˆ j 0 ,z ) T is the vectorized spin operator defined in Yb's electron spin basis {|↑〉 , |↓〉} (which is different from the {| 0 〉 , | 1 〉} basis).

$$\langle 1 | _ { j } \, \vec { \mu } _ { j } \, | 0 \rangle _ { j } & = \langle 0 | _ { j } \, \vec { \mu } _ { j } \, | 1 \rangle _ { j } = \mu _ { B } \begin{pmatrix} g _ { \perp } \, \langle 0 | _ { j } \, \hat { S } _ { 0, x } ^ { j } \, | 1 \rangle _ { j } \\ g _ { \perp } \, \langle 0 | _ { j } \, \hat { S } _ { 0, y } ^ { j } \, | 1 \rangle _ { j } \\ g _ { \| } \, \langle 0 | _ { j } \, \hat { S } _ { 0, z } ^ { j } \, | 1 \rangle _ { j } \end{pmatrix} = \frac { 1 } { 2 } \mu _ { B } \begin{pmatrix} 0 \\ 0 \\ g _ { \| } \end{pmatrix} \\ \langle 0 | _ { j } \, \vec { \mu } _ { j } \, | 0 \rangle _ { j } & = \langle 1 | _ { j } \, \vec { \mu } _ { j } \, | 1 \rangle _ { j } = 0. \\ \quad \dots$$

Note that both | 0 〉 and | 1 〉 states have no magnetic moments. This indicates that the qubit states are first-order insensitive to external magnetic field, defining our 'clock' transition with a zero first-order Zeeman (ZEFOZ) shift.

We then consider the matrix elements of the dipolar interaction Hamiltonian (Eq. S7) in a two-qubit basis, {| 00 〉 , | 01 〉 , | 10 〉 , | 11 〉} , and find that the ZEFOZ nature of our qubit leads to the following form:

$$\hat { H } _ { d d } ^ { i j } = \frac { J _ { i j } } { 2 } ( \hat { S } _ { + } ^ { i } \hat { S } _ { - } ^ { j } + \hat { S } _ { - } ^ { i } \hat { S } _ { + } ^ { j } + \hat { S } _ { + } ^ { i } \hat { S } _ { + } ^ { j } + \hat { S } _ { - } ^ { i } \hat { S } _ { - } ^ { j } ).$$

Here, ˆ S j ± = ˆ S j x ± iS ˆ j y are the raising and lowering spin operators for ion j , defined in the qubit {| 0 〉 , | 1 〉} basis, and J ij is defined in Eq. S3. The dipole-dipole interaction described above includes flip-flop, flop-flip, flip-flip, and flop-flop terms. However, the last two terms can be neglected due to energetic suppression, as the qubit transition frequency, f g ≈ 675 MHz, is much larger than the interaction strength, J ij / π 2 ∼ 0 1 . MHz. This approximation is known as the secular approximation, which leads to the following pure spin-exchange pairwise interaction:

$$\hat { H } _ { d d } ^ { i j } \approx \frac { J _ { i j } } { 2 } ( \hat { S } _ { + } ^ { i } \hat { S } _ { - } ^ { j } + \hat { S } _ { - } ^ { i } \hat { S } _ { + } ^ { j } ) = J _ { i j } ( \hat { S } _ { x } ^ { i } \hat { S } _ { x } ^ { j } + \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } ),$$

consistent with the interaction Hamiltonian between ions i and j , as presented in Eq. S2.

## 2.2 Tunable dipolar XXZ Hamiltonians

In the absence of an external magnetic field, our system Hamiltonian is described by the dipolar XY model, which contains only spin-exchange interactions (Eq. S2). However, applying a small magnetic field along the z -axis allows us to achieve a tunable spin model that includes Ising interactions, known as the XXZ model. Specifically, with a small magnetic field strength B , the original eigenstates | 0 〉 and | 1 〉 are perturbed as follows [55]:

Here, { ∣ ∣ ˜ 0 〉 , ∣ ∣ ˜ 1 〉 } are the perturbed eigenstates with an admixture ratio α B ≡ g µ ‖ B B/ω , and ω = 2 πf g = 2 π × 675 MHz is the angular qubit transition frequency. Under these perturbed eigenstates, the system's original Hamiltonian, shown in Eq. S2, can be modified as

$$0 \rangle \, \text{and} \, | 1 \rangle \, & \text{are perturbed as follows} \, [ \text{55} \colon \\ | \tilde { 0 } \rangle & = | 0 \rangle - \frac { \alpha _ { B } } { 2 } \, | 1 \rangle & ( S 1 1 ) \\ | \tilde { 1 } \rangle & = | 1 \rangle + \frac { \alpha _ { B } } { 2 } \, | 0 \rangle & ( S 1 2 ) \\ \text{:rurbed eigenstates with an admixture ratio } \, \alpha _ { B } \, \equiv \, g _ { \| } \mu _ { B } B / \omega, \, \text{and} \\ \text{Iz is the angular qubit transition frequency. } \, \text{Under these perturbed}$$

$$\text{system's original Hamiltonian, shown in $Eq.$ SZ$, can be modified as} \\ \hat { H } _ { \text{int} } & = \left ( 1 - \frac { \alpha _ { B } ^ { 2 } } { 2 } \right ) \hat { H } _ { \text{exchange} } + 2 \alpha _ { B } ^ { 2 } \hat { H } _ { \text{Iming} } ^ { z }$$

where the spin-exchange Hamiltonian, ˆ H exchange , and the Ising Hamiltonian, ˆ H z Ising are given by

$$\hat { H } _ { \text{exchange} } = \sum _ { i j,$$

We have used the spin operators for the new perturbed basis, { ∣ ∣ ˜ 0 〉 , ∣ ∣ ˜ 1 〉 } , and applied the secular approximation to preserve the energy-conserving term. From the new Hamiltonian in Eq. S13, one can observe that applying an external magnetic field tunes the magnetic field-dependent admixture ratio, α B , thereby controlling the relative strengths of the spin-exchange interaction, ˆ H exchange , and the Ising interaction ˆ H z Ising .

$$\hat { H } _ { \text{ising} } ^ { z } = \sum _ {$$

## 2.3 Analysis of influences from the | Aux 〉 states

In the main text, we consider only the effective interactions within the qubit manifold, despite the presence of the doubly-degenerate | Aux 〉 state in the ground state. Using the electron and nuclear spin notations introduced earlier, the two degenerate states are expressed as |↑⇑〉 and |↓⇓〉 . The populations in these | Aux 〉 states can also interact with the spins in the qubit manifold. In this section, we will analyze the influence of the | Aux 〉 state on the qubit manifold.

$$\text{it states. For instance, between} \left | \uparrow \uparrow \right$$

Specifically, the auxiliary states can interact with the qubit states via the full dipole-dipole Hamiltonian given in Eq. S7. Following a similar procedure presented in Sec. 2.1, we first calculate the matrix element of the magnetic moment operator between | Aux 〉 and one of the qubit states. For instance, between |↑⇑〉 and | 0 〉 , we have:

Next, we consider the spin-exchange process between ion i at | Aux 〉 and ion j at | 0 〉 , such as |↑⇑〉 i | 0 〉 j ↔ | 〉 |↑⇑〉 0 i j . The corresponding relevant matrix element of the Hamiltonian can be computed as

$$\langle 0 | _ { i } \, \langle \uparrow \uparrow | _$$

where J ij is the spin-exchange interaction strength within the qubit manifold defined in Eq. S3. Note that the spin-exchange interaction strength between | Aux 〉 and the qubit states is significantly smaller by a factor of ( g ⊥ /g ‖ ) 2 ≈ 0 02 . . Consequently, using Fermi's golden rule, the transition rate from | 0 〉 to | Aux 〉 can be compared against that from | 0 〉 to | 1 〉 :

$$\begin{array} {$$

This implies that the population exchange between the ions in the | Aux 〉 state and the ions in the qubit manifold occurs at a rate that is approximately 10 4 times slower than the population exchange within the qubit manifold. Additionally, the larger on-site disorder of the | Aux 〉 states, due to their sensitivity to external fields, further prevents this flip-flop process between | Aux 〉 and the qubit manifold, ruling out its contribution to the late-time slow decoherence observed in the spin-echo experiment.

## 3 Detailed analysis for /epsilon1 -CPMG sequence

Figure S7: Interaction-picture-based toggling-frame transformation of the ˆ S z operator for /epsilon1 -CPMG sequences. a, /epsilon1 = -π/ 2 is a preferable choice for decoupling spin-spin interactions in a strongly interacting spin system. b, /epsilon1 = 0 corresponds to the conventional CPMG sequence, which decouples on-site disorder through periodic π pulses. These base sequences are repeated in time stroboscopically with fixed pulse spacing τ in dynamic Hamiltonian engineering.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-6ae6b296-7606-4b20-8437-e2ffc5857289.jpg)

The /epsilon1 -CPMG sequence is employed to verify the presence of strong spin-spin interactions in our system due to its distinct sensitivity to interaction and on-site disorder (see Fig. 2 of the main text). When /epsilon1 = 0 , it is identical to the standard CPMG sequence, effectively de-

coupling time-independent on-site disorder while leaving pairwise spin-spin interactions unchanged (Fig. S7b). In contrast, when /epsilon1 = ± π/ 2 , the initial spin states polarized along the y -axis can be protected against both on-site disorder and spin-spin interactions to leading order. This is consistent with experimental observations showing that coherence at /epsilon1 = ± π/ 2 is higher than at /epsilon1 = 0 (see Fig. 2c of the main text).

In the following, we will derive the zeroth-order and first-order Hamiltonians for /epsilon1 -CPMG sequences at /epsilon1 = 0 and -π/ 2 via average Hamiltonian theory [56], starting from the parent system Hamiltonian ˆ H :

$$\hat { H } & = \hat { H } _ { \text{dis} } + \hat { H } _ { \text{int} } = \sum _ { i = 1 } ^ { N } \Delta _ { i } \hat { S } _ { z } ^ { i } + \sum _ { i j, \, i > j } ^ { N } J _ { i j } ( \hat { S } _ { x } ^ { i } \hat { S } _ { x } ^ { j } + \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } ) \\ \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \quad \. \.$$

The effective average Hamiltonians, ˆ H eff , dynamically engineered from ˆ H , can be analytically computed order-by-order using the framework based on the interaction-picture-based togglingframe transformation of the ˆ S z operator [56]. Specifically, the system Hamiltonian ˆ H transforms into an effective Hamiltonian ˆ H eff as follows:

$$\hat { H } \to \hat { H } _ { \text{eff} } = \sum _ { m = 0 } ^ { \infty } \hat { H } _ { m } \\ \vdots \quad \ddots \ \cdot \quad \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot$$

where ˆ H m denotes the average Hamiltonian of order m determined by a control pulse sequence.

$$\begin{matrix} 3. 1 & \epsilon = - \pi / 2 \\ \dashrightarrow & \dots \end{matrix}.$$

Here we consider the -π 2 -CPMG sequence with the base control sequence depicted in Fig. S7a. The zeroth-order Hamiltonian, ˆ H 0 , is given as

$$\hat { H } _ { 0 } & = \frac { 1 } { 2 } \sum _ { i j, \, i > j } ^ { N } J _ { i j } ( \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } + \vec { S } ^ { i } \cdot \vec { S } ^ { j } ) = \frac { 1 } { 2 } \hat { H } _ { \text{ising} } ^ { y } + \frac { 1 } { 2 } \hat { H } _ { \text{Heis} }, \\ \wedge \,. \quad & = \pi \, \quad \wedge \,. \quad \,. \quad \,. \quad \,. \quad \,.$$

where ˆ H y Ising = ∑ N ij, i&gt;j J ij ˆ S S i y ˆ j y is the Ising interaction along the y -axis, and ˆ H Heis = ∑ N ij, i&gt;j J ij /vector S i · /vector S j is the Heisenberg interaction. Note that the initially polarized spin state along the y -axis is the eigenstate of ˆ H 0 , resulting in the preservation of coherence.

The next first-order Hamiltonian, ˆ H 1 , is given as

/negationslash

$$\tilde { H } _ { 1 } & = \frac { \tau } { 4 } \left [ \sum _ { i } ^ { N } \Delta _ { i } ^ { 2 } \hat { S } _ { y } ^ { i } + \sum _ { i j, \, i \neq j } ^ { N } \Delta _ { i } J _ { i j } \left ( 2 \hat { S } _ { z } ^ { i } \hat { S } _ { y } ^ { j } - \hat { S } _ { y } ^ { i } \hat { S } _ { z } ^ { j } \right ) \right ], \\ & \cdot \tilde { \tau }$$

which gives rise to the dephasing of the polarized spin ensemble. In principle, symmetrizing the pulse sequence can eliminate all odd-order average Hamiltonians [56], potentially extending the coherence time. We leave it for future work.

## 3.2 /epsilon1 = 0

Here we consider the /epsilon1 = 0 case corresponding to the conventional CPMG sequence with the base control sequence depicted in Fig. S7b. Its zeroth-order Hamiltonian, ˆ H 0 , is given as

$$\hat { H } _ { 0 } = \sum _ { i j, \, i > j } ^ { N } J _ { i j } ( \hat { S } _ { x } ^ { i } \hat { S } _ { x } ^ { j } + \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } ) = \hat { H } _ { \text{Heis} } - \hat { H } _ { \text{Iss} } ^ { z }$$

where ˆ H z Ising is the Ising interaction along the z -axis. The initially polarized spin state along the y -axis undergoes dephasing with ˆ H 0 because the ˆ H z Ising term provides an interaction-induced effective magnetic field along the z -axis.

The next first-order Hamiltonian, ˆ H 1 = 0 . In fact, all odd-order average Hamiltonians of the CPMG sequence are zero due to the mirror symmetry of the toggling-frame transformations [56].

## 3.3 Comparsion between /epsilon1 = -π/ 2 and /epsilon1 = 0 cases

The time traces of spin ensemble coherence for the /epsilon1 -CPMG sequences at /epsilon1 = 0 and /epsilon1 = -π/ 2 are shown in Fig. S8a,b. Consistent with the average Hamiltonian analysis, the /epsilon1 -CPMG sequence with /epsilon1 = -π/ 2 exhibits a significant enhancement in coherence, attributed to the zerothorder Hamiltonian supporting the initially polarized state as an eigenstate. It is noteworthy that the rapid driving of the base sequence with a shorter periodicity leads to a much longer coherence time, which corroborates the effectiveness of the zeroth-order Hamiltonian in dominating the engineered many-body dynamics in the fast driving regime. Furthermore, when the rotation angle offset, /epsilon1 , is swept at a fixed sequence repetition number k , the /epsilon1 -CPMG sequence reveals an unconventional profile (Fig. S8c). This profile is in stark contrast to the behavior expected in a non-interacting spin system, providing further evidence of the strong spin-spin interactions within our spin system (see Fig. 2c of the main text).

Figure S8: Coherence characterization in /epsilon1 -CPMG measurements. a, b, We compare coherence decay profiles as a function of total interrogation time T for two different /epsilon1 -CPMG measurements at /epsilon1 ≈ 0 (yellow) and /epsilon1 ≈ -π/ 2 (purple). For both cases, the corresponding base sequence has a periodicity of τ + t p , where τ is the free evolution period and t p is the pulse duration, and is repeated k times to advance in time (see the insets). We measure coherence at stroboscopic times T = k τ ( + t p ) . a, τ/ 2 = 200 ns, t p = 45 ns for /epsilon1 ≈ 0 ; t p = 25 ns for /epsilon1 ≈ -π/ 2 . b, τ/ 2 = 50 ns, t p = 43 ns for /epsilon1 ≈ 0 ; t p = 22 ns for /epsilon1 ≈ -π/ 2 . Experimental data are obtained under the large J condition, and solid lines represent phenomenological stretched exponential fits. To facilitate comparison between the two different /epsilon1 cases, the coherence decay profiles are normalized by their respective maximum coherence values. c, Dependence of /epsilon1 -CPMGsequence coherence on /epsilon1 for the long τ (blue) and short τ (red) cases. When sweeping /epsilon1 , fixed base sequence parameters of ( τ/ 2 = 200 ns , k = 5; blue ) and ( τ/ 2 = 50 ns , k = 20; red ) are used respectively. Here, all experimental data across different /epsilon1 values are globally normalized by the initial state polarization, η pol (see Sec. 1.4 and Fig. S5d,e). The solid lines denote numerical simulation results using the experimentally calibrated system parameters. The simulated results were globally rescaled to facilitate comparison with the experimental data.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-b3623a09-af24-482b-9d53-6b50f2cf3248.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-d9dc6db9-cafe-4860-976c-c86f0ddc91a4.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-c8cbcf61-e98a-416d-be25-83a343fe9a14.jpg)

## 4 Models I and II for the spin-echo measurement

In this section, we provide a more detailed explanation of Models I and II mentioned in Fig. 3 of the main text.

In Model I, we simulate the spin-echo signal governed by only a single pair of resonant spins with the largest pairwise interaction strength in a given realization, denoted as J max . The spinexchange interaction induces periodic entanglement and disentanglement dynamics between this pair, resulting in coherence oscillations at a rate of J max / 2 . Concretely, after the initial π/ 2 pulse along the x -axis, the quantum state is prepared as | ψ (0) 〉 = 1 2 ( 0 | 〉 -i | 1 ) 〉 ⊗ | 〉 -( 0 i | 1 ) 〉 , which then evolves as:

$$| \psi ( \tau ) \rangle = \frac { 1 } { 2 } \left [ | 0 0 \rangle - | 1 1 \rangle - i e ^ { - i J _ { \max } \tau / 2 } ( | 0 1 \rangle + | 1 0 \rangle ) \right ]$$

under the two-spin Hamiltonian ˆ H = J max 2 ( ˆ S S + ˆ -+ ˆ S S -ˆ + ) . As a result, the corresponding coherence oscillates sinusoidally as cos( J max τ/ 2) . Note that J max is a random number determined by a stochastic realization of spin positions in the numerical simulation. Through repeated Monte Carlo sampling of different spin configurations, we calculate the ensemble average of coherence as a function of time, revealing rapid decay. Despite the simplicity of the model, we find that the early-time decay rate of Model I is consistent with that of the full simulation carried out in Model III.

In Model II, we expand upon Model I and simulate many-body dynamics involving more than two spins without disorder. Each spin can now become entangled with multiple neighboring spins simultaneously, undergoing intricate dephasing dynamics. In each simulation run, we calculate the coherence of the center spin and perform Monte Carlo averaging over different realizations of spin configurations. As seen in Fig. 3 of the main text, it is observed that adding more spins slows down the late-time decay compared to Model I, thereby improving the simulation's agreement with the experimental data.

To gain a better understanding of why Model II has slower dynamics compared to Model I at later times, we will derive and discuss two toy models-the two-spin and three-spin cases-in the following subsections.

## 4.1 Two-spin case

For a pair of spins with interaction strength J and relative detuning ∆ , we can derive the average spin polarization along the ˆ y -axis at time τ , P τ ( ) , under the spin-echo sequence as follows:

$$P ( \tau ) = \frac { \Delta ^ { 2 } } { \Delta ^ { 2 } + J ^ { 2 } } + \frac { J ^ { 2 } } { \Delta ^ { 2 } + J ^ { 2 } } \cos \left ( \frac { \sqrt { \Delta ^ { 2 } + J ^ { 2 } } } { 2 } \tau \right ).$$

Recall that the polarization signal represents the spin coherence. For the resonant spins with ∆ = 0 , as in Model I, P τ ( ) = cos( Jτ/ 2) . This implies that the two spins oscillate between the ˆ y and -ˆ y directions with a frequency of J/ 2 , similar to Rabi oscillation dynamics (blue curve, Fig. S9a). Ensemble averaging over a distribution of J leads to a decay in spin polarization, resulting from the incoherent averaging of sinusoidal oscillations with different frequencies (blue curve, Fig. S9b).

/negationslash

However, when a non-zero detuning ∆ = 0 is introduced, the spin polarization signal develops a non-zero positive offset of ∆ 2 ∆ + 2 J 2 with a reduced contrast of J 2 ∆ + 2 J 2 (red curve, Fig. S9a, and see Eq. S24). Similar to the earlier case, ensemble averaging over random distributions of

Figure S9: Two-spin model with interaction strength J and relative detuning ∆ . Two spins are initially polarized along the y -axis, and coherence is defined as the remnant polarization along the y -axis after evolving for a duration τ . a, Single realization (shot) dynamics of coherence for a fixed J = 2 π × 0 3 . MHz, but with ∆ = 0 (blue) and ∆ = J (red). b, Ensemble-averaged dynamics of coherence for different J and ∆ values. The distribution of J is assumed to follow a Gaussian distribution with a full width at half maximum (FWHM) of J FHWM = 2 π × 0 65 . MHz. Two cases are considered for a distribution of relative detuning ∆ , which is also assumed to follow a Gaussian distribution with a FWHM of W : one with no disorder ( W = 0 , blue) and another with disorder ( W = 2 π × 0 65 . MHz, red). c, Longer-time evolution is simulated for b . The presence of on-site disorder gives rise to a non-zero offset in the late-time saturated coherence.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-b5ef2eb1-5fb9-4c51-9c6c-8dff3339bffa.jpg)

J and ∆ leads to a decay in spin polarization (red curve, Fig. S9b). Notably, we observe that the early-time decay rate at small τ remains the same between the no-disorder and disordered cases (Fig. S9b). Specifically, the ensemble-averaged coherence decay rate at early times is given by

$$\xi \cdot \mathcal { V } \mathcal { P } \colon \text{pem any, and an uniou-sauveau} \xi \colon \text{conference away raw at any } \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \, \alpha \} \, \max \, \xi \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal{ V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal{\v } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal {V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { W } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { U } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { U } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { V } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { U } \, \mathcal { V } \, \mathcal { P } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { Y } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { V } \, \mathcal { T } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { T } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { T } \, \mathcal { T } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P } \, \mathcal { P$$

which depends only on 〈 J 2 〉 , the expectation value of the interaction strength squared. Here, 〈 ... 〉 denotes ensemble averaging.

At late times, however, the presence of on-site disorder gives rise to a non-zero offset in

the late-time saturated coherence, which does not decay in this two-spin model (red curve, Fig. S9c). This suggests that involving more ions is necessary to explain our experimental observation of a slow decay emerging in the late-time decoherence profile, as numerically demonstrated in Fig. S5b.

## 4.2 Three-spin case

Here we introduce a third spin to analyze how it affects the coherence dynamics. First, we examine a case with no on-site disorder for simplicity, where the three-spin Hamiltonian is given by

$$\hat { H } = J _ { 0 } ( \hat { S } _ { x } ^ { 1 } \hat { S } _ { x } ^ { 2 } + \hat { S } _ { y } ^ { 1 } \hat { S } _ { y } ^ { 2 } ) + J _ { 1 } ( \hat { S } _ { x } ^ { 1 } \hat { S } _ { x } ^ { 3 } + \hat { S } _ { y } ^ { 1 } \hat { S } _ { y } ^ { 3 } ) + J _ { 2 } ( \hat { S } _ { x } ^ { 2 } \hat { S } _ { x } ^ { 3 } + \hat { S } _ { y } ^ { 2 } \hat { S } _ { y } ^ { 3 } ).$$

We assume that the third spin interacts only weakly with both spins at sites 1 and 2, namely, | J 1 | , | J 2 | /lessmuch | J 0 | . By employing second-order perturbation theory with the dominant Hamiltonian ˆ H 0 = J 0 ( ˆ S S 1 x ˆ 2 x + ˆ ˆ ) S S 1 y 2 y and the perturbing Hamiltonian ˆ V = J 1 ( ˆ S S 1 x ˆ 3 x + ˆ ˆ ) + S S 1 y 3 y J 2 ( ˆ S S 2 x ˆ 3 x + ˆ S S 2 y ˆ ) 3 y , we can solve for the coherence dynamics of each spin under the spin-echo sequence. Without loss of generality, we focus on the coherence dynamics of spin 1, P 1 , whose analytical expression is shown to be the sum of oscillations at different frequencies with a DC offset:

$$P _ { 1 } = P _ { 1 } ^ { \text{Dc} } + P$$

where the individual terms on the right-hand side are defined as follows:

$$P _ { 1 } ^ { \text{DC} } = \frac { J _ { 2 }$$

$$P _ { 1 } ^ { J _ { 0 } + \nu _ { 0 } } & =$$

$$\begin{array} {$$

$$& \ell _ { 1 } & 2 J _ { 0 } ^ { 2 } \stackrel {$$

$$P _ { 1 } ^ { \frac { J _ { 0 } } { 2 } + \nu _ { 3 } } & = \frac { 1 } { 2 } \left ( \frac { J _ { 1 } - J _ { 2 } } { 2 J _ { 0 } } + \frac { 3 ( J _ { 2 } ^ { 2 } - J _ { 1 } ^ { 2 } ) } { 4 J _ { 0 } ^ { 2 } } \right ) \cos \left ( ( \frac { J _ { 0 } } { 2 } + \nu _ { 3 } ) \tau \right ) \\ \omega _ { J _ { 0 } / 2 + \nu _ { 4 } } & \quad 1 / J _ { 2 } - J _ { 1 } \quad J _ { 1 } ^ { 2 } - J _ { 2 } ^ { 2 } \ \quad ( J _ { 0 } \quad \cdots \ \quad \, \quad \, \dots \,.$$

$$P _ { 1 } ^ { \frac { J _ { 0 } } { 2 } + \nu _ { 2 } } & = \frac { 1 } { 2 } \left ( 1 + \frac { J _ { 1 } + J _ { 2 } } { 2 J _ { 0 } } - \frac { J _ { 1 } ^ { 2 } + 3 J _ { 2 } ^ { 2 } + 4 J _ { 1 } J _ { 2 } } { 4 J _ { 0 } ^ { 2 } } \right ) \cos \left ( ( \frac { J _ { 0 } } { 2 } + \nu _ { 2 } ) \tau \right ) \\ \sim \frac { J _ { 0 } } { 2 } + \nu _ { 3 } & = 1 / J _ { 1 } - J _ { 2 } \quad 3 ( J _ { 2 } ^ { 2 } - J _ { 1 } ^ { 2 } ) \, \bigwedge \quad ( \, J _ { 0 } \quad \cdots \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \$$

$$P _ { 1 } ^ { J _ { 0 } / 2 + \nu _ { 4 } } & = \frac { 1 } { 2 } \left ( \frac { J _ { 2 } - J _ { 1 } } { 2 J _ { 0 } } + \frac { J _ { 1 } ^ { 2 } - J _ { 2 } ^ { 2 } } { 4 J _ { 0 } ^ { 2 } } \right ) \cos \left ( ( \frac { J _ { 0 } } { 2 } + \nu _ { 4 } ) \tau \right )$$

Figure S10: Comparison between numerical simulation and perturbation theory for three resonant spins. We monitor the short-time ( a, b ) and long-time ( c, d ) coherence dynamics of a single spin interacting strongly with its nearest-neighbor spin with strength J 0 and with a farther spin with strength J 1 . The second and third spins interact with strength J 2 . All spins are initially polarized along the y -axis and subsequently subject to resonant spin-exchange Hamiltonian dynamics (Eq. S26). In a,c , J /J 1 0 = J /J 2 0 = 0 2 . , while in b,d , J /J 1 0 = 0 2 . , J 2 /J 0 = -0 3 . . In all cases, J 0 = 2 π × 0 3 . MHz is fixed. For the short-time dynamics ( a, b ), we compare the full numerical simulations (blue) with predictions by perturbation theory (orange). For the long-time dynamics ( c, d ), we show the full numerical simulations (blue) alongside the offset oscillation of the slow-frequency component predicted by perturbation theory (orange, Eqs. S28 and S29).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-0e7f25f0-fff6-4ded-b7b4-1b04c195c26f.jpg)

As shown in Fig. S10, the analytical perturbative solutions are validated by comparing them with numerical simulations, showing reasonable agreement. Interestingly, similar to the disordered two-spin case, the resonant three-spin scenario also exhibits an offset (Eq. S28) in coherence oscillation with reduced contrast, as the additional third spin prevents the spin-exchange dynamics between spins 1 and 2 from perfectly 'rephasing' back to their original condition (Fig. S10a,b). Moreover, we identify a very slow oscillation component with frequency J J 1 2 J 0 (Eq. S29) in the three-spin model (Fig. S10c,d), which results in the slow late-time relaxation upon ensemble averaging. The offset (Eq. S28), as well as the slow frequency component (Eq. S29) appearing in the large-spin model, are attributed to the long tail observed in the latetime decoherence profile in Model II (see Fig. 3 of the main text).

Finally, the early-time coherence decay rate of spin 1 is primarily dominated by J 0 , the strongest interaction strength between itself and spin 2, with a small correction due to its weak coupling to the third spin, J 1 :

$$\frac { d P _ { 1 } } { d \tau } \approx - \frac { J _ { 0 } ^ { 2 } + J _ { 1 } ^ { 2 } } { 4 } \tau \approx - \frac { J _ { 0 } ^ { 2 } } { 4 } \tau,$$

when τ /lessmuch 1 / J | 0 | and J 2 0 /greatermuch J 2 1 (Fig. S11).

Figure S11: Short-time behavior of the three-spin model for different interaction strengths. Three different conditions of ( J /J 1 0 , J /J 2 0 ) are considered for numerical simulation, as indicated in the figure legend of the left panel. The blue trace with J /J 1 0 = J /J 2 0 = 0 represents the resonant two-spin case. In all cases, J 0 = 2 π × 0 3 . MHzis fixed. The right panel shows that the three-spin cases (orange and green) exhibit reduced contrast with beating frequencies, due to the hindered perfect rephasing of coherence dynamics caused by the additional spin. Refer to Sec. 4.2 for discussion.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-247151b8-cad1-4859-a232-5f35befeb1b5.jpg)

## 5 Additional experimental data

In this section, we provide additional experimental data and corresponding analyses related to the experiments shown in Figures 4 and 5 of the main text.

## 5.1 WAHUHA-echo measurement

In the WAHUHA-echo sequence [56], the base pulse sequence consists of π/ 2 and π pulses with judiciously chosen rotation axes, which are repeatedly applied to the spin system. The base sequence has a periodicity of 6 τ , where τ here represents the time separation between the centers of two adjacent π/ 2 pulses (Fig. S12a). To extend the interrogation up to the total time T , the base sequence with a fixed τ is repeated k times, and coherence is measured stroboscopically at discrete times at the end of each repetition. Specifically, we sweep the rotation axis θ of the π/ 2 analyzer pulse at the end of the entire sequence, fit the resulting photoluminescence signal as a function of θ , and extract the contrast (see Fig. S4b for details). To quantify the decoherence rates under the WAHUHA-echo sequence, we fit the coherence decay profile as a function of T with a single exponential decay (blue line, Fig. S12b) and extract the 1 /e decay times (Fig. S12c). The 1 /e decay times under the WAHUHA-echo sequence are characterized as a function of the π/ 2 pulse separation τ for both the small J and large J cases (orange/blue markers, Fig. S12c).

From the theoretical perspective of Floquet Hamiltonian engineering, a shorter τ more ef-

a

c

Figure S12: WAHUHA-echosequence measurement. a, WAHUHA-echo sequence. The base pulse sequence consists of π/ 2 and π pulses with a periodicity of 6 τ , where τ denotes the time separation between the centers of two adjacent π/ 2 pulses. We maintain a spacing of τ -t p between each pulse, where t p is the duration of the π/ 2 pulse. Note that the duration of a π pulse is 2 t p . The base sequence is repeated k times to evolve over a total interrogation time T . The final π/ 2 analyzer pulse with a variable phase angle θ is employed for coherence extraction. b, Comparison of the WAHUHA-echo coherence dynamics between numerical simulation (red) and experiment (blue). The experiment is carried out under the small J condition, with error bars representing the standard deviation of the data (see Sec. 1.3 for details). The simulation is conducted using the experimentally calibrated interaction and disorder strengths. The blue line is a single exponential fit to extract the 1 /e decay time. c, The 1 /e decay times of the WAHUHA-echo coherence measurements are shown as a function of the pulse separation τ for both the large J (blue) and small J (orange) cases.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-b09d22b0-ca9c-401f-93de-c865837b5467.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-adede315-6c82-42a4-8718-c0117be64e0c.jpg)

fectively suppresses the contribution from higher m th -order Hamiltonians, thereby leading to the effective Hamiltonian being dominated by the zeroth-order Hamiltonian, ˆ H eff ≈ ˆ H 0 , in the resultant spin dynamics (Eq. S20). Since the zeroth-order Hamiltonian of the WAHUHA-echo sequence is the isotropic Heisenberg interaction Hamiltonian, i.e., ˆ H 0 = 2 3 ˆ H Heis , this enables us to protect the coherence of the polarized spin ensemble by effectively decoupling spin-spin interactions.

We experimentally confirm that coherence is indeed better preserved as τ decreases, indicating the dominance of the zeroth-order average Hamiltonian in the Floquet-engineered manybody dynamics (Fig. S12c). However, in experiments, control imperfections due to finite pulse

duration and rotation angle errors accumulate more for smaller τ because a larger number of pulses are needed to achieve a longer interrogation time T . We speculate that these imperfections led to the discrepancy between the numerical simulation (red line, Fig. S12b) and the experiment (blue markers, Fig. S12b), as well as the saturation of coherence improvement as τ approaches the finite pulse duration, i.e., τ ≈ 20 ns, for the large J case (Fig. S12c). In our experiments, the maximum coherence times are attained when τ = 22 ns and τ = 33 ns for the small J and large J cases, respectively. We anticipate that coherence times could improve with the implementation of shorter pulse durations and more robust pulse sequences [56].

## 5.2 Spin-locking measurement

As presented in Figure 4b of the main text, the spin-locking sequence is employed to investigate how long the polarized spin ensemble remains stably locked along the y -axis under a continuous drive with strength Ω y . Specifically, under the spin-lock drive, the zeroth-order average Hamiltonian can be shown to be

$$\hat { H } _ { 0 } = \Omega _ { y } \sum _ { i } ^ { N } \hat { S } _ { y } ^ { i } + \frac { 1 } { 2 } \sum _ { i j, \, i > j } ^ { N } J _ { i j } ( \hat { S } _ { y } ^ { i } \hat { S } _ { y } ^ { j } + \vec { S } ^ { i } \cdot \vec { S } ^ { j } ) = \hat { H } _ { \text{on-site} } ^ { y } + \frac { 1 } { 2 } \hat { H } _ { \text{In-sing} } ^ { y } + \frac { 1 } { 2 } \hat { H } _ { \text{Heis} }. \quad \ \quad \ \text{(3.36)} \\ \hat { \ } \ \hat { \ } \ \hat { \ } \. \quad \ \text{$\hat{ }$}.$$

Note that the zeroth-order Hamiltonian can be approximated as ˆ H 0 ≈ ˆ H y on-site = Ω y ∑ N i ˆ S i y when Ω y /greatermuch J ij . While in the leading order, the polarized spin ensemble can be protected against interaction-induced dephasing since the initial state is an eigenstate of ˆ H 0 , higher-order average Hamiltonians, as well as environmental noise from the spin bath creating AC noise at frequency Ω y , could induce dephasing of the spin-locked signal.

In Figure S13, we experimentally probe the decoherence dynamics of the spin-locking sequence for both small J and large J cases. The spin-lock Rabi frequency is fixed at Ω y ≈ 2 π × 10 MHz, significantly greater than other rates in the system, including large J = 2 π × 0 35 . MHz, small J = 2 π × 0 19 . MHz, and W = 2 π × 0 65 . MHz. The decay time traces are fitted with a single exponential function, yielding 1 /e decay times of ≈ 50 µ s and ≈ 73 µ s for the large J and small J cases, respectively. We attribute the different decay times observed at different J values to contributions from the aforementioned higher-order average Hamiltonians and varying levels of device heating induced by optical power during spin density tuning.

Figure S13: Spin-locking sequence measurement. a, Spin-locking sequence. After initializing the spin ensemble along the y -axis via the initial π/ 2 pulse, a continuous drive along the y -axis with strength Ω y is applied for duration T to lock the spin orientation. The final π/ 2 analyzer pulse with a variable phase angle θ is employed for coherence extraction. b, Decoherence dynamics of the spin ensemble under the spin-locking sequence for large J (blue) and small J (orange) cases. Error bars represent the standard deviation of the experimental data (see Sec. 1.3 for details). The solid lines are single exponential fits used to extract the 1 /e decay times, T 1 /e , yielding T 1 /e ≈ 50 µ s and T 1 /e ≈ 73 µ s for the large J and small J cases, respectively.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-75604ea8-9ad5-4a54-a60b-02bd823d1432.jpg)

## 5.3 Robustness of DTC phases to initial states

In Figure S14, we investigate the robustness of discrete time-crystalline (DTC) phases relative to initial spin states. The Floquet pulse sequence used to investigate the robustness against initial states is similar to Fig. 5a of the main text. Here, however, we apply a variable angle φ ranging from φ = 0 to φ = π/ 2 to the initialization pulse, preparing initial spin orientations along the z -axis and y -axis, respectively (Fig. S14a). For each initial state specified by a rotation angle φ , we sweep the perturbation strength, /epsilon1 , and the interaction time, τ , of the base control sequence to generate a DTC phase diagram (Fig. S14b-f). Details on the construction procedures for DTC phase diagrams can be found in the main text.

We experimentally observe that when the initial spin orientation is near the y -axis, for instance, when | φ -π/ 2 | /lessorsimilar π/ 4 (Fig. S14d-f), the DTC phases robustly exhibit characteristic linear phase boundaries (Fig. S14g). We observe that the slope of the phase boundaries remains consistent for different φ values, implying the robustness of the observed subharmonic DTC

a

Figure S14: Robustness of DTC phases to initial spin states a, Floquet control pulse sequence for probing DTC phases with different initial spin orientations, controlled by the first pulse with a variable rotation angle φ ; otherwise, the sequence is the same as in Fig. 5a of the main text. b-f, DTC phase diagrams for different initial rotation angles ranging from φ = 0 to φ = π/ 2 . See the main text for details on the phase diagram reconstruction. g, Phase boundaries for d-f where the initial spin orientations are near the y -axis. The phase boundaries for each interaction time τ are determined by identifying a threshold in the value of /epsilon1 where the subharmonic peak intensity at frequency ν = 1 2 / , | S ν ( = 1 2) / | 2 , falls below 0.4, i.e., | S ν ( = 1 2) / | 2 &lt; 0 4 . . We note a lateral 1% offset in the swept angle /epsilon1 at the center position of the DTC phase diagram, i.e., /epsilon1/π ≈ 0 01 . , attributed to an experimental calibration error of the rotation pulse.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00252v1-figure-140cd912-9fb1-46e3-8b96-dda264b3265a.jpg)

signals relative to the initial spin states.

In contrast, when the initial spin orientation is near the z -axis (Fig. S14a,b), the subharmonic DTC signals are visible only within a small region near τ ≈ 0 and /epsilon1 ≈ 0 , indicating rapid decoherence under the applied Floquet sequence. We attribute this to the degradation of the underlying spin-locking performance, as the spin ensemble can be dephased by the dominant on-site field along the y -axis when the spin orientation deviates too much from the spinlocking field direction. In other words, when the initial spins are aligned close to the z -axis, the continuous drive induces Rabi precession instead of spin locking.

## References

- [51] Lei, M. et al. Many-body cavity quantum electrodynamics with driven inhomogeneous emitters. Nature 617 , 271-276 (2023).
- [52] Kindem, J. M. et al. Characterization of 171 Yb 3+ :YVO 4 for photonic quantum technologies. Phys. Rev. B 98 , 024404 (2018).
- [53] Davis, E. J. et al. Probing many-body dynamics in a two-dimensional dipolar spin ensemble. Nature Physics 19 , 836-844 (2023).
- [54] Bartholomew, J. G. et al. On-chip coherent microwave-to-optical transduction mediated by ytterbium in YVO . 4 Nature Communications 11 , 3266 (2020).
- [55] Ruskuc, A., Wu, C.-J., Rochman, J., Choi, J. &amp; Faraon, A. Nuclear spin-wave quantum register for a solid-state qubit. Nature 602 , 408-413 (2022).
- [56] Choi, J. et al. Robust dynamic hamiltonian engineering of many-body spin systems. Physical Review X 10 , 031002 (2020).
Mi Lei Rikuto Fukumori Chun-Ju Wu Edwin Barnes Sophia Economou Joonhee Choi Andrei Faraon
Quantum thermalization and Floquet engineering in a spin ensemble with a clock transition
Studying and controlling quantum many-body interactions is fundamentally
important for quantum science and related emerging technologies. Optically
addressable solid-state spins offer a promising platform for exploring various
quantum many-body phenomena due to their scalability to a large Hilbert space.
However, it is often challenging to probe many-body dynamics in solid-state
spin systems due to large on-site disorder and undesired coupling to the
environment. Here, we investigate an optically addressable solid-state spin
system comprising a strongly interacting ensemble of millions of ytterbium-171
ions in a crystal. Notably, this platform features a clock transition that
gives rise to pure long-range spin-exchange interactions, termed the dipolar XY
model. Leveraging this unique feature, we investigate quantum thermalization by
varying the relative ratio of interaction strength to disorder, dynamically
engineering the XY model into other many-body Hamiltonian models, and realizing
a time-crystalline phase of matter through periodic driving. Our findings
indicate that an ensemble of rare-earth ions serves as a versatile testbed for
many-body physics and offers valuable insights for advancing quantum
technologies.