Decoherence of nitrogen-vacancy spin ensembles in a nitrogen electron-nuclear spin bath in diamond

Nitrogen-vacancy (NV) centers in diamond have been developed into essential hardware units for a wide range of solid-state-based quantum technology applications. While such applications require the long spin lifetimes of the NV centers, they are often limited due to decoherence. In this study, we theoretically investigate the decoherence of NV-spin ensembles induced by nitrogen impurities (P1 centers), which are one of the most dominant and inevitable magnetic field noise sources in diamond. We combined cluster correlation expansion and density functional theory to compute the Hahn-echo spin-coherence time of the NV centers for a broad range of P1 concentrations. Results indicate a clear linear dependence of T2 on P1 concentrations on a log scale with a slope of −1.06, which is in excellent agreement with previous experimental results. The interplay between the Jahn–Teller effect and the hyperfine interaction in the P1 center plays a critical role in determining the bath dynamics and the resulting NV decoherence. Our results provide a theoretical upper bound for the NV-spin T2 over a wide range of P1 densities, serving as a key reference for materials optimization and spin bath characterization to develop highly coherent NV-based devices for quantum information technology.


INTRODUCTION
Nitrogen-vacancy (NV) centers in diamond are one of the leading solid-state-based qubit platforms, enabling diverse quantum applications including quantum information processing [1][2][3][4] , quantum sensing [5][6][7][8] , and quantum networks 9,10 . For a spin qubit, a long spin lifetime is one of the key properties that an NV center should possess because it plays a critical role in determining the performance of a spin qubit in quantum applications 6 . Combined with spin bath decoupling methods or bath engineering techniques, the Hahn-echo spin-coherence time (T 2 ) of an NV center has been shown to reach~1 s [11][12][13] . In practical applications, however, the spin lifetime of an NV center could be significantly degraded due to intensified decoherence caused by paramagnetic spin defects 14 , which in turn could limit the range of NV-spin qubit-based quantum applications 14,15 .
NV-spin decoherence has been extensively investigated and the dominant mechanism has been found to involve the coupling of the NV-spin dipolar to the fluctuating spin bath in diamond 16 . Two dominant types of spin baths identified are those with 13 C nuclear spins (natural abundance 1.1%, I = 1/2) and those considered as electronic spin baths with substitutional nitrogen electron-nuclear spins (the so-called P1 bath, 14 N 99.6%, S = 1/2, and I = 1). In a 1.1% 13 C nuclear-spin-dominant bath, the NV-spin's Hahn-echo coherence time (T 2 ) has been reported to be larger than 600 μs [17][18][19] and its decoherence dynamics is theoretically well understood 20 . Building upon this established knowledge, several methods have also been developed to suppress the 13 C-induced decoherence of the NV spin, including a dynamical decoupling method 21 and 12 C isotopic purification of a diamond sample 18 . Meanwhile, it is well known that the paramagnetic spins in diamond, or P1 centers 22 , are one of the dominant sources of NVspin decoherence [23][24][25][26][27] . It has been shown that the NV-spin lifetime is much more fragile under the presence of P1 spins than in a 13 C nuclear-spin bath owing to the larger gyromagnetic ratio of P1 spins 24,25,27,28 . Nevertheless, P1-induced NV-spin decoherence is less understood than 13 C-induced decoherence.
Recently, understanding the NV-spin decoherence mechanism due to various P1 concentrations ([P1]) has gained much interest [29][30][31][32][33][34] in several important applications such as nuclearspin hyperpolarization 35 , quantum sensing 5,6 , and quantum information applications 1,36,37 . Such applications may require either high NV densities or the controlled creation of P1 centers. In the NV center generation process, however, P1 centers are inevitably created with different concentrations (1-100 ppm) owing to the low N-to-NV center yield rate (<25%) 35,38 , which, in turn, may lead to the degraded coherence of the NV center due to the P1 centers 7 . Thus, it is of paramount importance to systematically understand P1-induced NV decoherence as a function of [P1].
Recent studies have made significant progress towards the quantitative analysis of P1-induced NV-spin decoherence. In particular, it has become possible to determine [P1] in a diamond sample with high precision 30,31 , enabling measurements of the NV's T 2 as a function of [P1] [29][30][31][32] . Recently, Bauch et al. 29 reported the T 2 of an NV-spin ensemble as a function of [P1], showing a linear dependence of T 2 on [P1] on a log scale with a stretched exponential parameter of~1.37. In another recent work by Li et al. 30 a local spin bath measurement technique is developed to extract the position-dependent concentration of various paramagnetic defects in a diamond sample. Importantly, however, the linear dependence between T 2 and [P1] found in the semiclassical simulation of the previous study 29 was much less clear. In addition, the stretching exponent measured varied significantly depending on [P1], which also deviates from that of the previous result 29 . Li et al. 30 and other researchers have suggested that the deviation may be caused by extra paramagnetic defects other than the P1 centers in the sample.
The lack of a full understanding of the NV's T 2 vs.
[P1] and the role of extra electronic spins can benefit from the accurate theoretical prediction of T 2 ([P1]), which can provide an upper bound for T 2 and a reference for the stretched exponential parameter of a P1dominant bath. Early theoretical studies have shown that the basic property of P1-induced NV-spin decoherence can be analyzed by using semi-classical theories based on the Ornstein-Uhlenbeck (O-U) process. Notably, a similar semi-classical method was used to understand recent experimental measurements 29,30 , and the overall qualitative tendency of T 2 as a function of [P1] was reproduced. The semi-classical results, however, underestimate the T 2 and overestimate the slope over a wide range of P1 concentrations 29 . In this regard, first-principles calculations with a predictive power are necessary to provide a quantitative theoretical guidance to P1induced NV decoherence.
In this study, we perform fully quantum mechanical calculations by combining cluster correlation expansion and density functional theory (DFT) to investigate the P1-induced decoherence dynamics of a negatively charged NV ensemble for a broad range of [P1]. We compute the Hahn-echo T 2 as a function of [P1] and our results provide a theoretical upper bound for T 2 from 1 to 100 ppm. We show that the T 2 ([P1]) exhibits a clear linear dependence on a log scale with a slope of −1.06, which is in excellent agreement with the experimental result of −1.07 [29][30][31] . We found that the stretched exponential parameter is close to 1 for the [P1] from 1 to 100 ppm. Building upon these results, we show that the Jahn-Teller-derived anisotropy in the hyperfine interaction of the P1 center plays a crucial role in determining P1 bath dynamics. Our study offers a complete understanding of the P1-induced spin decoherence dynamics of an NV ensemble, which can provide a key reference for developing diamond samples with optimized NV-spin lifetimes for different applications.

