${}^{31}\mathrm{P}$ NMR study of discrete time-crystalline signatures in an ordered crystal of ammonium dihydrogen phosphate

The rich dynamics and phase structure of driven systems include the recently described phenomenon of the “discrete time crystal” (DTC), a robust phase which spontaneously breaks the discrete time translation symmetry of its driving Hamiltonian. Experiments in trapped ions and diamond nitrogen vacancy centers have recently shown evidence for this DTC order. Here, we show nuclear magnetic resonance (NMR) data of DTC behavior in a third, strikingly different, system: a highly ordered spatial crystal in three dimensions. We devise a DTC echo experiment to probe the coherence of the driven system. We examine potential decay mechanisms for the DTC oscillations, and demonstrate the important effect of the internal Hamiltonian during nonzero duration pulses.

Figure 1

(a) Atoms in the unit cell of ammonium dihydrogen phosphate (ADP), which has chemical formula ${\mathrm{NH}}_4{\mathrm{H}}_2{\mathrm{PO}}_4$. ADP is an ionic, tetragonal crystal with space group $I\overline{4}2d$. At room temperature, the ${\mathrm{NH}}_4$ groups experience rapid in-place rotation, such that the time-averaged location of the four ${}^1\mathrm{H}$ is at the nitrogen site. The ${}^1\mathrm{H}$ in ${\mathrm{NH}}_4$ are shown in a distributed manner to reflect this. We also place the remaining so-called “acid” protons (${}^1\mathrm{H}$) in their time-averaged positions, between the lattice sites of the nearest oxygens [32, 33]. (b) The ADP crystal sample studied here, shown in a 5-mm-diameter NMR sample tube, held in place by rolled teflon tape (white).

Figure 2

(a) Magnetization decay from a Hahn echo experiment with ${}^1\mathrm{H}$ off (circles), where each data point is acquired with a Hahn echo sequence for a different value of $\tau$. We compare this to the simulated decay from an Ising-type approximation, both before (dashed line) and after (solid line) scaling $B_{ij}$ by $\frac{3}{2}$ to approximate the actual dipolar Hamiltonian. (b) ${}^{31}\mathrm{P}$ spectra as acquired by an FID with ${}^1\mathrm{H}$ on (blue squares), and by a Hahn echo with ${}^1\mathrm{H}$ off [red circles, FT of Hahn echo data in (a)], with the results of a numerical model at a single-crystal orientation (lines). (c) Comparison of the ${}^{31}\mathrm{P}$ spectrum from an FID (closed circles) to the line shape from an altered Hahn echo (open circles). The Hahn echo spectrum has been broadened using a Gaussian with FWHM 280 Hz, to account for the ${}^{31}\mathrm{P}-{}^{14}\mathrm{N}$ coupling.

Figure 3

DTC pulse sequence. After ${}^1\mathrm{H}$ spins are rotated with an $X_{\pi /2}$ pulse (tall blue block), ${}^{31}\mathrm{P}$ magnetization is created along $\widehat{y}$ via cross polarization with the ${}^1\mathrm{H}$ spins, and is then rotated into $\widehat{z}$ with an ${\overline{X}}_{\pi /2}$ pulse (tall orange block) to prepare the initial state of the system. We then apply repeated Floquet cycles consisting of a delay $\tau$ followed by a pulse $X_{\theta }$ (wide blue block). After $N$ cycles, an $X_{\pi /2}$ pulse is applied, and the magnetization is immediately measured, producing a single data point in $S\left(t\right)$. We increase $N$ by 1 and repeat the sequence, following a 3-s recycle delay. This sequence is applied for $N=1,2,...,128$. After cross polarization, continuous rf decoupling (red) can be applied (${}^1\mathrm{H}$ off) to remove the effect of the ${}^1\mathrm{H}$ (a), or decoupling can be omitted (${}^1\mathrm{H}$ on), allowing the ${}^1\mathrm{H}$ to act on the ${}^{31}\mathrm{P}$ spins (b).

Figure 4

