Signatures of discrete time crystalline order in dissipative spin ensembles

Discrete time-translational symmetry in a periodically driven many-body system can be spontaneously broken to form a discrete time crystal, an exotic new phase of matter. We present observations characteristic of discrete time crystalline order in a driven system of paramagnetic P-donor impurities in isotopically enriched 28Si cooled below 10 K. The observations exhibit a stable subharmonic peak at half the drive frequency which remains pinned even in the presence of pulse error, a signature of discrete time crystalline order. This signal has a finite lifetime of ∼100 Floquet periods, but this effect is long-lived relative to coherent spin–spin interaction timescales, lasting ∼104 times longer. We present simulations of the system based on the paradigmatic central spin model and show good agreement with experiment. We investigate the role of dissipation and interactions within this model, and show that both are capable of giving rise to discrete time crystal-like behaviour.


Introduction
A central paradigm in condensed matter physics has been to classify phases of matter by their symmetries. Indeed, spontaneous symmetry breaking describes many known phase transitions. A common symmetry-breaking phase is a crystal in real space, where the symmetry under continuous spatial translation is broken to a lower discrete one. A natural question is then to ask whether it is possible to analogously break time-translation symmetry [1]. Breaking of continuous time-translation symmetry has been shown to be impossible in thermal equilibrium [2,3], but periodically-driven ('Floquet') systems provide the means to break discrete time-translation symmetry, thereby forming a discrete time crystal [4][5][6][7][8][9][10]. Observations consistent with discrete time-crystalline behaviour were reported soon after in experiments [11][12][13][14].
In this article we report the observation of signatures of a discrete time crystal (DTC) in silicon doped with phosphorus. Silicon provides an ideal platform for implementing the dynamic pulse sequences crucial for realizing time crystals, owing to its ability to be isotopically engineered to an exceptionally high purity [15], having the longest coherence times of any solid state system [16,17], and its versatility with dopant species and concentration.
In phosphorus-doped silicon, the dopant electron spins interact via dipolar interactions. At donor concentrations below about 10 16 cm −3 these interactions are weak compared to dissipative effcts due to magnetic impurities and charge noise. Dissipation and nonlinearities provide a natural route to produce robust period-doubling, as has been observed with AC-driven charge density waves [18] and Faraday wave instabilities [19]. For intuition, one could imagine a dissipative system with two basins of attraction, and a periodic drive which flips between these two attractors. It is clear when the flip is perfect that, starting from one of the attractors, this system will exhibit period-doubling. Even when the flip is not perfect, provided the state is sufficiently close to one of the attractors, the dissipation will cancel out any imperfections and yield robust breaking of time-translation symmetry [20][21][22].
We observe signatures of robust DTC order over 200 Floquet periods using a sample of phosphorus-doped silicon in the strong dissipation and weak interactions regime, indicating that this order is indeed stable to weak interactions. We go on to produce a theoretical phase diagram of this dissipation-stabilized order as a function of dissipation rate and interaction strength, and experimentally probe the DTC and trivial phases. We produced the phase diagram using the driven central-spin model [17,23,24] coupled to a dissipative bath, direct simulations of which using experimental parameters agree well with our experimental observations. As a point of independent interest, in section 3.1 we comment on the use of the crystalline fraction as a DTC order parameter. In the experiment we used composite BB1 pulses [25,26] to enable extremely uniform spin rotations. This is important because in experiments DTC order is usually observed by fixing some nonzero rotation error and probing how robust the DTC order is to this error. One proposed experimental probe is the crystalline fraction [11,13,27], which roughly measures the strength of the subharmonic Fourier peak compared with the rest of the spectrum. A 'plateau' in the crystalline fraction across a finite window of rotation errors has been used as an indicator of DTC order. However, we argue that this plateau can arise simply from non-uniform spin rotations, since the plateau disappears when we use BB1 pulses to ensure uniform rotations. This indicates that using the crystalline fraction as a witness for DTC order can obfuscate the different effects.
We note here that we do not expect this system to be many-body localized (MBL) [28,29]. MBL has been argued to be a necessary condition for the existence of time-crystallinity [30], but recent observations of a discrete time crystal have been reported in systems which may not be many-body localized [11][12][13], due to the presence of a bath, and dipolar interactions [31,32], which are long-ranged in 3D. Though it has been suggested that MBL can survive in systems with long-range interactions [33], we do not expect those arguments to hold here, as they rely upon the Anderson-Higgs mechanism.

