Coherent dynamics of strongly interacting electronic spin defects in hexagonal boron nitride

Abstract

Optically active spin defects in van der Waals materials are promising platforms for modern quantum technologies. Here we investigate the coherent dynamics of strongly interacting ensembles of negatively charged boron-vacancy (${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$) centers in hexagonal boron nitride (hBN) with varying defect density. By employing advanced dynamical decoupling sequences to selectively isolate different dephasing sources, we observe more than 5-fold improvement in the measured coherence times across all hBN samples. Crucially, we identify that the many-body interaction within the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble plays a substantial role in the coherent dynamics, which is then used to directly estimate the concentration of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. We find that at high ion implantation dosage, only a small portion of the created boron vacancy defects are in the desired negatively charged state. Finally, we investigate the spin response of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ to the local charged defects induced electric field signals, and estimate its ground state transverse electric field susceptibility. Our results provide new insights on the spin and charge properties of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, which are important for future use of defects in hBN as quantum sensors and simulators.

Introduction

Solid-state point defects with optically addressable electronic spin states have become some of the most fertile playgrounds for new quantum technologies^{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18}. Significant recent progress has been made in creation and control of such spin-active quantum emitters in atomic-thin van der Waals materials. The two-dimensional (2D) nature of the host materials can enable seamless integration with heterogeneous, optoelectronic, and nanophotoic devices, providing a pathway to investigating light-matter interactions at the nanoscale^19,20,21,22.

From a wide range of contestant spin defects in 2D materials, the negatively charged boron vacancy center, ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, in hexagonal boron nitride (hBN) has particularly attracted substantial research interest in the past few years^{23,24,25,26,27,28,29,30,31,32}. Importantly, it has been demonstrated that the spin degree of freedom of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ can be optically initialized and readout, as well as coherently manipulated at room temperature. Compared to conventional spin qubits in three-dimensional materials, such as nitrogen-vacancy (NV) center in diamond, ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ features several unique advantages in quantum sensing and simulation.

From the perspective of quantum sensing, the atomically-thin structure of hBN can allow the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ sensor to be positioned in close proximity with the target materials, facilitating the imaging of inter-facial phenomena with unprecedented spatial resolution and sensitivity^25,33,34,35. Moreover, since hBN has been widely employed as the encapsulation and gating dielectric material in 2D heterostructure devices, introducing the embedded ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ sensors does not require any additional complexity in the fabrication process^{36,37,38,39,40}. On the quantum simulation front, the ability to prepare and control strongly interacting, two-dimensional spin ensembles opens the door to exploring a number of intriguing many-body quantum phenomena^41,42,43. For instance, dipolar interaction in 2D is particularly prominent from the perspective of localization and thermalization, allowing one to experimentally investigate the effect of many-body resonances^{44,45,46,47,48,49,50,51}.

${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ in hBN, like solid-state spin defects in general, suffers from decoherence. To this end, research effort has been devoted to characterizing the coherence time of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. However, the measured spin echo timescale, ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}}$, in several studies varies from tens of nanoseconds to a few microseconds^24,52,53,54. This immediately begs the question that where does such discrepancy originate from, and what are the different decoherence mechanisms in dense ensemble of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ?

In this letter, we present three main results. First, we introduce a robust differential measurement scheme to reliably characterize the spin coherent dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble (Figs. 1 and  2). We observe spin-echo ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}} \, \approx \, 70$ ns across three hBN samples with distinct ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities (created via ion implantation with dosages spanning two orders of magnitude), consistent with the expectation that the spin-echo coherence time is dominated by the Ising coupling to the nearby nuclear spin and dark electronic spin bath^52,55. By applying a more advanced dynamical decoupling sequence, XY-8, to better isolate ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ from the bath spin environment^56,57,58, we observe substantial extensions in the measured coherent timescales, ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$. Interestingly, the extracted ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ decreases with increasing ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, indicating that the dipolar interaction within the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble is critical for understanding the coherent dynamics. To further corroborate this, we utilize the DROID pulse sequence to decouple the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ −${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ dipolar interaction^59,60, and achieve an additional ~2-fold improvement in the measured coherence time, ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$. Second, by comparing the experimentally measured ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ to numerical simulations, we directly esimtate the spin density of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ across three hBN samples. We find that the ratio of negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ to total created boron vacancy defects (V_B) decreases significantly with increasing ion implantation dosage (Fig. 3). Third, based on the extracted ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, we introduce a microscopic model of local charges surrounding a spin defect to account for the observed energy splitting between $\left|{m}_{s}=\pm 1\right\rangle$ states at zero magnetic field^61,62, and estimate the transverse electric field susceptibility of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ to be around d_⊥ ≈ 40 Hz/(V ⋅ cm⁻¹) (Fig. 4).

Fig. 1: Spin dynamic of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble.

a Schematic of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin ensemble (red spins) inside hBN crystal lattice (Nitrogen--blue; Boron--white); z is defined along the c-axis (perpendicular to the lattice plane). x and y lie in the lattice plane, with x oriented along one of the three ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ Nitrogen bonds. Two types of decoherence sources are presented here for ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin ensemble: the Ising coupling (gray wavy lines) to the bath spins (gray), and the dipolar interaction within ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ themselves (red wavy lines). b Energy level diagram of the defect spin ground-state. In the absence of any external perturbation, the $\left|{m}_{s}=\pm 1\right\rangle$ states are degenerate and separated by D_gs ≈ 3.48 GHz from the $\left|{m}_{s}=0\right\rangle$ state. Under an external magnetic field B along the c-axis of hBN, the degeneracy between $\left|{m}_{s}=\pm 1\right\rangle$ states are lifted via the Zeeman effect, with a splitting ∝ 2B. We choose $\left|{m}_{s}=0\right\rangle$ and $\left|{m}_{s}=-1\right\rangle$ states as our two-level system. c Experimental pulse sequences for XY-8 (top) and DROID (bottom). The rotations along the positive x and y axes are plotted above the line, while the rotations along the negative axes are plotted below the line. d Differential measurement sequence for spin echo. I: 20 μs wait time to reach charge state equilibration. II: 10 μs laser pulse to initialize the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin to $\left|{m}_{s}=0\right\rangle$, with the reference signal, S_R(t), collected at the end of the laser pulse. III: microwave wave pulses for spin echo measurement; for the bright signal, a final $\frac{\pi }{2}$ pulse along the -y axis is applied; while for the dark signal, a final $\frac{\pi }{2}$ pulse along the +y axis is applied to rotate the spin to an orthogonal state. IV: laser pulse to detect the spin state. e Spin echo measurement on sample S3 at two different laser powers. Without differential measurement, the measured signal, S_B/S_R exhibits a laser power dependence which comes from charge relaxation dynamics (inset). Using differential measurement, the measured contrast, C(t), is independent of the laser power. Error bars represent 1 s.d. accounting statistical uncertainties.

Fig. 2: Spin coherent and relaxation dynamics.

a The spin coherent and relaxation timescales measured on sample S3 with the highest ion implantation dosage. Dashed lines are data fitting with single exponential decays. Error bars represent 1 s.d. accounting statistical uncertainties. b The extracted coherence timescales T₂ and relaxation timescales T₁ for the three hBN samples. Error bars in time represents 1 s.d. accounting fitting error, and error bars in implantation dosage represents the current fluctuation induced uncertainties (~10%) during implantation.

Fig. 3: Characterizing ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density.

a Comparison between the experimentally measured and numerically simulated coherent timescales, T₂, for DROID and XY-8 pulse sequences. The solid lines show the timescales extracted from simulations with error bars plotted as semi-transparent colored areas. To determine ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities for the three hBN samples, we minimize the relative squared residuals of ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ between simulations and experiments. Inset: fluorescence counts (PL) versus extracted densities after contrast adjustment (see Methods). Error bars in coherence time account for 1 s.d. of fitting error, and error bars of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities represent the range of densities whose residuals lie within 5% (see Methods). b The measured ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ charge state ratio $\eta={\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}/{\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}}$ for three hBN samples with different ion implantation dosages. Error bars account for both the current fluctuation induced implantation uncertainties and the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities uncertainties from (a).