System and model
We adopt a central spin model to compute the decoherence dynamics of an NV-spin interacting with a substitutional nitrogen impurity spin bath. We use the Hahn-echo pulse sequence, which has a π pulse in between two free evolution time periods τ to compute the homogeneous dephasing time T 2 of NV-spin ensembles 39,40 . The degree of NV-spin decoherence is obtained by computing the off-diagonal element of the reduced density matrix of the NV spin after tracing out the bath degrees of freedom at the end of the 2τ free evolution time 39 (see Methods for details). Figure 1 shows our central spin model in detail. The NV center is a spin-1 electronic defect, which consists of a C vacancy and substitutional N atom paired along the [111] crystal direction. An external magnetic field of 500 G is applied in the same direction as the symmetry axis of the NV spin. The P1 center is an isolated substitutional nitrogen impurity that replaces a carbon atom in the diamond lattice. In our model, P1 spins are randomly distributed around the NV spins, and we consider various [P1] from 1 to 100 ppm. An NV spin is coupled to the P1 centers via the magnetic dipole-dipole interaction (see Supplementary Note 1). The P1 center accompanies an electron spin-1/2 and a nuclear spins-1 (spin-1/2) for 14 N ( 15 N, natural abundance of 0.37%), which are strongly coupled to each other via a hyperfine interaction (see Supplementary Note 2). The electronic spin of the P1 center is highly localized at the lattice site due to the Jahn-Teller (JT) effect 41 .
The JT effect creates strong anisotropy in the hyperfine interaction of the P1 center 42,43 and induces a strong lattice relaxation in such a way that the N atom is distorted away from one of the four nearest neighboring C atoms, leading to four possible hyperfine anisotropy axes. In Supplementary Information, we provide a detailed description of the anisotropic hyperfine tensors. We compared our DFT calculations with experimental results 37 , the comparison of which demonstrated excellent agreement for the four different hyperfine anisotropy axes (see Supplementary Note 2). In our central spin model (see Fig. 1), the four possible anisotropy axes are randomly assigned to the P1 centers. In addition, the given anisotropy axis for the P1 centers is fixed during the simulation as the time scale of the JT axis change is known to be from 10 3 to 10 5 s in low temperatures (≤200 K) 44 , which is much larger than the time scale of NV-spin decoherence, which is in the order of hundreds of microseconds.
To compute the NV-spin's T 2 induced by a large number of P1 centers, we employ the cluster correlation expansion method (CCE) 39,40,45 , which enables a predictive quantum mechanical computation of T 2 without assuming any adjustable theoretical parameters. The CCE method was successfully applied to a wide range of central spin models including both nuclear 20,46,47 and electron spin baths 48,49 , yielding an excellent agreement with experimental results 47,49 . Further details on the spin Hamiltonian, theoretical methods, and the numerical convergence tests can be found in the Methods section and Supplementary Notes 3 and 4.  consistent with recent experimental measurements, which yielded n ¼ 1:37 ± 0:23 29 . It is worth noting that O-U based semiclassical theory predictions of n could vary from 1 (motionalnarrowing regime) to 3 (quasi-static regime) 29,33 . Several early experimental studies reported scattered results for n from 1 to 3 26,31,33 . For the NV-spin ensemble in a pure P1 bath, our CCE calculations show that the shape of the NV's spin-echo decay curve is close to the curve of a single exponential function.
The large deviation of the experimental T 2 results from the theoretical calculations obtained for a pure P1 bath indicates the important potential role of parasitic electron spins in addition to P1 centers that act as extra decoherence sources 30,34,50 . In Fig. 3, the effect of adding extra electronic spins (S = 1/2 with g = 2) in a P1 bath is examined and results indicate that the experimental T 2 time can be reproduced in theory when the concentration of the extra spins ([Ex]) is similar to [P1]. Figure 3a, Fig. 3c, d), the computed T 2 for [P1] = 1 ppm is 60.7 μs, which matches well with the experimental T 2 time of 65 μs. n is computed to be~0.74 and~0.92 for [Ex] = 0 ppm and 2 ppm, respectively. For a P1 spin bath with [P1] = 5 ppm (see Fig. 3e, f), the experimentally measured T 2 time of~30 μs is theoretically reproduced when [Ex] is 5 ppm, yielding T 2 = 28.  Figure 4 compares the NV-spin coherence computed for a pure P1 bath, whose JT axes are random, with two hypothetical P1 bath models for [P1] = 1, 10, and 100 ppm: the P1 centers with nohyperfine interaction and the P1 centers whose JT axes are all fixed in one direction. The P1 spin bath with no-hyperfine interaction is effectively the same as a bare electron spin bath in the time scale of a few hundreds of microseconds because the slow N nuclear-spin bath is decoupled from the electron spins. As shown in Fig. 4, the NV coherence in such "no-hyperfine" P1 bath model quickly decays due to the flip-flop transitions between the Fig. 2 Coherence times of diamond NV ensembles in a P1 bath. Hahn-echo coherence times (T 2 ) as a function of [P1] with random P1 JT axes (red dots) compared with previous experimental data of T 2 (black 29 and brown 30,31 crosses). An external magnetic field of 500 G was applied along the [111] direction for the computation. Two hypothetical P1 bath models were considered: a P1 bath in which all the JT axes are polarized in the [111] direction (blue dots) and a P1 bath in which the P1 hyperfine interaction is ignored (green dots). The gray shaded line is the fitted line for the experimental data to be compared with the red shaded line, which is the fitted line for our computational results. Inset: stretched exponential parameters (n) for the P1 bath models with nohyperfine interaction (green marker), polarized P1 JT axes (blue marker), and randomized P1 JT axes (red marker). bare electron spins, yielding T 2 times of 91.15 μs, 8.85 μs, and 0.87 μs for [P1] = 1, 10, and 100 ppm, respectively. Interestingly, however, we observed that the NV's coherence is significantly enhanced when we include the hyperfine interaction between the electron and nuclear spins of the P1 center. For [P1] = 1, 10, and 100 ppm, the T 2 times are increased to 249.62 μs, 24.83 μs, and 2.47 μs, respectively, when the JT axes are fixed to one direction. The T 2 times are further increased to 406.61 μs, 38.73 μs, and 3.32 μs for [P1] = 1, 10, and 100 ppm, respectively, when the JT axes in the P1 bath are randomized. In addition to the P1 concentrations considered in Fig. 4, we also compared the T 2 times computed using the three bath models for other [P1] from 1 ppm to 100 ppm, the results of which are reported in Fig. 2 and from which the same conclusions are drawn. Our calculation results show that the electron spin flip-flop dynamics are significantly suppressed in the presence of the N nuclear spin due to the hyperfine interaction and JT-induced anisotropy in the hyperfine interaction.
To elucidate the hyperfine-induced suppression of the P1 electron spin flip-flop dynamics, we analyzed the energy levels of a P1-P1 spin pair for the three P1 bath models in Fig. 5. Figure 5a shows the energy levels in the absence of the hyperfine interaction. We only considered the states whose electron spins are up and down ( "#i j or #"i j ) because they are coupled by the electron spin dipolar flip-flop interaction. Considering the nuclearspin state ( þi j ; 0i j ; Ài j ), there are 18 possible states whose electron spins are up and down as listed in Fig. 5b. In the absence of the hyperfine interaction, there are 9 active flip-flop transition channels as shown in Fig. 5a due to the nuclear Zeeman energy splitting (Δ n $ 154 kHz) and the mismatch in the nuclear-spin states (i.e., " þ# 0i j $ # þ" 0i j is allowed, " À# Ài j $ # 0 " Ài j is not allowed). However, when the hyperfine interaction is turned on in the bath, the number of these flip-flop channels is greatly reduced. For two P1 centers having the same JT axes, the number of flip-flop channels is three due to the energy level mismatch induced by the hyperfine interaction (Δ HF % 1 2 A on jj % 57 MHz). In addition, for a P1-P1 pair with different JT axes, the degenerated states obtain an additional energy shift (δ HF % A k À A 0 k % 29 MHz), which reduces the  ; 9i j for the " 1 # 2 i j -derived states, and 1 0 i j ; ; 9 0 i j for the # 1 " 2 i j -derived states as specified in (b). In the absence of the hyperfine interaction (a), the spin states are split by nuclear Zeeman splitting (Δ n ) and 9 electron spin flip-flop transition channels are allowed (indicated by blue arrows). c Energy levels of a P1 pair having the same JT axes parallel to the [111] direction. The energy levels are largely shifted by the hyperfine interaction on the order of tens of MHz (Δ HF ). The electron spin flip-flop transitions are allowed only between 1i $ j j1 0 i; 5i $ j j5 0 i and 9i j $ 9 0 i j due to the energy level shifts. d Similar energy level diagram for a P1 pair whose JT axes are not parallel (see inset). The energy levels are further split by the difference in the hyperfine interactions (δ HF ). The electron spin flip-flop transition is only allowed between 5i j $ 5 0 i j . Additionally, the computed T 2 times are consistently larger than previous experimental results, which indicates the importance of other decoherence sources in describing the NV-spin decoherence in a P1 spin bath. After implementing extra electron spins as an additional source of decoherence, we found that the calculated T 2 time becomes similar to the experimental T 2 time when the concentration of the extra spins is similar to the P1 concentration. It should be noted, however, that other paramagnetic defects or complexes may be generated in different diamond sample preparation conditions, acting as extra parasitic spins in the diamond lattice. Further systematic investigation on other paramagnetic spins will help to better understand NV-spin decoherence in practical applications.
Electron dipolar spin ensembles in diamond are already being exploited for practical applications such as quantum sensing or quantum simulations. Understanding the bath spin ensemble dynamics derived from the results obtained in our study could be a vital resource in the engineering of NV-spin-based scalable quantum systems for quantum-enhanced technologies.

