More broadly, our work sheds light on several fundamental challenges for the realization of spin squeezing in optically active, solid-state sensors. Perhaps most formidable is the presence of strongly coupled clusters of spins, which contribute an anomalously fast decay to the collective spin length. Our approach to overcoming this challenge, namely, spectrally resolved lattice engineering, is widely applicable to other disordered dipolar platforms, including magnetic atoms and polar molecules in optical lattices with low filling59, as well as ultracold gases of Rydberg atoms60. Finally, our work also introduces and implements a new technique for reading out the quantum variance of a state, without the need for projection-noise-limited measurements.
Our sample is grown using plasma-enhanced chemical vapour deposition (PECVD) with nitrogen delta-doping. The growth substrate is sliced from an electronic grade (001)-oriented CVD substrate (Element Six) along the (111) plane and then polished by Applied Diamond to a surface roughness less than 500 pm. Before growth, the miscut angle is measured using X-ray diffractometry rocking curves about the (111) omega peak and surface roughness is verified with AFM. PECVD growth is performed in a SEKI SDS6300 reactor at 770 °C, which finally produces a 99.998% C isotopically purified epilayer. In the nitrogen delta-doping process, N gas is introduced into the PECVD chamber for 1 min. After growth, secondary ion mass spectroscopy (SIMS) is performed in-house with conditions outlined in refs. to estimate the isotopic purity, epilayer thickness and the properties of the delta-doped layer (Extended Data Fig. 1a, i-iii).
Because the growth substrate is a diamond with natural carbon abundance while the epilayer is isotopically purified, the epilayer thickness is given by the low C density layer, that is, about 270 nm (Extended Data Fig. 1a, iii). A significantly higher N density is observed in the delta-doped layer 125 nm below the diamond surface compared with the background within the epilayer (Extended Data Fig. 1a, ii). The thickness of the delta-doped layer extracted from SIMS is estimated from the full-width at half-maximum (FWHM) of the peak to be . The error is calculated from the C SIMS step function used to calculate the growth rate of the epitaxial layer, and it represents the data broadening due to surface roughness. Layer thicknesses extracted from SIMS should be taken only as an estimate and represent an upper bound due to the roughness-induced broadening.
To further increase the density of NV centres, the sample is then irradiated with electrons at 200 keV from a transmission electron microscope (TEM, ThermoFisher Talos F200X G2 TEM). The irradiation dosage is 1.29 × 10 e cm for the spot used in this experiment. The sample is then annealed at 400 °C for 2 h and 850 °C for 4 h in a vacuum furnace (4 × 10 torr) to promote vacancy diffusion and NV formation. Finally, the sample is further cleaned by boiling in a tri-acid solution (1:1:1 HSO:HNO:HClO) and annealed at 450 °C in air for 4 h to remove surface contaminants and form oxygen termination for the diamond surface. Further growth details can be found in ref. .
The (111)-oriented diamond is crucial for this experiment. The NV-NV interaction is
where the interaction strength , in general, has an angular dependence on θ, that is, the angle between the NV quantization axis and the vector connecting two spins. In a conventional three-dimensional (3D) sample or even a delta-doped 2D (001)-oriented diamond sample, the NV-NV interaction J averages to zero because of the angular dependence of each of the four possible NV lattice orientations. The vanishing average interaction cannot generate twisting dynamics, which makes spin squeezing impossible. In our 2D (111)-oriented sample, the NV group quantized in the out-of-plane direction satisfies θ = π/2, so . In this case, the NV-NV interactions all have the same sign, providing a non-vanishing average interaction for twisting dynamics.
The delta-doped sample is measured with a custom-built scanning confocal microscope at room temperature. For optical initialization and readout, 75 μW of 532 nm light is sent through an air objective (Nikon CFI Plan Apo 60X, numerical aperture 0.95). The NV fluorescence is separated from the green excitation light with a dichroic filter and collected on fibre-coupled single-photon counters (Excelitas SPCM-AQRH-12-FC). A typical scan of our irradiated spot is shown in Extended Data Fig. 1b. A magnetic field is generated by three sets of orthogonal electromagnets with tunable current sources. The field is aligned to the out-of-the-plane direction of the sample and set to 393 G to maximize the nuclear spin hyperpolarization. Microwave pulses at 1.742 GHz (|m = 0⟩ ↔ |m = -1⟩ transition) and 3.998 GHz (|m = 0⟩ ↔ |m = 1⟩ transition) at 393 G are generated by an arbitrary waveform generator (AWG, Tektronix AWG7122C), amplified by a broadband microwave amplifier (Mini-Circuits ZHL-50W-63+) and delivered through a stripline to achieve a typical Rabi frequency of about (2π) × 25 MHz. Various microwave filters and circulators are inserted into the signal chain to suppress reflections and improve the microwave pulse quality. To generate both strong and weak pulses simultaneously, we use two channels of the AWG and combine the signals using a resistive power combiner (Pulsar RP2-05-411). The amplitude of the microwave pulse is further tuned with a programmable attenuator (Mini-Circuits RUDAT-6000-90).
At the magnetic field B = 393 G used in our experiment, the nuclear spin polarization is 89(1)%, measured by relative amplitudes of the optically detected magnetic resonances (ODMR) (Extended Data Fig. 1e). In our experiments, we work with the majority m = +1/2 nuclear spin subgroup and shelve the residual m = -1/2 population in the state |m = 1⟩ to avoid off-resonant driving of these spins. The Rabi frequency of the shelving pulse is chosen to be about (2π) × 600 kHz, which is smaller than the hyperfine splitting and larger than the inhomogeneous broadening. The ODMR spectra before and after shelving are shown in Extended Data Fig. 1e.
Substitutional nitrogen (P1) centres are the densest defects in the delta-doped layer.
These paramagnetic spin-1/2 impurities generate a quasi-static onsite field for each NV centre, which we cancel with dynamical decoupling pulses (see the next section).
We estimate the P1 density from ODMR spectroscopy, using the linewidth of the NV |m = 0⟩ → |m = -1⟩ transition (Extended Data Fig. 1e). From the FWHM, the NV-P1 interaction strength is around (2π) × 360 kHz (setting ħ to 1 throughout), leading to an estimated P1 density of about 75 ppm nm.
To characterize the NV areal density in the irradiation spot used for our experiments, we measure the NV decoherence timescale under an XY-8 pulse sequence. To show that this timescale is limited by interactions between NV centres rather than by the stronger interactions between NV and P1 centres, we measure as a function of the interpulse spacing τ (Extended Data Fig. 1c,d). Because smaller interpulse spacings result in improved decoupling of the NV-P1 interactions, we observe an increase in the decoherence time when the interpulse spacing τ decreases. We choose an optimal interpulse spacing τ = 50 ns for our measurements.
We then extract the NV density by comparing cluster DTWA numerical simulations with the measured decoherence (Extended Data Fig. 1f,g), from which we infer an NV density of 8 ppm nm for the out-of-the-plane NV group. The total NV density in the irradiation spot is thus 32 ppm nm. The NV-NV interaction strength in the out-of-plane group, with which we work exclusively, is then about (2π) × 80 kHz.
We first discuss the pulse sequence used to measure the mean-field interactions shown in Fig. 1. Then, we elaborate on the pulse sequence used to generate and the readout squeezing. As shown in Extended Data Fig. 2, the full experiment pulse sequence comprises a lattice-engineering sequence, followed by squeezing generation and readout steps. Each is discussed in detail in the following sections.
After optically pumping to the |m = 0⟩ state, we initialize a spin-polarized state, |x⟩, which is offset above (or below) the x-axis by an angle φ (Fig. 1c). This state has an expectation value and will precess about the z-axis such that and . We measure ⟨S⟩ by applying a rotation and ⟨S⟩ by applying a rotation before reading out the population difference between the spin states. We then calculate the precession angle in the equatorial plane as .
The first step in the experiment is to prepare a lattice-like geometry with relatively homogeneous spin-spin couplings by removing closely spaced NV clusters from the system. Removal of these clusters is accomplished by harnessing their most prominent feature, that is, the large interaction energy induced by the small distance between spins. The Ising component of this interaction shifts the bare splitting of the |m = 0⟩ → |m = ±1⟩ transition for a spin i by ±J, defined as
For sufficiently large system sizes, the mean-field coupling J can be treated as a continuous random variable J with probability P(J). For disordered spins in 2D, clusters -- in particular, dimers formed by nearby pairs of spins -- give rise to the heavy tail in the distribution P(J) (Fig. 3a). Because the dimer energy is larger than the average mean-field interaction , dimers can be isolated from the remainder of the system using a frequency-selective control pulse in either the z or x basis as described below.
Our first method for removing dimers from the system is frequency-resolved shelving (Extended Data Fig. 2b, i). After optically initializing the spins in the state |m = 0⟩, we apply a weak π-pulse whose amplitude Ω is set to the average NV coupling strength . This pulse drives typical NV centres to the |m = -1⟩ state while dimers remain in |m = 0⟩. Because the relevant NV-NV interactions are dominated by larger NV-P1 interactions, we decouple the latter by interleaving XY-8 pulses with the weak control pulse. After the weak π-pulse, we apply two strong π-pulses: the first shelves the dimers in the |m = 1⟩ state, and the second transfers the target spins back to the |m = 0⟩ state. At the end of the state preparation sequence, the dimers are shelved in |m = 1⟩, and the remaining spins are prepared in |m = 0⟩.
Because of the large difference between the amplitudes of the XY-8 and frequency-selective control pulses, we use two channels of the AWG and attenuate the microwaves separately before sending them to the microwave amplifier. We calibrate the relative phase between the microwaves generated by the two channels such that the +x-axis defined by the XY-8 pulses coincides with that defined by the weak continuous drive.
In principle, the dimers shelved in |m = 1⟩ can still interact with the rest of the system. However, the XY-8 pulse sequence is always applied during the subsequent quench dynamics (see sections 'Generation of spin squeezing' and 'Readout of quantum spin projection noise'), so the Ising interaction between the dimers and the rest of the system is decoupled. Furthermore, the NV-P1 Ising interactions suppress NV spin exchange within the {|m = 0⟩, |m = 1⟩} subspace. As a result, the shelved spins do not participate in the squeezing or readout dynamics.
Several effects can limit the performance of this protocol, including XY-8 pulse imperfections and dimer dynamics during the slow pulse. In practice, we optimize the shelving parameters empirically to maximize the average twisting rate compared with the initial decay of the collective spin ⟨S⟩.
Our second method for removing NV dimers from the system is adiabatic depolarization (Extended Data Fig. 2b, ii). After optical initialization, we prepare the spin-polarized state |x⟩ with a pulse. We then apply a weak, time-dependent transverse field h(t) while simultaneously performing XY-8 decoupling pulses. The initial field h(0) is much larger than NV dimer interactions, so the spin-polarized state |x⟩ is an eigenstate. Then, we ramp the field down to a final value to depolarize dimers by two complementary mechanisms. First, spins with J ≫ h undergo fast decoherence, whereas spins with J ≪ h remain spin-locked along the +x quantization axis. Second, the temporal profile of the field h(t) improves the efficiency of depolarization by constraining the spins to adiabatically follow the eigenstates of the Hamiltonian. As a result of these two effects, this protocol removes the S polarization for the dimers without significantly depolarizing the other spins.
We determine the ramp-down profile of the transverse field h(t) by optimizing for adiabaticity. Specifically, for a two-body XXZ interaction with a strong transverse field h ≫ J,
The dimer eigenstate to first order in perturbation theory in J is
Using equation (7), we can rewrite the adiabatic condition in terms of the transverse field. The optimal ramp-down profile h(t) then obeys the differential equation , whose solution is given by
where k is the parameter controlling the adiabaticity. Note that we assume that h(t) is much larger than J in this derivation, but the same ramp-down profile still holds even when h(t) is no longer large compared with h(t) as verified by the numerical diagonalization of the Hamiltonian. We empirically optimize the parameter k and the total ramp duration by maximizing the twisting rate χ compared with the initial decay of the spin during the squeezing dynamics. After the state preparation protocol, we rotate the remaining polarized spins back to |m = 0⟩ and wait for a time to fully decohere the dimers.
After lattice engineering, we apply a pulse to prepare the spin-polarized state |x⟩ for the remaining system spins. To generate the squeezed state, we allow the system to evolve under the native NV-NV interactions H for a time t by applying an XY-8 pulse sequence with interpulse spacing τ = 50 ns to suppress the larger NV-P1 interactions (Extended Data Fig. 2c).
To measure the variance of the operator , we rotate the squeezed state about the +x-axis by an angle θ with an X pulse, such that the S operator for the pre-rotated state becomes the S operator for the post-rotated state . We then measure the decay profile S(t) for the post-rotated state, while continuously decoupling NV-P1 interactions with an XY-8 pulse sequence with interpulse spacing τ = 50 ns (Extended Data Fig. 2d).
A spin-polarized state |x⟩ evolving under H decays as a stretched exponential . In particular, a disordered dipolar ensemble in 2D decays with stretch power p = 2/3, whereas spins on a regular lattice decay with p = 2 (ref. ). The stretch power p thus characterizes the spatial distribution of the spins and shows the efficacy of our dimer removal protocols (see section 'NV cluster removal by lattice engineering'). The removal of dimers imposes a minimum distance between the remaining spins, such that the local geometry is lattice-like while the global geometry is still disordered. We thus expect Gaussian decay (p = 2) at early times, crossing over to disordered dynamics (p = 2/3) when the effective interaction length scale is much larger than the minimum spacing .
To verify the above intuition, we simulate the decoherence profile as a function of (Extended Data Fig. 3). Beginning with spins randomly distributed in a 2D plane with areal density 8 ppm nm, we then remove spins that have a neighbour within . When no spins are removed from the system, the stretch power of the decay profile is p = 2/3 for three decades in time (Extended Data Fig. 3, purple). As the minimum distance is increased (Extended Data Fig. 3, blue to green curves), we observe the predicted crossover of the stretch power from p = 2 (early time lattice regime) to p = 2/3 (late time disordered regime). Furthermore, the crossover timescale t increases as increases, which is consistent with the analytical prediction (Supplemental Information).
In the analysis below, we describe how we extract an effective minimum radius imposed by the shelving and depolarization protocols from XY-8 decoherence timescales. These models for spectral tailoring of the spatial geometry are then used to compute the input state for the squeezing and readout numerics described in section 'Mapping between decay time scale and quantum variance'.
To quantitatively model dimer shelving, we first calculate the distribution P(J) for NV centres in a 2D plane. With the sequence in section 'Dimer removal method I', the probability that the frequency-selective control π-pulse drives a spin to the |m = -1⟩ state is
As shown in Extended Data Fig. 4a (ii and iii), we measure the XY-8 decoherence profile after the shelving protocol, and compare the experimental result (black points) with cluster DTWA numerical simulations (coloured curves) assuming the above shelving probability and areal NV density 8 ppm nm. Different curves correspond to simulated decoherence profiles assuming different shelving radii r. We thus extract a dimer removal radius r = 7 nm and r = 6 nm for the shelving protocol for the pink and purple data, respectively, shown in Fig. 4c.
To quantitatively model dimer depolarization, we estimate the depolarization rate for dimers with spacings r by calculating the expectation value ⟨S⟩ for the dimer eigenstate. Assuming a transverse field strength , the dimers whose inter-spin distance r < r are mostly removed. Here, r parametrizes the dimer removal radius for the depolarization protocol. Rewriting the dimer Hamiltonian (equation (6)) as
As shown in Extended Data Fig. 4b (ii and iii), we measure the XY-8 decoherence profile after the depolarization protocol and compare the experimental result with the cluster DTWA numerical simulation assuming the above depolarization probability. The curves are the decoherence profiles assuming different cut radius r. The NV density is fixed at 8 ppm nm. We extract a dimer removal radius r = 14 nm for the depolarization protocol corresponding to the data in Fig. 4f.
Briefly, we outline our procedure for extracting the variance Var(S) from readout quench dynamics ⟨S(t)⟩. Each step is then detailed individually in the following sections. After evolving the input state for a time t to generate squeezing, we measure the spin length decay ⟨S(t)⟩ as a function of the readout time t for a set of rotation angles θ. We then shift each measured curve by an offset time t to obtain (see section 'Offset time'). After shifting, we fit the set of curves with a stretched exponential decay within a window . The timescale T(θ, t) is fit for each curve individually, and the stretch power p is globally optimized over all curves (see section 'Curve fitting'). Next, we calculate a one-to-one mapping T(θ, t) → Var(S(t)) numerically (see section 'Mapping between decay time scale and quantum variance'). We finally use this mapping to calculate the normalized squeezing parameter for our data (see section 'Extraction of squeezing parameter').
The essence of our detection protocol relies on the insight that the timescale T for the decay of the spin length ⟨S⟩ is inversely related to the variance Var(S). To fit the decay timescale correctly, we want to start with maximal ⟨S⟩; in other words, starting the fit at some arbitrary time during the decay of ⟨S⟩ will not yield the correct timescale, especially when compared across measurements. The ideal scenario is shown in Fig. 2c: three fiducial states with ⟨S⟩ = ⟨S⟩ = 0, maximal⟨S⟩ ~ S, and different Var(S) are allowed to evolve under H for a time t. During these readout quench dynamics, the spin vector ⟨S⟩ for the state with the largest variance Var(S) decays the fastest, and that with the smallest variance Var(S) decays the slowest. Because each state satisfies maximal ⟨S⟩ at t = 0, the curves S(t) can be fitted for T directly starting at t = 0. The result is shown in Fig. 2d.
In practice, we cannot achieve the above ideal scenario in the experiment. After squeezing generation for time t and a rotation X, states with different Var(S) are prepared at t = 0. However, these states do not generically satisfy maximalS. Rather, the corresponding fiducial state with maximal ⟨S⟩ exists at some offset time t that depends on θ. To accurately fit T across different rotation angles θ, the readout time needs to be shifted t - t. This idea is shown in Extended Data Fig. 6 for four example rotation angles θ. From an initial spin-polarized state, a squeezed state is generated after time t. Then, a rotation X is applied, and the readout quench begins. Two rotation angles θ do yield a state of maximal S at t = 0; the first rotates the anti-squeezed quadrature to +z (red), and the second rotates the squeezed quadrature to +z (blue). For all other rotation angles θ, the fiducial state occurs at some other time (which can be either positive or negative). Take the θ = 0 (black): because no rotation is applied, Var(S) is the same as the initial spin-polarized state. The fiducial state for θ = 0 is then exactly the spin-polarized initial state, with t = -t. In other words, the T fitting for the black curve should start at t - t = -t, not at t = 0. Similarly, to accurately fit T for the purple curve, we should fit from t - t = t.
In general, the offset time t(θ, t) depends on both the squeezing duration t and the rotation X. We determine the offset for each measurement ⟨S(t)⟩ analytically (Supplementary Information) and benchmark against numerics for the disordered system, finding good agreement (Extended Data Fig. 5c). We then shift each measurement ⟨S(t - t)⟩ to obtain and fit the decay timescale as described in the next section.
We fit each shifted decoherence profile to the functional form , where T(θ) and A(θ) are local fit parameters and the stretch power p is optimized globally across all measurements {θ, t}. We fix a fitting window . We choose the lower limit of the fitting window by ensuring that the decay profile for each rotation angle θ has the same number of data points in the fitting window. In the experiment, data are only available for the detection quench time t > 0, that is, . Because t can be negative, we set the lower limit of the fitting window , where the maximum is reached at θ = 0.
We choose the upper limit of the fitting window by evaluating how well the decay profile can be effectively described as a stretched exponential decay with a single stretch power p. As the decay profile has a crossover behaviour for the stretch exponent (see section 'Stretch power as a measure of spatial geometry'), the crossover timescale sets . We perform the fits for a range of reasonable fitting windows with (Extended Data Fig. 5d,e) and include the resulting fluctuations in the error bars of Fig. 4.
To map a measured decay time T to Var(S), we establish a one-to-one mapping between the extracted decay timescale and the quantum variance. A quantitative study for the mapping is provided as follows.
To build the mapping, we place N = 80 spins in a 2D plane and remove dimers according to r or r as detailed in section 'Stretch power as a measure of spatial geometry'. We perform cluster DTWA numerics with 50 realizations of positional disorder to simulate quench dynamics under H: generation of squeezing for a time t, rotation X and subsequent decay ⟨S(t)⟩. The input parameters for our simulations are determined from independent measurements (Extended Data Table 1). Thus, the numerics used to generate the mapping contain no free parameters.
A comparison between the resulting raw numerical and experimental data is shown in Extended Data Fig. 5a. We note that the solid curves ⟨S(t)⟩ are directly output from our numerics -- that is, there are no fit parameters. We shift and fit the numerical and experimental data following the procedures in sections 'Offset time' and 'Curve fitting'. We can thus compare the timescales T obtained by the numerics and the experiment, as shown in Extended Data Fig. 5d-f. From the same numerical data, we compute the variance Var(S) as a function of t and θ. Finally, plotting Var(S) and T exhibits a clear one-to-one correspondence, with consistent results across the measured squeezing durations from 0 μs to 4 μs. From the plot in Extended Data Fig. 5i, we can thus measure the variance Var(S) by the decay timescale T.
With the protocol explained in section 'Mapping between decay time scale and quantum variance', we extract the relative variance Var(S(t))/Var(S(t)) by measuring the S decay in the readout quench. Meanwhile, we directly measure S(t)/S(0) in the generation quench.
We emphasize that the normalized squeezing parameter is not equivalent to the Wineland squeezing parameter
except when S(0) = N/2 (assuming spin-1/2), because it normalizes away any imperfections in the initial state preparation (arising in our case from imperfect spin polarization). Hence, it cannot be used as an entanglement witness. Nevertheless, the normalized squeezing parameter does capture the metrological enhancement that would be obtained if ideal Ramsey spectroscopy were performed using the state at t = 0 compared to t = t (ref. ). The normalized squeezing parameter captures the squeezing dynamics under H while the imperfect state preparation, for example, finite initial polarization, is normalized out.