Fig. 4: Imaging the local electric field signals.

a Top: Schematic depicting the charged defects environment surrounding a ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ electronic spin. Nearby negatively and positively charged defects create a local transverse electric field E on ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. Bottom: Energy level diagram of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin state in the presence of the electric field: the E-field mixes the $\left|{m}_{s}=\pm 1\right\rangle$ states to new eigenstates $\left|\pm \right\rangle$, with a splitting, δ ∝ 2d_⊥E. b Measured ESR spectra of sample S1 and sample S3 at zero magnetic field. Dashed lines are the simulated results from our microscopic charged model using d_⊥ = 40 Hz/(V cm⁻¹) and ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities extracted from the previously measured coherent dynamics. Fluorescences are shifted vertically for comparison. c Numerically simulated ESR splitting δ using different electric susceptibilities, d_⊥. The red dots are the experimental results with error bars representing the range of densities whose residuals lie within 5% (see Methods).

Results

Experimental system

To investigate the coherent spin dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble at various defect densities, we prepare three hBN samples with different implantation dosages. Specifically, we irradiate hBN flakes (thickness ~ 100 nm) using 3 keV He⁺ ion beams with dose densities, 0.30 ± 0.03 nm⁻² (sample S1), 1.1 ± 0.1 nm⁻² (sample S2), and 10 ± 1 nm⁻² (sample S3), respectively, to create ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ defects^27,53. Here error bars on the implantation dosages account for the current fluctuations during the implantation process. We remark that, given an ion implantation dosage, the total created V_B concentration can be estimated via SRIM simulation (see Methods)⁶³, yet the actual density of the negatively-charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ has remained unknown.

The ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ center has a spin triplet ground state ($|{m}_{s}=0,\pm 1\rangle$), which can be initialized and read out via optical excitation and coherently manipulated using microwave fields^23,30. In the absence of any external perturbations, the $\left|{m}_{s}=\pm 1\right\rangle$ states are degenerate and separated from $\left|{m}_{s}=0\right\rangle$ by D_gs ≈ 3.48 GHz (Fig. 1b). In the experiment, we apply an external magnetic field B ≈ 250 G along the c-axis of the hBN lattice to separate the $\left|{m}_{s}=\pm 1\right\rangle$ states via the Zeeman effect and isolate an effective two-level system $\left|{m}_{s}=0,-1\right\rangle$. A microwave field is used to coherently manipulate the spin ensemble with a Rabi frequency Ω ≈ 83 MHz (π-pulse length τ_π = 6 ns). We note that such a strong Rabi drive is crucial for the high fidelity control of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, as the spin transition is largely broadened by the hyperfine interaction to the nearby nuclear spin bath (see Methods).

Robust measurement scheme

To reliably probe the spin dynamics of a dense ensemble of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, we introduce a robust differential measurement scheme illustrated in Fig. 1d^64,65. Specifically, after letting the spin system reach charge state equilibration for 20 μs without any laser illumination (I), we apply a 10 μs laser pulse (532 nm) to initialize the spin state of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ (II), followed by the measurement pulse sequences (III). Taking spin echo coherent measurement as an example, we first apply a $\frac{\pi }{2}$-pulse along the y axis to prepare the system in a superposition state ${\otimes }_{i}\frac{{\left|0\right\rangle }_{i}+{\left|-1\right\rangle }_{i}}{\sqrt{2}}$, and then let it evolve for time t. A refocusing π-pulse along the x axis at time t/2 is used to decouple the spin ensemble from static magnetic noise. A final $\frac{\pi }{2}$-pulse along the -y direction rotates the spin back to the z axis for fluorescence detection (IV), and the measured photon count is designated as the bright signal, S_B(t). By repeating the same sequence but with a final $\frac{\pi }{2}$-pulse along the positive +y axis before readout, we measure the fluorescence of an orthogonal spin state to be the dark signal, S_D(t). The difference between the two measurements, C(t) = [S_B(t) − S_D(t)]/S_R(t), can faithfully represent the measured spin coherent dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, where S_R(t) is a reference signal we measure at the end of the initialization laser pulse (II).