Cluster correlation expansion computation
To compute the T 2 time of diamond NV ensembles, we implemented the single-sample CCE theory in our in-house CCE code, and several important features of our CCE computation can be summarized as follows. The P1 bath was clusterized in terms of the number of P1 centers instead of the number of spins in the original CCE formulation. For example, in our CCE-2 method, two P1 centers are considered in a cluster. In this way, we treat both the electron spin and the 14 N nuclear spin on an equal footing and fully consider the anisotropy of the hyperfine tensor for individual JT axes. In addition, we employed the so-called single-sample approach, which treats each bath state independently 39 . We found that our single-sample CCE calculations give a converged T 2 time when averaged over 100 bath states. For an ensemble average, we used 100 different bath spin configurations. We note that in the single-sample CCE method, one can take into account the interaction between clustered bath spins and external spins outside the cluster by adding their Ising-like interaction as a mean-field effect. The single-sample CCE method combined with the mean-field method was shown to produce a robust solution for sparse electron spin baths at the CCE-2 level of theory 48,49 , which we adopted in this study. In Supplementary Note 2, details of our CCE computation and convergence test results are provided.

Density functional theory calculations
To compute the spin Hamiltonian parameters of a P1 center, we conducted DFT calculations using a plane-wave basis set with an energy cutoff of 85 Ry along with projector-augmented wave (PAW) pseudopotentials (QE PAW v0.3.1 set) 51,52 as implemented in the QUANTUM ESPRESSO code 53 . We employed the Perdew-Burke-Ernzerhof semi-local functional to describe the exchange-correlation potential 54,55 . To simulate the presence of an isolated P1 center in diamond, we used the supercell method. The supercell is oriented in the [111] direction and it contains 1007 atoms. The Brillouin zone is sampled with the Γ point only. We used the GIPAW module of the QE code to compute the hyperfine tensor and the quadrupole moment of a P1 center. The core polarization effects were included in the evaluation of spin densities near the nuclei. Our DFT calculations were in excellent agreement with the experimental results for the hyperfine tensor and the quadrupole moment as described in detail in Supplementary Note 2.

DATA AVAILABILITY
The data that support the findings of this study are available upon reasonable request to the corresponding author.