(a) Applying repeated $\pi$ pulses with ${}^1\mathrm{H}$ decoupling at $\theta \approx \pi$ and small $\tau$, we see oscillations in the time-domain signal, corresponding to a single peak in the FT signal at $\tilde{\nu }=\frac{1}{2}$. Each data point is acquired in a separate experiment, using the DTC sequence for a given $N$. $N$ odd are in green (starting negative) and $N$ even are in blue (starting positive), with black lines between them to guide the eye. (b) Decreasing $\theta$, we observe beating in the time-domain signal, corresponding to a splitting of the Fourier peak. (c) Near the same $\theta \approx 0.962\pi$ but for increased $\tau$, the oscillations are restored, once again producing a single peak in the Fourier spectrum. (d) Significantly increasing $\tau$, we still see the same behavior. (e)–(h) We observe qualitatively similar behavior in the absence of ${}^1\mathrm{H}$ decoupling. Note that for (h), $NT$ becomes comparable to the ${}^{31}\mathrm{P}$ lattice relaxation time, $T_1^{\text{P}}=103\mathrm{s}$. In each case (a)–(h), $T=\tau +t_p$ with $t_p=7.5\mu \mathrm{s}$. Data in (a)–(d) were acquired with $\times$ the number of scans as (e)–(h), doubling the maximum value of $S\left(t\right)$.

Figure 5

Waterfall plots showing the spectra $|S\left(\tilde{\nu }\right){|}^2$ at different $\theta$ with ${}^1\mathrm{H}$ off. Near $\theta =\pi$, a prominent subharmonic response is observed at $\tilde{\nu }=\frac{1}{2}$. (a) At short drive periods $T$, the subharmonic response splits into two peaks, which diverge almost immediately as $\theta$ is adjusted away from $\pi$. (b) For long drive periods $T$, the subharmonic peak lowers in amplitude as $\theta$ is adjusted away from $\pi$, but remains rigidly locked at $\tilde{\nu }=\frac{1}{2}$.

Figure 6

${}^1\mathrm{H}$ off: crystalline fractions $f$ across a broad range of drive periods $T$ (labeled, where the $f$ data for the $m\mathrm{th}$ value of $T$ are vertically offset by $m-1$ for clarity). The crystalline fractions are well fit by Gaussians. Over the duration of the many experiments, the tuning of the NMR tank circuit can drift, leading to poorly calibrated $\theta$ (left). The black squares represent well-calibrated benchmarks which we use to correct the data to match the actual $\theta$ values (as described in the main text), resulting in the data on the right. Because of heating, the $100-\mathrm{ms}$ data will not appear in Fig. 8. Error bars (not shown) are much smaller than the markers.

Figure 7

${}^1\mathrm{H}$ on: in the absence of ${}^1\mathrm{H}$ decoupling, we are free to explore an even greater expanse in $T$ (labeled, where the $f$ data for the $m\mathrm{th}$ value of $T$ are vertically offset by $m-1$ for clarity) without the danger of the circuit heating effects seen in Fig. 6. As with ${}^1\mathrm{H}$ off, we observe Gaussian shapes in the crystalline fraction with ${}^1\mathrm{H}$ on, and the width of the Gaussians increases with the drive period. The black squares are used in the same correction procedure as those shown in Fig. 6, to correct for miscalibrations of the actual $\theta$ from the expected $\theta$ (a very minor effect here). Error bars (not shown) are much smaller than the markers.

Figure 8

(a)–(d) Probing $f$ with ${}^1\mathrm{H}$ off. (a) We establish a cutoff (dotted line) in the Gaussian fits to the crystalline fraction at $f=0.1$. Crystalline fractions from $T=20\mu \mathrm{s}$ (triangles) and $400\mu \mathrm{s}$ (circles) are shown. The intersection of $f$ with the cutoff defines a boundary point. (b) Cutoff at $f=0.1$ (red circles), corresponding to the boundaries within which we observe persistent oscillations at $\tilde{\nu }=\frac{1}{2}$ (the “DTC region”). We also show cutoffs at $f=0.05$ and 0.15 (dotted lines). We compare this to an effective dipolar interaction angle by plotting $|\theta -\pi |=W\tau$, with $W=W^{\text{P,P}}$ [(b), (c), black dashed lines] and $W=W^{\text{P,PN}}$ [(b), (c), gray dashed lines, very close to $W^{\text{P,P}}]$. (c) DTC region on a semilog scale. For $\tau =100\mathrm{ms}$, the results in Fig. 6 become unreliable because of tank circuit heating from rf decoupling, so they are not plotted here. (d) $f$ versus $\tau$ for ${}^1\mathrm{H}$ off at $\theta =1.067\pi$ [angle marked in (c)]. (e)–(h) Probing $f$ with ${}^1\mathrm{H}$ on. In (f) and (g) we also include $|\theta -\pi |=W^{\text{P,HPN}}\tau$ (blue dotted-dashed lines). In (g), the data span the range $0.03<W^{\text{P,P}}\tau <3200$ radians. Error bars (not shown) are much smaller than the markers in (a)–(h).