Figure 1e shows the measured spin echo dynamics of the highest dosage hBN sample S3. We find that the measured fluorescence contrast, S_B(t)/S_R(t) [S_D(t)/S_R(t)], changes dramatically with different laser powers (inset), originating from the charge state relaxation dynamics after the laser pumping. This is particularly prominent at high laser power, where the optical ionization of the defect charge state is enhanced. This effect can lead to an artifact in the extracted spin echo timescales, which may explain the previous discrepancy in the measured ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}}$. However, the obtained fluorescence contrast from differential measurement, C(t), is consistent across different laser powers, enabling an accurate extraction of the spin coherent timescales.

A few remarks are in order. First, this differential measurement scheme has been widely employed in the studies of the dense ensemble of NV centers in diamond to counter the ionization process^{9,41,65,66,67}. Secondly, previous theoretical studies predict that the ionization of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ requires significantly higher energy ( ~ 4.46 eV) than the ionization of NV centers ( ~ 2.7 eV)^66,68,69. This may explain why our experimental observation that the two-photon ionization process for ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ only becomes evident under strong laser power ( ~ 10 mW); while the ionization of NV centers happens at ~ 10 − 20 μW laser^65,66. Third, we note that unlike neutral NV⁰ centers which emit fluorescence starting at 575 nm, neutral boron-vacancy ${V}_{{{{{{{{\rm{B}}}}}}}}}^{0}$ has not been directly observed from photo-luminescence signals. Therefore the proposed ionization process only offers a potential explanation of the experiment.

Coherent dynamics

Equipped with the robust differential measurement scheme, we now turn to the investigation of coherent dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble at various defect densities. The decoherence mechanism of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ consists of two major contributions: (1) the Ising coupling to the bath spins in the environment; (2) the dipolar interaction between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble themselves (Fig. 1a and Methods)⁵⁹. To isolate the effect of each component, we measure the coherent dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ using three different dynamical decoupling pulse sequences.

We start with the spin echo pulse sequence, which is commonly used to characterize the coherent properties of a quantum system. Spin echo can decouple the static components of the Ising coupling between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ and the spin bath. By fitting the measured dynamics to a single exponential decay, $\sim {e}^{-(t/{T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}})}$, we extract ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}} \, \approx \, 70$ ns across all three hBN samples (Fig. 2b). This observation indicates that the spin echo decoherence of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ is predominantly limited by the spin fluctuation within the environmental spin bath, which does not depend on the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration. Indeed, a previous study has shown that the Ising coupling to the local nuclear spin bath (nitrogen-14, boron-10, and boron-11), as well as the dark electronic spins, can accurately account for the measured spin echo timescales⁵².

Next, we apply a more advanced dynamical decoupling pulse sequence, XY-8, to better decouple the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble from the environment. Instead of a single refocusing π-pulse, XY-8 employs a series of π-pulses with alternating phases (Fig. 1c). We fix the time intervals between pulses, τ₀ = 4 ns, sufficiently smaller than the correlation timescale of the local spin bath (estimated from the spin echo timescale)^41,45. As a result, XY-8 is expected to further suppress the fluctuations within the local spin noise and improve the measured spin coherent timescales. This is indeed borne out by our data. As shown in Fig. 2, the extracted coherence times, ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$, are significantly extended in all three samples. In contrast to the previous spin echo measurement where ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}}$ does not depend on ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, here we observe that ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}=[250\pm 35]$ ns of sample S1 is longer than sample S3, ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}=[167\pm 10]$ ns. This suggests that ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$  − ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ interaction plays a key role in the measured XY-8 coherent timescales. Indeed, in XY-8 measurement, since the refocusing π-pulses flip all ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spins together, there is no suppression of the dipolar interaction between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ (see Methods).

To this end, we introduce DROID pulse sequence to further decouple the dipolar interaction within ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ themselves (Fig. 1c)⁵⁹. By applying a series of π/2 rotations along different spin axes to change the frames of interaction (also known as toggling frames), DROID modifies the dipolar Hamiltonian to an isotropic Heisenberg interaction, where the initial state, ${\otimes }_{i}\frac{{\left|0\right\rangle }_{i}+{\left|-1\right\rangle }_{i}}{\sqrt{2}}$, constitutes an eigenstate of the Heisenberg interaction, and consequently does not dephase (see Methods). As shown in Fig. 2, the measured coherent timescales, ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$, indeed exhibit an approximate two-fold increase compared to ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ across all three samples, agreeing with the cancellation of dipolar-induced decoherence.

Interestingly, we also observe that the spin relaxation time, T₁, and spin-locking time, ${T}_{1}^{\rho }$, both decrease with increasing ion implantation dosages (Fig. 2b). In principle, the dipolar interaction between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ will not lead to a decrease of T₁ due to the conservation of total spin polarization during the flip-flop process (see Supplementary Note 2.2). This T₁ related trend may be attributed to the presence of lattice damage during the implantation process or local charge state fluctuations⁶⁵. We note that the spin relaxation process will introduce an additional decay to the coherent dynamics. However, the measured T₁ and ${T}_{1}^{\rho }$ are much longer than T₂ across all three samples at room temperature (Fig. 2). Nevertheless, we fix the duration between the polarization (II) and the read-out (IV) laser pulses to account for the effect of T₁ relaxation on the T₂ measurement (see Methods).

Extracting ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density