Materials and methods
A periodic pulse sequence is implemented on the electronic spins as shown in figure 1(a). We use a sample of 28 Si enriched to 99.995%; this provides a magnetically clean environment which gives the dopant spins an exceptionally narrow linewidth of less than 10 μT. The sample is placed in a Bruker X-band microwave resonator and microwave driving is applied at a frequency of 9.775 GHz. Full details of the apparatus used are given in section 6.1 (https://stacks.iop.org/NJP/22/085001/mmedia). A global magnetic field B 0 is applied along the z direction in the vicinity of the sample, about which the spins within the ensemble undergo Larmor precession. We use a Bloch sphere representation with the ground state defined in the +z direction and consider all spin dynamics in the rotating frame, which rotates with the spins at the Larmor frequency ω about the z axis. Microwave pulses with frequency ω appear in the rotating frame as a static magnetic field B 1 in the x-y plane. The phase of the pulse determines its orientation within the x-y plane. The spins precess about this B 1 field in the rotating frame and thus we can use this to drive the spins around the Bloch sphere. We begin with a spin ensemble polarised in the +z direction, parallel to the applied magnetic field with a setpoint of 3512 G. The system is initialised with a π/2-pulse, rotating the spins into the x-y plane. Once rotated away from the ground state, the spins immediately begin to undergo phase decoherence, which causes the spin ensemble magnetisation to decay towards the centre of the Bloch sphere. Relaxation along the z axis towards the ground state also occurs but this process is much slower. While interactions between spins are always present, at long times the spins all relax towards +z and the dynamics are uninteresting. Following the initial π/2-pulse, a long microwave 'spin-locking' pulse is applied along the x direction for a duration τ lock , which serves to allow the spins to interact while decoupling them from their environment. We then immediately rotate about the y axis by π with a deliberate error . This rotation is achieved via the use of a composite BB1 pulse [25,26], which compensates for variations in the Rabi frequency across the sample. This was crucial to removing artefacts in the data that can otherwise be mistaken for features of a DTC (see section 3.1). Together, the spin-locking and BB1 pulse form a single Floquet unitary. After applying N Floquet unitaries we rotate the ensemble back to the z axis with a π/2-pulse. We then read out the resulting polarisation using a π/2-π Hahn echo sequence. The time τ lock is chosen to correspond to an integer multiple of 2π-rotations of the spins about the x axis to avoid any dynamical decoupling effects. For larger values of , we return to the trivial symmetry-unbroken phase. (d) Comparison of the Fourier peak location as a function of rotation error , using a double Lorentzian fit, for both the DTC and trivial phases, corresponding to the white and red stars in section 3.2 respectively. In the DTC phase, there is a finite window around = 0 where the peak remains rigidly locked at half the drive frequency, indicating the robustness of the DTC phase. In the trivial phase, the peak splits for any > 0, providing a clear contrast to the DTC phase. The asymmetry around = 0 is because for this data, corresponding to the τ lock = 10π data in figure 2(c), there was some variation in the π-rotation duration, which meant that data was only gathered in the region −0.02 0.08. We attribute this to small variations in the pulse amplifier power output depending on the duty cycle. Parameters (c, d):

Results
We sweep the the number of Floquet cycles n from 0 to 200 with error and map out the integrated echo amplitude as a function of n. The result is a decaying oscillating signal with frequency ν = 1/2. We take the Fourier transform of these oscillations and examine what happens to the ν = 1/2 peak as we change . In the absence of interactions between electron spins, when we increase on the π-pulse in each Floquet cycle we expect the cumulative error caused by successive over-or under-rotation of the spin ensemble to cause modulation of the oscillations. This is apparent in the Fourier transform as a splitting of the ν = 1/2 peak, as shown in figure 1(c). Increasing , we expect to see these peaks split further apart as the error in the π rotations accumulates faster leading to faster modulation. However, in the presence of sufficiently strong dissipation, we do not see an immediate splitting of the ν = 1/2 peak, but instead find a region where the oscillations remain resilient to error.
We use a sample of 28 Si doped with 1 × 10 15 cm −3 phosphorus. The phase coherence time T 2 was measured as 260 μs, spin-lattice relaxation time T 1 was measured as 1097 μs and the Rabi frequency was 1.650 MHz, corresponding to a π-pulse duration of 303 ns. The decay rate of the driven oscillations T 1ρ was measured at 190 μs from an exponential fit.
We can tune the spin-locking time τ lock . Increasing τ lock increases both the dissipation and interaction time per Floquet cycle of the spins. For sufficiently long τ lock we observe a range of values of for which the central ν = 1/2 peak does not split. Further increasing τ lock increases the range of for which we have a Note that this diagram describes the dynamics in the rotating frame. (c) Experimental data of the full width at half maximum of the Fourier peak of the net magnetization as a function of rotation duration, plotted for different interaction times τ lock . The lines are quadratic fits; solid line region indicates the points that were used for data fitting, beyond which the data began to deviate from a quadratic dependence. The fit indicates that for small rotation errors the peak width follows the expected 2 dependence for a dissipative DTC. The π-rotation duration varied slightly for different spin-locking times, which we attribute to variations in the pulse amplifier power output depending on the duty cycle. This results in slight variation in the position of the minima of the quadratic fits. (inset) Crystalline fraction against rotation duration. This is defined as the ratio of the ν = 1/2 peak to the total spectral power, f = |S(ν = 1/2)| 2 / ν |S(ν)| 2 [11], and is a measure of the total fraction of spins in the DTC phase. The fraction decreases as the rotation duration deviates further from π as expected; the rate of decrease is slower for longer spin-locking. single peak, indicating that the oscillations have increased resilience to errors. A clear comparison between the phenomenology of the DTC and the trivial phases can be found in figure 1(d), which shows the value of at which the subharmonic peak splits in the two phases, corresponding to the white and red stars in section 3.2. In the trivial phase, the peak splits for any > 0, but upon entering the DTC phase the peak remains pinned at ν = 0.5ν F for small but nonzero, indicating the resilience to perturbations of the DTC phase. There is also good agreement in the critical value of at which the peak splits in the DTC phase between the experiment and the simulation using the model described in section 3.3.
One indication that the signatures of DTC order in this experimental system are effectively described by a dephasing model can be found by looking at the dependence of the Fourier peak width on the rotation error . In the strong-dephasing limit, analysis of a simple single-spin model indicates that the YZ-dephasing results in exponential decay of | σ x | at a rate Γ ∼ −log[cos π] ∼ ( π) 2 /2 per Floquet period, which is quadratic in for 1 (see section 6.2). Figure 2(c) shows the Fourier peak width as a function of for the experimental data, which shows the predicted 2 dependence. There are additional contributions to the peak width from other sources, such as errors from finite rotation times, but the 2 contribution can be seen as a hallmark of the dissipative phase.
Further, in figure 6.2 we investigate another sample doped with 3 × 10 15 cm −3 P in 28 Si and show that despite a 200% increase in spin concentration, the key features remain the same, indicating that we are still in the dephasing dominated regime.

Effects of inhomogeneous spin rotations on the crystalline fraction
The crystalline fraction, defined in terms of the Fourier transform S(ν) as |S(ν = 1/2)| 2 / ν |S(ν)| 2 , has been recently utilised as an experimental probe of Z 2 discrete time crystal robustness [11,13,27]. A 'plateau' in the crystalline fraction in a finite window around zero rotation error has been used an indicator for DTC order. However, we have observed that this plateau can emerge simply as a result of having nonuniform rotation pulses across the extent of the sample.
Experimentally, this nonuniform rotation comes in our case from variations in the microwave frequency magnetic field generated by the split ring cavity used for driving and detection. To facilitate a comparison of the effect of these nonuniform rotations, we used a composite rotation pulse known as a BB1 pulse, which allows us to correct variations across our sample, estimated to be on the order of ±10%, up to 5th order. The BB1 pulse takes the form of a simple θ rotation pulse, immediately followed by a π(φ) − 2π(3φ) − π(φ) corrective pulse sequence, where the phase of each pulse φ = arccos(−θ/4π). We used a BB1 pulse during the main experiment to ensure an extremely uniform rotation of the spin ensemble, which enables us to properly diagnose the robustness of the DTC phase to deliberate rotation errors. Figure 3(a) shows a comparison of the crystalline fraction as a function of θ rotation error with and without BB1 pulses. It should be noted that while the two experiments were essentially the same, the experimental apparatus used was different between these two runs. Without BB1, the nonuniform rotations result in a flattening of the crystalline fraction, which could be seen as a 'false-positive' for DTC order. However, if we do use BB1 pulses to ensure uniform rotations, the 'plateau' disappears, and the crystalline fraction is simply peaked around the point of zero pulse error. To corroborate these results, we performed simulations of the driven-dissipative central spin model using nonuniform rotation pulses ( figure 3(b)). For a given disorder realization, the rotation error is drawn from a Gaussian distribution with mean¯ and variance σ 2 . The maximum deviation from¯ is fixed at 10% to model the finite extent of the sample. The same rotation error is used periodically up to 200 Floquet cycles, from which we calculate the crystalline fraction. This is repeated 500 times and averaged to produce a given curve in figure 3(b). Here we observe results consistent with the experiment: larger variance σ 2 in the rotation error results in a distinct flattening of the crystalline fraction.

Phase diagram
Discrete time crystals can be characterized by their robust subharmonic response to a periodic drive. Intuitively, the dissipation can cancel out any small perturbations to the state or the drive which might otherwise break the subharmonic response (see figure 2(b)). What is more non-trivial is whether this dissipation-stabilized subharmonic response can persist in the presence of destabilizing interactions. In this section we map out the phase diagram of a discrete time crystal (DTC) with competing dissipation and interactions, and numerically show that a dissipation-driven DTC is indeed stable to weak interactions ( figure 4).
In the experiment outlined in sections 2 and 3, we explore the strong dissipation/weak interaction regime of this phase diagram, and find that our observations align with the theoretically predicted phase boundary (see section 3.2). At the other extreme, the weak dissipation/strong interaction regime could be probed by looking at samples with higher concentrations of phosphorus dopants.
In an experiment, a DTC is typically detected by measuring the time-dependence of local observables (or averages thereof), and looking for oscillations at an integer fraction 1/n of the drive frequency ν 0 which are robust against perturbations. This robustness can be detected by looking at the Fourier spectra of these local observables; the DTC phase will have a strong peak at ν 0 /n which remains unsplit for a finite window of rotation errors. To use this to produce a phase diagram, we fix a nonzero rotation error , and then for each set of parameters calculate the Fourier spectrum of σ x 0 in the driven-dissipative central spin model

Theoretical model
To produce the phase diagram, we focus on the driven central spin model (CSM) coupled to a dissipative bath (see figure 2(a)). The CSM has been successful as a semiclassical effective model for describing decoherence in solid state systems [34][35][36][37][38][39]. In the experiment described in sections 2 and 3, the strongest coherent interaction is between the electronic spins of the phosphorus dopants, which form the spins in the CSM. However, in the dephasing-dominated regime relevant to this experiment, we do not expect the spatial interaction structure to be important, and view the infinite-range central spin model as an approximation to the long-range dipolar interactions of the phosphorus spins. This theoretical approximation has been validated for describing the decoherence of donor spins in experiments [17,23,24]. Finally, nuclear spins from residual 29 Si contribute to dissipation rates, and provide a small source of random field for the electron spins.
To have a notion of competing interactions and dephasing, we use σ y σ y interactions and σ x dephasing, which in principle should produce DTC-like phenomena in σ y and σ x respectively. Interestingly, in the weakly interacting case, we find that strong σ x dephasing can also stabilize a DTC-like response in σ y .
The Hamiltonian of the central spin model with σ y σ y interactions is given by where the central spin has index 0, and the outer spins have indices 1 to N. J 0i and h j are modelled as random variables taken from the uniform distributions over [−h, h] and [−J y , J y ] respectively. To incorporate the driving, the Hamiltonian follows a two-stage protocol given by For = 0, the second stage of the pulse protocol exactly flips the spins. The discrete time crystal phase can be defined operationally by its robustness against nonzero rotation errors . Note that we have τ flip = (1 + )t π + 4t π , where t π is the time to do an uncorrected π-rotation, because the BB1 pulse is made up of the uncorrected rotation by (1 + )π, followed by a π-2π-π corrective sequence. For simplicity, we do not simulate the final rotation from the x-axis to z-axis before readout mentioned in section 2, since this is needed only for experimental purposes; we measure all our magnetizations along the x-axis for the simulations.
In addition to the interactions present within the central spin model, we model the effects of external dephasing using the Lindblad master equation given by Note that this dissipation model enacts YZ-dephasing, i.e. it draws the state of a single qubit to the x-axis of the Bloch sphere. This can be seen as being in competition with the σ y σ y interactions present within the central spin model (equation (1)).
In the experiment, T 1ρ and T 1 were measured to be 193 μs and 1097 μs respectively for the higher density sample, giving κ 1 = 1/2T 1ρ − 1/4T 1 ∼ 2.3 kHz, while J and h can be estimated from the phosphorus concentration and are of the order 300 Hz and 10 Hz respectively. Hence we expect YZ-dephasing to dominate the dynamics of the experimental system.
In producing figure 1(c), we also included T 1 -type dissipation via the Lindblad operators where κ 2 = 1/T 1 ∼ 0.9kHz using the experimental T 1 . We checked numerically that this effect does not significantly affect the phase boundary, and serves merely to slightly broaden the Fourier peaks.
In the strong dephasing limit, both the central and the outer spins exhibit robust period-doubling. The same is also true in the strong interactions limit, where interestingly the outer spins have Fourier spectra which split for larger values of than the central spin.

Comparison with experiment
The experiment described in section 3 probes the dephasing-driven DTC transition. The detail to section 3.2 shows experimental Fourier spectra for experiments with τ lock = 10t π and τ lock = 100t π . Increasing τ lock lengthens the Floquet period and hence reduces the Floquet frequency ν F . Since the axes in the phase diagrams are measured in units of ν F , this allows us to explore regions of the phase diagram with effectively stronger interactions and dephasing. These experimental Fourier spectra thus probe either side of the dissipative DTC phase boundary, and demonstrate that the dissipative DTC is stable in the predicted region, even with potentially destabilizing interactions. Furthermore, the location of the simulated phase boundary matches that observed in the experiment, as demonstrated in figure 1(d).

Discussion
In summary, we have shown experimental signatures of the formation of a discrete time-crystal phase in naturally purified silicon doped with phosphorus, using composite BB1 pulses for exceptionally uniform rotation of the spin ensemble. We observe the formation of discrete time-crystalline order by driving the electron spin ensemble at frequency ν 0 , while producing a response at a sub-harmonic frequency ν 0 /2. We show that this peak remains pinned and is robust to perturbations in the periodic pulse protocol. Motivated by the experimental system, we investigate the dissipative central spin model as a phenomenological description for time-crystalline behaviour in solid state systems with long-range interactions. We show that the model, with no free parameters, is in remarkably good agreement with experiment. Furthermore, we investigate the role of interactions and dissipation stabilizing the DTC phase in the central spin model and the crossover regime in which both effects are significant.
The high level of quantum control of electron and nuclear spins in phosphorus doped silicon provides a promising platform for studying many-body quantum coherence in driven systems. With over 10 13 spins, the system allows us to work in the true disordered many-body thermodynamic limit in which DTC was originally predicted to exist.
The system has numerous advantageous properties that could be used to further investigate out-of-equilibrium effects, such as the presence of two nearly degenerate electron spin transitions, and additional nuclear spin transitions. By driving the nuclear and electronic spins independently, proximity effects in time-crystalline behaviour and effects of dynamic nuclear polarization can be explored. The difference in the dynamic time scales of electron and nuclear spins could serve as a useful tool for manipulating DTC order and its long time coherence. Furthermore, the lifetime of the DTC can be exploited to probe the dephasing and thermalizing properties of long range systems. The DTC order exhibited in the central spin model poses several interesting questions on non-thermal steady states in Floquet systems which are ripe for further investigation [40].