Figure 9

Using Fourier transforms of only 50 late-time points in $S\left(NT\right), N=51$–100, the crystalline fractions become flatter around $\theta =\pi$, for both ${}^1\mathrm{H}$ off (left, red triangles) and ${}^1\mathrm{H}$ on (right, blue circles). We fit these to symmetrical super-Gaussians (lines): $F\left(\theta \right)=Aexp[-(|\theta -{\theta }_0{|/\sigma )}^p/2]$, where we fix ${\theta }_0$ using the Gaussians in Figs. 6 and 7.

Figure 10

(a) DTC echo sequence, designed to approximately reverse the effect of the original DTC sequence. The “approximate reversal” block consists of a rotation ${\overline{X}}_{\theta }$ (wide orange block), followed by a duration $2\tau$ during which a strong pulse of phase $y$ is applied to the ${}^{31}\mathrm{P}$. We apply “wrapper” pulses $X_{\pi /2}$ and ${\overline{X}}_{\pi /2}$ (tall blue and orange blocks, respectively) to rotate $-{\mathcal{H}}_{yy}^{\text{P,P}}$ into $-{\mathcal{H}}_{zz}^{\text{P,P}}$. Since the last two pulses of the sequence negate one another, neither is applied in practice. ${}^1\mathrm{H}$ decoupling is used throughout. (b), (c) DTC echoes for $T=200\mu \mathrm{s}$ and $\theta =1.08\pi$ (b) and $1.16\pi$ (c). For $N$ cycles of the “forward” block, we see the signal decay in red closed triangles. After $N=6$, the reversal sequence is applied for $N^{\prime }$ cycles (green open triangles), where we expect an echo to appear at $N^{\prime }=N=6$ (filled point and arrow). (d), (e) DTC echoes for $N=\{2,6,10\}$ (open blue circles, green triangles, yellow diamonds), where we show the absolute values of each signal for easier inspection. Expected echo locations are marked with filled points and arrows. In (b)–(e), blue dots show the DTC signal decay for $\theta \approx \pi$.

Figure 11

$S\left(t\right)$ for $\theta =\pi$, at $T=200\mu \mathrm{s}$, identical to the blue dots in Figs. 10–10. This decay cannot be explained by dipolar interactions in an ideal delta-function pulse model.

Figure 12

(a) Envelope of a $\pi$ pulse applied at frequency ${\omega }_0$, as measured by a pickup coil placed near the resonator. The phase transients (small out-of-phase signal, dark green) are much smaller than the in-phase (red) pulse amplitude. (b) Cumulative integral of the out-of-phase pulse amplitude, scaled by setting the integral of the in-phase signal to ${180}^{\circ }$. Each transient produces less than a degree of rotation. (c) Incorporating the effect of phase transients into the ${cos}^N\left(\epsilon \right)$ decay model (black, identically 1 for $\theta =\pi$) leads to a modified decay model (dashed line). The modified model decays too slowly to account for the decay envelope of the measured $\left|S\right(t\left)\right|$ at $\theta =\pi$, shown here for $T=200\mu \mathrm{s}$ (blue circles).

Figure 13