The difference between ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ originates from the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$  − ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ dipolar interaction, which can be used to estimate the density of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ directly. In particular, by randomly positioning 12 electronic spins at different sampling concentrations, we construct the dipolar interacting Hamiltonian of the system,

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}=\mathop{\sum}\limits_{i < j}-\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{i,j}}{{r}_{i,j}^{3}}\left({S}_{i}^{z}{S}_{j}^{z}-{S}_{i}^{x}{S}_{j}^{x}-{S}_{i}^{y}{S}_{j}^{y}\right),$$

(1)

where J₀ = 52 MHz ⋅ nm³, ${{{{{{{{\mathcal{A}}}}}}}}}_{i,j}$ and r_i,j represent the angular dependence and the distance between the i^th and j^th${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spins, and $\{{S}_{i}^{x},{S}_{i}^{y},{S}_{i}^{z}\}$ are the spin-1/2 operators for i^th${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ center (see Methods). By evolving the many-body system under different pulse sequences and averaging the spin coherent signals across random spin positional configurations, we obtain the simulated results of the corresponding XY-8 and DROID coherent timescales (Fig. 3a, see Methods)^9,65,70. We observe from our simulation that both ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ indeed decrease with increasing ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, while ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ exhibits a longer timescale than ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ across the density range surveyed. By minimizing the relative squared residuals of ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ between simulation and experiment, we estimate the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration to be ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}^{{{{{{{{\rm{S1}}}}}}}}} \, \approx \, 123\begin{array}{c}+8\\ -8\end{array}$ ppm, ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}^{{{{{{{{\rm{S2}}}}}}}}} \, \approx \, 149\begin{array}{c}+25\\ -21\end{array}$ ppm, and ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}^{{{{{{{{\rm{S3}}}}}}}}} \, \approx \, 236\begin{array}{c}+35\\ -31\end{array}$ ppm. The discrepancy between the measured and simulated timescales may stem from imperfect spin rotations in the experiment, as well as finite-size effects from the simulations (see Methods). To further validate our ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density estimation, we measure the fluorescence count rates for the three hBN samples and find them to be proportional to the estimated ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$ (Fig. 3a inset, and Supplementary Table S1).

We highlight that although the ion implantation dosage spans nearly two orders of magnitude across three hBN samples, the estimated ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density only differs approximately by a factor of 2. This indicates that with larger implantation dosage, one may create more V_B defects, but most of them remain charge neutral^61,62,66,71. Using SRIM (Stopping and Range of Ions in Matter) program, we estimate the created V_B defect density in the experiment to be ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}}^{{{{{{{{\rm{S1}}}}}}}}} \, \approx \, (5.4\pm 0.5)\times 1{0}^{2}$ ppm, ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}}^{{{{{{{{\rm{S2}}}}}}}}} \, \approx \, (2.0\pm 0.2)\times 1{0}^{3}$ ppm and ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}}^{{{{{{{{\rm{S3}}}}}}}}} \, \approx \, (1.8\pm 0.2)\times 1{0}^{4}$ ppm, increasing linearly with the implantation dosage (see Methods). Figure 3b shows the negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ratio, $\eta \equiv {\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}/{\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{}}$, which exhibits a substantial drop with increasing implantation dosages. This suggests that one may need to seek alternative solutions other than simply cranking up the irradiation dosage to achieve higher ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration for future applications in quantum information. We note that if one directly uses ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}}$ from SRIM to represent the negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, the simulated coherent timescales ${T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}$ and ${T}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ will be significantly shorter than the experimental results (see Supplementary Fig. 3b).

Probing the local charged defect environment

The presence of negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble in hBN also leads to a local electric field signal that can be directly probed using the spin degree of freedom of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ (Fig. 4a). Given the mirror symmetry of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ lattice structure respect to the x − y plane, its electric field susceptibility along z vanishes, and one only needs to consider the transverse component of the local electric field. Without any external magnetic field, a transverse electric field to the z-axis of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ (c-axis of hBN), E_⊥, will mix the original $\left|{m}_{s}=\pm 1\right\rangle$ states of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, and split them into two new eigenstates, $\left|\pm \right\rangle$^61,62,72,73. To the leading order, the energy splitting, δ, between $\left|\pm \right\rangle$ is proportional to the strength of the transverse electric field, δ ∝ 2d_⊥E_⊥, where d_⊥ is the ground state transverse electric field susceptibility of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ (Fig. 4a). In reality, the presence of the three first-shell ¹⁴N nuclear spins as well as the intrinsic broadening of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ transitions will lead to additional modification to the measured energy splitting δ, and a detailed discussion of such effect can be found in Methods.

The splitting δ can be probed via the electron spin resonance (ESR) measurement: by sweeping the microwave field frequency and monitoring the fluorescence signals of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, one observes a fluorescence drop when the microwave is resonant with one of the spin transitions. Figure 4b shows the measured ESR spectra for sample S1 and S3 at zero magnetic field. Crucially, we observe that the splitting increases with ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration, consistent with the expectation that a higher charged defect density can generate a stronger local electric field. We also notice a small shift of the ESR center frequencies with increasing implantation dosages, which may originate from the implantation-induced strain effect^25,74,75,76.

To quantitatively understand the density dependence of the measured splitting, we utilize a microscopic model based upon randomly positioned electrical charges inside the hBN lattice. Such model has been successfully applied to capture the measured energy splitting between $\left|{m}_{s}=\pm 1\right\rangle$ sublevels of NV centers in diamond before^61,62. Specifically, we randomly position charged defects surrounding a ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ center at a density ρ_c, and calculate the corresponding transverse electric field E_⊥ at the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ site. Here we assume that these charges consist primarily of the negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ centers themselves (which are electron acceptors) and their associated donors — as a result, the local charged defect density ${\rho }_{{{{{{{{\rm{c}}}}}}}}} \, \approx \, 2{\rho }_{{V}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$. By diagonalizing the lab frame spin Hamiltonian in the absence of an external magnetic field (see Methods), we calculate the transition frequencies of the ESR experiment. The final simulated ESR spectrum is obtained via averaging over different charge defect configurations, as well as the spin states of the three closest hyperfine-coupled ¹⁴N nuclear spins. Since d_⊥ of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ has not been determined before, we survey a range of different d_⊥ in our numerics to obtain a series of simulated ESR splitting at a variety of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density (Fig. 4c). Comparing the experimentally measured ESR splitting δ to the simulated results from our model, we are able to get a rough estimation of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ground state transverse electric field susceptibility, d_⊥ ≈ 40 Hz/(V ⋅ cm⁻¹). We note that the estimated d_⊥ of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ is on the same order of NV center in diamond, ${d}_{\perp }^{{{{{{{{\rm{NV}}}}}}}}}\approx 17\,{{{{{{{\rm{Hz}}}}}}}}/({{{{{{{\rm{V}}}}}}}}\cdot {{{{{{{{\rm{cm}}}}}}}}}^{-1})$⁷⁷.

Outlook

Looking forward, our work opens the door to a number of intriguing directions. First, the characterization and control of coherent dipolar interaction in dense ensembles of spin defects in 2D materials represent the first step to using such platforms for exploring exotic many-body quantum dynamics. One particularly interesting example is to investigate the stability of phenomena such as many-body localization and Floquet thermalization in two and three dimensions. In fact, in long-range interacting systems, the precise criteria for delocalization remain an open question; whereas in Floquet systems, the thermalization dynamics involve a complex interplay between interaction and dimensionality^44,46,49,51. Secondly, the measured low negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ratio at high ion implantation dosage suggests that one may be able to use external electric gating to substantially tune and enhance the portion of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration. Indeed, electric gating has been recently demonstrated as a powerful tool to engineer the charge state of optical spin defects in solid-state materials^78,79,80,81. Finally, the estimated transverse electric field susceptibility highlights the potential use of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ as an embedded electric field sensor for in-situ characterization of heterogeneous materials^62,72,82,83.

Methods

hBN device fabrication

The sample consists of ion-irradiated hBN flakes and a titanium/gold (10/300 nm thick) coplanar waveguide (CPW) with a 50 μm wide central stripline, on a sapphire substrate. hBN flakes were tap-exfoliated from a commercially available hBN single crystal and transferred onto Si substrates. Boron vacancy defects were generated using He⁺ ion implantation with an energy of 3 keV with dose densities, 0.30 ± 0.03 nm⁻² (sample S1), 1.1 ± 0.1 nm⁻² (sample S2), and 10 ± 1 nm⁻² (sample S3), respectively. Error bars on the dosages account for measured current fluctuations during the implantation process. After ion irradiation, hBN flakes were transferred on top of the CPW using the PC/PDMS transfer method^31,53.

Sources of decoherence

The sources of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ decoherence have two major contributions: (1) the Ising coupling to the environmental spin bath, such as nuclear spins and dark electronic spins, and (2) the dipolar interaction within the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble⁵². In particular, the interaction between a ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ center and the local off-resonant spin bath takes the form of Ising coupling (under rotating-wave approximation)

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{Ising}}}}}}}}}=\mathop{\sum}\limits_{k}{A}_{k}^{zz}\,{S}^{z}{{\mathfrak{S}}}_{k}^{z}=\left(\mathop{\sum}\limits_{k}{A}_{k}^{zz}{{\mathfrak{S}}}_{k}^{z}\right){S}^{z},$$

(2)

where ${A}_{k}^{zz}$ represents the strength of the Ising coupling between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ and the kth bath spin, ${{\mathfrak{S}}}_{k}^{z}$ is the spin operator for the kth bath spin, and S^z is the spin-1/2 operator for ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ when restricting to the spin subspace $\left|{m}_{s}=0\right\rangle,\left|{m}_{s}=-1\right\rangle$ (See Supplementary Note 2.1).

A few remarks are in order. First, summing over the Ising coupling to bath spins results in a random on-site field disorder on each ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, $h={\sum }_{k}-\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{k}}{{r}_{k}^{3}}{{\mathfrak{S}}}_{k}^{z}$, which would effectively broaden the electronic spin transition. Second, the strength of the Ising coupling consists of two different contributions: Fermi contact interaction and dipolar interaction. For the closest three spin-1 ¹⁴N nuclear spins, A^zz ≈ 47 MHz is dominated by Fermi contact, which leads to the previously reported seven-peaks hyperfine structure in the measured ESR spectrum³¹. On top of that, the interaction with other far away bath spins like ¹⁰B, ¹¹B and ¹⁴N nuclear spins, as well as dark electronic spins⁵², would further broaden the ESR spectrum. Indeed in our experiment, the measured ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ESR peak associated with the spin transition $\left|{m}_{s}=0\right\rangle \iff \left|{m}_{s}=-1\right\rangle$ (under an external magnetic field ~ 250 G) exhibits a large linewidth which can be well captured using a Gaussian function. The extracted standard deviation of the ESR resonance is around 80 MHz across all three hBN samples studied in this work (Supplementary Fig. 2a), independent of the ion implantation dosages. Last, given the large broadening of the transition, one requires a strong microwave field to efficiently drive the spin state of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. In our experiment, we utilize microwave pulses with a Rabi frequency Ω ≈ 83 MHz (π − pulse duration t_π = 6 ns), similar to the measured ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ linewidth, to achieve rapid spin manipulation with reasonable fidelity (Supplementary Fig. 2b).

The other source for decoherence is the dipolar interaction within the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin ensemble. In the rotating frame, the Hamiltonian that governs the dipolar interaction of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ can be written as

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}=\mathop{\sum}\limits_{i < j}-\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{i,j}}{{r}_{i,j}^{3}}\left({S}_{i}^{z}{S}_{j}^{z}-{S}_{i}^{x}{S}_{j}^{x}-{S}_{i}^{y}{S}_{j}^{y}\right),$$