(a) Nutation of ${}^{31}\mathrm{P}$ (black dots), after removing the effect of the homonuclear dipolar coupling. We do this by dividing the measured nutation data (not shown) by the simulated Hahn echo, scaled by a factor of 2 in time (red solid line). A rotary echo experiment approximately undoes the effect of the $H_1$ inhomogeneity, and results in data (small blue dots) which match closely to the scaled Hahn echo simulation. We fit (gray) the modified nutation data to a model of a plausible field profile in the coil, with good results (a, inset). (b) The $H_1$ field profile corresponding to the fit parameters (from the fit to the scaled nutation data) for a two-Gaussian model; this represents a histogram of applied frequencies $\gamma H_1/2\pi$ during an applied pulse. (c) The histogram of applied frequencies $\gamma H_1/2\pi$ can be used to deduce a spread in the applied pulse angle $\theta$, which we use to modify the original product-of-cosines decay model (black, identically 1 for perfect $\pi$ pulses) to the corrected decay model (gray dashed line). Including the effects of both $H_1$ inhomogeneity and phase transients (red dotted-dashed line) only slightly modifies the modeled decay envelope. The magnitude of $S\left(t\right)$ at $\theta \approx \pi$ and $T=200\mu \mathrm{s}$ (blue circles) decays much faster than the product-of-cosines model, even after including the effect of these pulse imperfections.

Figure 14

(a) Significant differences in the decay rate between sequences that are identical in the delta-function $\pi$-pulse approximation. At $\tau =20\mu \mathrm{s}$, the pulse sequences $\{X,X\}$ (black open squares) and $\{Y,Y\}$ (red open circles) produce very different lifetimes than $\{X,Y\}$ (green open triangles). The effect of the internal Hamiltonian during the pulse time $t_p$ creates differences between these sequences, which gives the latter sequence a much longer lifetime (see text). Because the signal is only observed every two cycles, the oscillations in the signal are not seen here. (b) Results of the pulse sequence ${\{\tau -X_{\pi }-Y_{\pi }-X_{\pi }-Y_{\pi }\}}^N$ (closed blue diamonds), which again exhibits an extended lifetime compared to the original DTC sequence, even at long $\tau$. (c), (d) The difference in lifetimes as a function of absolute time is significant, but displaying the pulse sequences as functions of the number of applied $\pi$ pulses or repeated blocks shows even more dramatic differences. Here, we define $T^{*}$ as the shortest repeated period, ignoring the phase of the pulses. For ${\{\tau -X_{\pi }-Y_{\pi }-X_{\pi }-Y_{\pi }\}}^N, T^{*}=\tau +4t_p$, while for $\{X,Y\}, T^{*}=\tau +t_p$.

Figure 15

Computed ${}^{31}\mathrm{P}$ spectra $S\left(\nu \right)$, with marked rms frequencies $W/2\pi$. (a) The pure ${}^{31}\mathrm{P}-{}^1\mathrm{H}$ dipolar spectrum $S^{\text{P,H}}\left(\nu \right)$ (solid line), with marked rms frequency (dashed line) $W^{\text{P,H}}/2\pi =3500\mathrm{Hz}$. (b) The pure ${}^{31}\mathrm{P}-{}^{14}\mathrm{N}$ dipolar spectrum $S^{\text{P,N}}\left(\nu \right)$ (red solid line) is narrower than the pure ${}^{31}\mathrm{P}-{}^{31}\mathrm{P}$ dipolar spectrum $S^{\text{P,P}}\left(\nu \right)$ (green solid line). The rms frequencies (dashed lines) $W^{\text{P,N}}/2\pi =97\mathrm{Hz}$ and $W^{\text{P,P}}/2\pi =508\mathrm{Hz}$ are marked. Not shown here are $S^{\text{P,PN}}\left(\nu \right)$ [similar to Fig. 2, red circles] with rms frequency $W^{\mathrm{P},\mathrm{PN}}/2\pi =517\mathrm{Hz}$, and $S^{\text{P,HPN}}\left(\nu \right)$ [see Fig. 2, blue squares] with rms frequency $W^{\text{P,HPN}}/2\pi =3538\mathrm{Hz}$.

Figure 16

(a) Lorentzian decay constant $N^{*}$ as a function of $\theta$ used in this example. (b) Crystalline fraction calculated for the distribution in (a), using $N=1$–128 (blue dashed line), $N=1$–50 (red solid line), and $N=1$–20 (green dotted line).