(3)

where J₀ = 52 MHz nm³, ${{{{{{{{\mathcal{A}}}}}}}}}_{i,j}$ and r_i,j represent the angular dependence and the distance between the i^th and j^th${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spins, and ${S}_{i}^{x}$, ${S}_{i}^{y}$, ${S}_{i}^{z}$ are the spin-1/2 operators for i^th${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ centers. We note that ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}$ corresponds to the energy-conserving terms of the dipolar interaction, i.e., the rotating-wave approximation, when restricting our attention to the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spin subspace $\left|{m}_{s}=0\right\rangle$, $\left|{m}_{s}=-1\right\rangle$ (See Supplementary Note 2.1). For the highest ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density sample S3 in this work, ${\rho }_{{{{{{{{\rm{S3}}}}}}}}}\approx 236\begin{array}{c}+35\\ -31\end{array}$ ppm, we estimate the average dipolar interaction strength between nearby spins to be 〈J〉 ~ 1.2 MHz.

Dynamical decoupling sequences

To selectively isolate the effect of each decoherence source, we introduce three different dynamical decoupling sequences to investigate the coherent dynamics of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble. We start with the most basic sequence, Spin Echo, which applies a single refocusing π-pulse at the center of the time evolution (Supplementary Fig. 3b). This refocusing pulse reverses the on-site field disorder (Ising coupling in Eq. S1) for the second half of the time evolution, ${\sum }_{i}{h}_{i}{S}_{i}^{z}\longrightarrow {\sum }_{i}-{h}_{i}{S}_{i}^{z}$, thus negating the accumulated phase from the static component of the on-site disorder. However, Spin Echo cannot decouple the decoherence that arises from the time-dependent fluctuation of the on-site random field. For the case of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, Spin Echo decay originates from the internal spin flip-flops within the environmental spin bath. As a result, the Spin Echo timescale ${T}_{2}^{{{{{{{{\rm{Echo}}}}}}}}} \, \approx \, 70$ ns offers a good estimation of the correlation time, τ_c, of the spin bath.

Next, the more advanced decoupling sequence XY-8 also employs refocusing π-pulse, but, instead of one, XY-8 applies a series of π-pulses to further isolate the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ from the environment (Supplementary Fig. 3c). The alternating phase of π-pulses along different spin axes (X-Y-X-Y-Y-X-Y-X) is designed to suppress intrinsic pulse errors to higher order. We fix the interval time between π-pulses to be 4 ns, much shorter than the correlation time τ_c of the spin bath, and sweep the total pulse number. As expected, using XY-8 pulse sequence, the local bath spin fluctuation is further decoupled from the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble and we observe significant extensions of the measured coherence times across all three hBN samples.

However, the decoupling mechanics of XY-8 does not apply to the dipolar interaction within the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble, ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}$. The intuition is simple. The microwave π-pulses flip all the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ spins together, and the two-body terms in ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}$ remains unchanged. For instance, under a rotation along the x axis, ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}$ becomes ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}^{{\prime} }={\sum }_{i < j}-\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{i,j}}{{r}_{i,j}^{3}}[(-{S}_{i}^{z})(-{S}_{j}^{z})-{S}_{i}^{x}{S}_{j}^{x}-(-{S}_{i}^{y})(-{S}_{j}^{y})]={{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{dip}}}}}}}}}$.

This is where the specifically designed interaction decoupling sequence, DROID (Disorder-RObust Interaction-Decoupling), comes in⁵⁹. In particular, by applying a series of π/2 rotations along different spin axes to periodically change the frames of interaction (also known as the toggling frame), the leading-order effective Hamiltonian can be described by a simple weighted average of each toggling-frame Hamiltonian (Supplementary Fig. 3d). For spin-1 ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, the resulting effective Hamiltonian, ${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{eff}}}}}}}}}$, takes the form of an isotropic long-range Heisenberg interaction,

$${{{{{{{{\mathcal{H}}}}}}}}}_{{{{{{{{\rm{eff}}}}}}}}}=\mathop{\sum}\limits_{i < j}\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{i,j}}{{r}_{i,j}^{3}}\frac{1}{3}\left({S}_{i}^{z}{S}_{j}^{z}+{S}_{i}^{x}{S}_{j}^{x}+{S}_{i}^{y}{S}_{j}^{y}\right)=\mathop{\sum}\limits_{i < j}\frac{{J}_{0}{{{{{{{{\mathcal{A}}}}}}}}}_{i,j}}{{r}_{i,j}^{3}}\frac{1}{3}{\hat{{{{{{{{\boldsymbol{S}}}}}}}}}}_{i}\cdot {\hat{{{{{{{{\boldsymbol{S}}}}}}}}}}_{j}$$

(4)

Since our initial spin state, ${\otimes }_{i}\frac{{\left|0\right\rangle }_{i}+{\left|-1\right\rangle }_{i}}{\sqrt{2}}$, constitutes an eigenstate of the above Heisenberg Hamiltonian, DROID sequence further suppresses the decoherence effect originating from the dipolar interaction within ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ensemble. This is indeed borne out by our data. The measured coherent timescales under DROID sequence exhibit an additional 2-fold increment compared with XY-8 sequence across all three hBN samples.

Numerical simulation of the coherent dynamics

To quantitatively analyze the coherent timescales measured from the three different samples, we perform numerical simulation using 12 quantum spins randomly positioned onto the hBN lattice with varying density ρ. For the dipolar interaction between ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, we build up the Hamiltonian for each pair of spins using the form in Eqn. (3). As for the on-site random field disorder, we introduce an additional term ${\sum }_{i}{h}_{i}{S}_{i}^{z}$ into the Hamiltonian, with h_i drawn from a Gaussian distribution with standard deviation 80 MHz which is independently characterized from the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ESR resonance spectrum (Supplementary Fig. 2a). To account for the effect of finite pulse duration in the experimental sequence, a microwave driving term $\pm \Omega {\sum }_{i}{S}_{i}^{x(y)}$ is included into the simulation whenever a pulse with specific phase is applied, where Ω = 83 MHz is directly determined using the corresponding Rabi oscillation measurement (Supplementary Fig. 2b). For each pulse sequence, respective π- and $\frac{\pi }{2}$-pulses and intervals are applied as in the experimental procedure, and we also utilize the differential measurement scheme in the simulation to faithfully capture the experimental details. After time evolution, we only use the central spin’s polarization to represent the coherent dynamics, as the far away spins would suffer from significant finite-size effects due to the small system size that one can compute in numerics. At a given ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$, we average for 1000 disorder realizations of random spin positions and on-site fields to obtain the final simulated curve.

To efficiently compute the corresponding quantum dynamics, we employ DYNAMITE, a powerful package providing a simple interface to fast evolution of quantum dynamics⁸⁴. In contrast to the traditional Hamiltonian diagonalization method that requires exponentially increasing time as spin number goes up, DYNAMITE uses the PETSc/SLEPc implementations of Krylov subspace exponentiation and eigensolving, which drastically lower the computational resources for simulations involving large spin number.

${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ Density extraction

Calculation of residuals

In order to extract the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, we compare the experimentally measured XY-8 and DROID coherent timescales to the numerical simulations. In particular, we first fit the coherence time ${{{{{{{{\mathcal{T}}}}}}}}}_{2}$ to ${{{{{{{{\mathcal{T}}}}}}}}}_{2}\propto {\rho }^{-\alpha }$, where ρ is the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density, and m is a free fitting parameter. This fitted line should be a straight line in our log-log plot (Supplementary Fig. 4a), where the slope corresponds to the value of α. We find the fitted value of α ~ 0.8 ± 0.1 for XY-8 to be fairly close to one, and this validates that the main decoherence source for XY8 is from the dipolar interaction whose strength scales linearly with ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density. Note that ${{{{{{{{\mathcal{T}}}}}}}}}_{2}^{{{{{{{{\rm{D}}}}}}}}}$ theoretically should be unchanged with increasing ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ density since the DROID sequence should decouple the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ dipolar interaction, but due to other factors like imperfect drive the decoherence timescale still decrease. To extract the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities, we calculate the sum of squared relative residuals of XY-8 and DROID between the experimental values and the fitted ${{{{{{{\mathcal{T}}}}}}}}$.

$${{{{{{{\rm{Residuals}}}}}}}}={\left(\frac{\log {{{{{{{{\mathcal{T}}}}}}}}}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}-\log {T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}}{\log {T}_{2}^{{{{{{{{\rm{XY8}}}}}}}}}}\right)}^{2}+{\left(\frac{\log {{{{{{{{\mathcal{T}}}}}}}}}_{2}^{{{{{{{{\rm{D}}}}}}}}}-\log {T}_{2}^{{{{{{{{\rm{D}}}}}}}}}}{\log {T}_{2}^{{{{{{{{\rm{D}}}}}}}}}}\right)}^{2}$$

(5)

Here, we adopt the minimums of the residual curves as our estimated density and estimate their errors from the range values whose residuals lie within 5% of the minimum (Supplementary Fig. 4b blue shaded regions). This put the estimated densities for the three samples to be $123\begin{array}{c}+8\\ -8\end{array}$, $149\begin{array}{c}+25\\ -21\end{array}$, and $236\begin{array}{c}+35\\ -31\end{array}$ ppm respectively. We note the discrepancy between the measured and simulated timescales may stem from (1) pulse imperfections in the experiment, and (2) finite-size effect, i.e., how many spins one can simulate, in numerical simulations. The presence of nearby bath spins (nuclear spins and dark electronic spins) leads to a large broadening of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ transition with measured resonance linewidth σ ~ 80 MHz (standard deviation). Therefore, with Rabi driving frequency Ω ≈ 83 MHz, which is comparable to the broadening σ, the applied pulses will have significant imperfections. To better capture this effect in simulation, we incorporate a random static on-site field on each spin drawn from the measured resonance distribution, and average across all different configurations. However, the static on-site field does not portray the whole picture of the real experiment, where the bath spins also have dynamics. Due to the limited system size we can compute, it is particularly difficult to simulation this effect.

The measured T₂ decay in the experiment is a combination of both decoherence and imperfect pulse effect. Therefore, for DROID sequence where one decouples a large portion of the interactions, such pulse imperfection becomes more evident compared to XY8 sequence where the T₂ decay is governed by the dipolar interactions within ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. This also explains why the discrepancy between experiment and simulation is larger in DROID T₂ compared to XY-8 T₂.

Fluorescence counts

To verify our estimation of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities, we also measure the fluorescence counts N_tot of the three samples with excitation laser of the same power (~0.4 mW). The resulting counts are (0.64 ± 0.03) × 10⁶, (0.84 ± 0.03) × 10⁶, and (1.44 ± 0.04) × 10⁶ photons/s respectively for sample S1, S2, and S3 (see Supplementary Table 1). It is important to note that this total counts N_tot also include the noise background.

To account for this non-negligible background, we define the adjusted counts by

$${N}_{A}={N}_{{{{{{{{\rm{tot}}}}}}}}}\times C,$$

(6)

where C is the measured ESR contrast for the three samples. Here we also normalize N_A to 1 for sample S1 for better comparison, and the resulting adjusted counts across three samples do closely follow the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ densities estimated from spin coherent dynamics (see Main Text Fig. 3a inset).

SRIM simulation

To estimate the total created boron-vacancy V_B defect density at different ion implantation dosages, we perform a detailed calculation with full damage cascades using Stopping and Range of Ions in Matter (SRIM)⁶³. In particular, we choose the incident beam to be Helium ion with the energy of 3 keV, targeting hBN layer with a thickness 100 nm. The simulated damage events distribution is shown in Supplementary Fig. 5a. We find that the created vacancies are distributed within the initial ~ 60 nm of the hBN sample. By integrating target vacancies at different depths, we conclude that each ion on average creates ~ 11 boron vacancies. For sample S1, S2, S3 with ion implantation dosages of 0.30 ± 0.03 nm⁻², 1.1 ± 0.1 nm⁻², and 10 ± 1 nm⁻², the created area density of boron vacancy is around 3.3 ± 0.3, 12.1 ± 1.2, and 110  ± 10 nm⁻² respectively. Using the atomic number density for hBN, 101.9 atom nm⁻³, we estimate the total created boron vacancy density in sample S1, S2 and S3 to be ~ (5.4 ± 0.5) × 10², (2.0 ± 0.2) × 10³ and (1.8 ± 0.2) × 10⁴ ppm.

In main text Fig. 3, by comparing the experimentally measured coherent timescales to our numerical simulations, we conclude that only a small portion of the created boron-vacancy centers (V_B) from ion implantation process are in the desired negatively charged state (${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$). If one naively uses the V_B density estimated from SRIM to represent ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration, we find that the simulated timescales are much shorter than the experimental results for all three hBN samples investigated in this work (Supplementary Figure 5), indicating significant overestimation of the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ concentration. Note that when the electronic spin density exceeds ~1000 ppm in our simulation, due to the strong dipolar interaction within the spin system, our finite-duration microwave pulses cannot faithfully drive the spin anymore. As a result, the simulation exhibits nearly vanishing spin polarization even after the first set of driving pulses. Here we use a dashed line to extend the DROID simulation above ~ 1000 ppm for better comparison.

Numerical simulation of ESR splittings and extraction of transverse electric field susceptibility

In this section, motivated by pioneer studies of nitrogen-vacancy centers in diamond^61,62,73, we model the measured ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ ESR splitting in the absence of any external magnetic field using a microscopic model based upon the local charged defects in hBN. In particular, the presence of randomly distributed charges will create a vector electric field $\vec{E}$ at the site of each ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$. Intuitively, we expect these charges to consist primarily of the negatively charged ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ centers themselves (which are electron acceptors) and their associated donors — as a result, the local charged defect density ${\rho }_{{{{{{{{\rm{c}}}}}}}}} \, \approx \, 2{\rho }_{{V}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$ with ${\rho }_{{V}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$ the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ defect density. The corresponding ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ electronic spin-1 ground state Hamiltonian can be written as:

$${{{{{{{\mathcal{H}}}}}}}}={D}_{{{{{{{{\rm{gs}}}}}}}}}{S}_{z}^{2}+{\Pi }_{x}\left({S}_{y}^{2}-{S}_{x}^{2}\right)+{\Pi }_{y}({S}_{x}{S}_{y}+{S}_{y}{S}_{x})+\mathop{\sum }\limits_{i=1}^{3}{A}_{zz}{I}_{z}^{i}{S}_{z},$$

(7)

where D_gs ≈ 3.48 GHz is the zero-field splitting between $\left|{m}_{s}=0\right\rangle$ and $\left|{m}_{s}=\pm 1\right\rangle$ spin levels, Π_{x, y} = d_⊥E_{x, y} characterizes the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ coupling to local perpendicular electric field with susceptibility d_⊥, S and I ⁱ are the spin-1 operators for ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ electronic spin and the closest three ¹⁴N nuclear spins respectively, and A_zz ≈ 47 MHz is the hyperfine coupling strength^31,52.

A few remarks are in order. First, the z is defined along the c-axis of hBN (perpendicular to the lattice plane), x and y lie in the lattice plane, with x oriented along one of the three vacancy-nitrogen bonds. Due to the mirror symmetry of ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ with respect to the x − y plane, its electric field susceptibility in z vanishes, d_∥ = 0⁸⁵. Second, the presence of ¹⁴N nuclear spins generate an effective local magnetic field along z direction on the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$, whose strength depends on the specific nuclear spin states. Since the three nuclear spins are all in fully mixed states, we can simply set each ${I}_{z}^{i}=\{-1,0,+1\}$ with equal probabilities, and average across different nuclear spin configurations to account for their contributions.

To quantitatively capture the experimentally measured ESR spectrum at zero magnetic field, we randomly sample the charged defects in space at a given density ${\rho }_{{{{{{{{\rm{c}}}}}}}}} \, \approx \, 2{\rho }_{{V}_{{{{{{{{\rm{B}}}}}}}}}^{-}}$ and calculate the perpendicular electric field, E_{x, y}, at the ${{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}$ site from the closest 10 sampled charges. By diagonalizing the Hamiltonian in Eqn. (7), we obtain the corresponding eigenstates and the associated eigenenergies, with which we can back out the resonances in ESR. To account for the intrinsic broadening of the spin resonance, we convolve each resonance with a Gaussian profile with standard deviation σ ≈ 25 MHz to better capture the experimental data. The final simulated ESR spectra are obtained by averaging over 1000 different randomly positioned charge configurations.

Supplementary Fig. 6a shows the simulated ESR spectra on top of the experimental data using ${\rho }_{{{{{{{{{\rm{V}}}}}}}}}_{{{{{{{{\rm{B}}}}}}}}}^{-}}=\{123\begin{array}{c}+8\\ -8\end{array},149\begin{array}{c}+25\\ -21\end{array},236\begin{array}{c}+35\\ -31\end{array}\}\,{{{{{{{\rm{ppm}}}}}}}}$ for sample S1, S2, and S3, and d_⊥ = 40 Hz/(V ⋅ cm⁻¹). The excellent agreement between the numerical simulation and experiment highlights the validity of our microscopic charge model. We remark that, given the large intrinsic broadening of the spin transitions, the measured ESR spectra can be also well captured with a sum of two Lorentzian profiles, as evinced in Supplementary Fig. 6b, with which we extract the values of the ESR splitting δ in the main text Fig. 4c.