Non-flipping 13C spins near an NV center in diamond: hyperfine and spatial characteristics by density functional theory simulation of the C510[NV]H252 cluster

Single NV centers in diamond coupled by hyperfine interaction (hfi) to neighboring 13C nuclear spins are now widely used in emerging quantum technologies as elements of quantum memory adjusted to a nitrogen-vacancy (NV) center electron spin qubit. For nuclear spins with low flip-flop rate, single shot readout was demonstrated under ambient conditions. Here we report on a systematic search for such stable NV−13C systems using density functional theory to simulate the hfi and spatial characteristics of all possible NV−13C complexes in the H-terminated cluster C510[NV]-H252 hosting the NV center. Along with the expected stable ‘NV-axial−13C’ systems wherein the 13C nuclear spin is located on the NV axis, we found for the first time new families of positions for the 13C nuclear spin exhibiting negligible hfi-induced flipping rates due to near-symmetric local spin density distribution. Spatially, these positions are located in the diamond bilayer passing through the vacancy of the NV center and being perpendicular to the NV axis. Analysis of available publications showed that, apparently, some of the predicted non-axial near-stable NV−13C systems have already been observed experimentally. A special experiment performed on one of these systems confirmed the prediction made.


Introduction
Hybrid spin systems consisting of the electronic spin (e-spin) S=1 of single nitrogen-vacancy (NV) centers in diamond coupled to the intrinsic nuclear spin (n-spin) of its own nitrogen atom and, potentially, to nearby n-spins of isotopic 13 C atoms presenting in the diamond lattice, are now widely used to implement roomtemperature quantum technologies for quantum information processing, sensing and metrology (see, e.g., recent reviews [1][2][3][4][5]). In these systems, the 14 N/ 15 N or 13 C n-spins, with their excellent coherence times, serve as quantum memories accessed via the more easily-controllable e-spin of the NV center that possesses the property of optical pumping and e-spin-projection-dependent fluorescence, allowing its initialization and readout even at room temperature. Currently, the techniques for creating a given spin state of e−n 14 NV/ 15 NV or NV− 13 C spin complexes as well as coherent manipulation of their states to implement one-and two-qubit gates are well established [2,3]. An essential prerequisite for high-fidelity spin manipulation with tailored-control pulse sequences is a complete knowledge of the hyperfine interactions (hfi) in such spin systems, which split the NV e-spin sublevels m S =±1 and can induce stochastic n-spin flipping that violate its coherence. The absence of the last process is of great importance for many quantum-technology applications of e−n spin systems in diamond benefiting from the long coherence time of n-spins. In particular, this is the case for the 14 NV/ 15 NV systems because in these the quantization axes for the NV e-spin and the 14 N/ 15 N n-spins coincide, the respective hfi matrices are intrinsically diagonal, and hence these n-spins are not subjected to hfi-induced flip-flops.
In turn, the I = 1/2 spins of 13 C atoms, presenting in the diamond lattice with the probability 1.1% under natural conditions and with much smaller probability in isotopically engineered diamond samples [6], allow an increase in the number of 'working' n-spins. This is important for many NV applications, in particular, for implementation of small multi-qubit quantum registers [1][2][3] and quantum memories [7], or for quantum error correction [8] which protects quantum states by encoding a logical qubit in multiple physical qubits. Distant 13 C n-spins located somewhat far from the NV center are especially suitable for these purposes since they do not dephase under the NV center spin. Such weakly-coupled n-spins can be detected and characterized experimentally [9][10][11][12][13] using dynamical decoupling methods [14][15][16] that, in particular, allow direct probing of their stochastic flip-flop dynamics resulting from hfi with the NV center [17].
Among various hfi-coupled NV− 13 C spin systems with different positions of the 13 C n-spin in the diamond lattice, there are several positions located at the NV axis [9,[18][19][20] wherein, by analogy with the 14 N/ 15 N e−n spin systems, the 13 C n-spin does not undergo hfi-induced flip-flops due to the coincidence of the e-and n-spins' quantization axes. The nearest-to-the-vacancy 'axial' 13 C position for such an 'NV-axial− 13 C' system exhibiting strongest hfi with the NV center is disposed at the distance of about 6.5 Å from the N atom on the vacancy side at the NV axis [18]. In [18], for this specific 'NV-axial− 13 C' spin system, the characteristic zero-field hfi-induced splitting = 187 kHz of the substates m S = ±1, was predicted by density functional theory (DFT) simulation of the H-terminated cluster C 291 [NV] -H 172 hosting the NV center. More recently [19,20], analogous DFT simulation of hfi characteristics has been done for the larger C 510 [NV] -H 252 clusters, and the spatial and hfi characteristics for eight more distant 'NV-axial− 13 C' systems were found [19]. In any case, the number of available 'axial' 13 C positions is insignificant and the question arises of whether there are other stability (or nearstability) positions for 13 C spins near the NV center in diamond. Considering the great interest in experimental research into non-flipping e−n spin systems, here we did a systematic study of all possible NV− 13 C systems in the C 510 [NV] -H 252 cluster, to find among them the stable (or near-stable) n-spins. We showed that, along with the above-mentioned stable 'NV-axial− 13 C' systems, there are many other near-stable non-axial NV− 13 C spin systems wherein the 13 C n-spins are located basically in the diamond bilayer being perpendicular to the NV axis and containing the vacancy of the NV center.

C system
We start with an analysis of the internal dynamics of an arbitrary NV− 13 C spin system comprising the e-spin  S =1 of the ground-state NV center, coupled to the n-spin  I =1/2 of a 13 C atom disposed somewhere in a diamond lattice, and consider the spin Hamiltonian H of the system written in the principal axes coordinate system of the NV center (NV-PACS), where the Z axis is taken along the С 3V symmetry axis of the center while the X and Y axes are fixed arbitrarily. The spin operators  S and  I have the components S X , S Y , S Z , I X , I Y , I Z , being the standard S=1 and I=1/2 spin matrices in the Cartesian representation. The first term in (1) takes into account the zero-field splitting of the m S = 0 and m S = ±1 e-spin substates due to dipole−dipole interaction between two unpaired electrons of the center (D = 2870 MHz [21]) underlying the triplet e-spin  S = 1 of the center. The second and third terms describe the Zeeman interactions of the two spins  S and  I with the external magnetic field  B.γ e =2.803 MHz gauss −1 and γ n =1.071 kHz gauss −1 are the NV center e-spin and the 13 C n-spin gyromagnetic ratios, respectively. The last term in (1) describes the hfi of the two spins  S ,  I and can be written as A is the hfi tensor presented in the NV-PACS by the symmetric matrix A KL (K, L=X, Y, Z), whose elements depend on the position of the 13 C atom with respect to the NV center. As we are focusing here on NV− 13 C systems, for simplicity we do not take into account in (1) the presence of the internal n-spin I ( N) = 1 of the 14 N nucleus of the center and its intrinsic hfi with the NV e-spin as well as the quadrupole moment Q = −5.01 MHz [21] of the 14 N nucleus. (Note in this regard that, experimentally, the 14 N n-spin projection of the single NV center can be kept fixed, while the effect of the quadrupole moment results only in the total shift of the energy levels of the NV− 13 С spin system under consideration.) As is well known [22], the hfi tensor  A can be subdivided into the isotropic (Fermi contact) and the anisotropic (dipolar) parts: e n  2   2  2   5  e n  2  5 where I 3 is the unit matrix 3×3,  R is the vector connecting the electron and the nucleus, having spatial components R K(L) (K, L=X, Y, Z) and length R. Y  ( ) R is the electron wavefunction at the position of the 13 C nucleus, and the angle brackets are used to indicate the average over electron spatial coordinates. Evidently, the isotropic part can be determined from the full hfi matrix A KL as A iso =Tr(A)/3, and one can check that Tr(T)=0.
In the general case the matrix A KL of the tensor (2a) is non-diagonal in the NV-PACS due to the presence of non-diagonal elements in the dipole−dipole term T KL at K ≠ L. However, symmetry properties allow the zeroing out of some of them. In particular, for each NV− 13 C spin system the T KL matrix can be simplified by choosing the X axis in such a way that the XZ plane passes through the 13 C atom of the system [23]. In this specific NV− 13 C-PACS, all elements of the respective hfi matrix ¢ A KL being non-invariant with respect to the inversion of the Y coordinate are zero due to the C S symmetry of the system. The hfi term then takes the form [23]: . q U is the rotation matrix about the Z axis by some angle θ, which can be determined straightforwardly when one knows the coordinates of the particular 13 C atom.
In many practical cases (excluding those at B∼1027 gauss, where avoided-crossing of sublevels with m S =0 and m S =−1 takes place), a good approximation for the energies, eigenstates and transition rates for the NV− 13 C spin system is provided by the secular approximation where only the terms with S Z are kept in the hfi term of the Hamiltonian (1): ZZ Z Using this approximation, one can find that the energy levels of an arbitrary spin system NV− 13 C in a magnetic field B, aligned along the Z axis, are  Additionally, according to equation (3c), the magnetic field B||OZ splits the e-spin substate m S =0 by the value d g = B n n due to the nuclear Zeeman effect. Respectively, in the secular approximation, the eigenfunctions Y ñ a | , corresponding to the eigenvalues (3), can be written explicitly through the hfi matrix elements as  From equations (4a) and (4b) one can see that the substates Y ñ | 1 and Y ñ | 2 of an arbitrary spin system NV− 13 C in the external magnetic field B||OZ are the basis states ñ º ñ ñ º ñ | | | | III IV 0 , 0 wherein the 13 C n-spin is quantized along the NV symmetry axis and has n-spin projections ( ) m I C =½,−½. In contract, in the general case the substates (4c)-(4f) belonging to the m S =±1 manifolds are the mixtures of the basis states ñ = ñ ñ = ñ | | | | I I I , and ñ = ßñ , with coefficients depending on the hfi matrix elements and on the magnetic field. From this it follows that the 13 C n-spin does not change its projection when the NV e-spin is in the m S =0 state, while it will oscillate between the basis states ñ = ßñ due to the hfi with the NV e-spin having m S =−1 or m S =+1 projections, respectively. One can find within the secular approximation that, for example, in the case of the m S = −1 e-spin state, the 13 C n-spin prepared in the ñ = ßñ | | V state at time τ = 0 can be found in the state Analogously, in the m S =+1 manifold the probability of the 13 C n-spin flipping from the initial state ñ = ñ In an experimental setting, however, one usually measures the average n-spin flip rate by observing the spin system over rather a long time. Under these conditions the time-averaged n-spin flip rates are proportional to G = D B nd 2 2 (see also [24]). At zero magnetic field the 13 C n-spin flips with a rate proportional to the parameter G = 18,24]. Note that by applying a strong external magnetic field B||OZ (γ n B ? A ZZ ) one can substantially reduce the flipping 'rate' G  B in comparison with Γ 0 [17].
It follows from the above consideration that the most important parameter determining the value of the flipping rates Γ 0 , G  B (or the inverse quantities, the Clearly, if the quantity T nd is zero then the stochastic 13 C n-spin flipping dynamics will be absent and the 13 C n-spin in such an NV− 13 C system will keep its state for a long time. As pointed out in the introduction, this property is of great importance for many quantum-technology applications of NV− 13 C spin systems in diamond. In particular, recently, some relatively stable NV− 13 C spin systems were found both in natural [25] and in isotopically-engineered [26][27][28] diamond and were used to implement error correction codes. In all cases, a search was rather time-consuming because it was carried out by a routine systematic study of a very large number of NV− 13 C systems, to find only a few stable ones among them. Evidently, it would be preferable to have information about such systems in advance and we hope we provide it below, using computational chemistry to simulate hfi in the C 510 [NV] -H 252 cluster hosting the NV center in its central part (see figure 1(a)).

Methods, results and discussion
The geometric structure of the C 510 [NV] -H 252 cluster was optimized, and spin density distribution over the cluster was calculated using the DFT/B3LYP/UKS/MINI/3-21 G level of the theory which has been previously [18] used for calculating the hfi characteristics of smaller diamond clusters, for which it showed good agreement with the available experimental data. Calculations have been performed for the singly negatively charged cluster in the triplet ground state (S=1). For geometry optimization, we used the Firefly QC package [29], which is partially based on the GAMESS (US) [30] source code. The full hfi matrices A KL for all possible positions of the 13 C atom in the cluster have been calculated using the ORCA software package [31]. Calculations have been done in the NV-PACS, wherein the Z axis is aligned along the NV center axis while the X and Y axes are chosen arbitrarily. The NV-PACS origin was located at the position of the N atom of the center. Each possible position of the 13 C atom in the cluster was assigned its own number. For example, the positions 1-3 and 4-6 were the nearest-neighbors of the N atom and the vacancy of the NV center, respectively, the positions 7, 8, 469 and 505 were 'axial' positions disposed from the vacancy and N sides, respectively, and so on.
Focusing here on the search for stability positions for the 13 C n-spin in the cluster we calculated the flipping , or respective 'lifetimes' τ 0 =1/Γ 0 , for all possible NV− 13 C systems in the cluster. The results are shown in figure 1 (see also the table in the supplementary information available online at stacks. iop.org/NJP/20/023022/mmedia) which presents the y-logarithmic bar graph of calculated 13 C n-spin 'lifetimes' as dependent on the number of the respective 13 C position in the cluster. One can see that the calculated 'lifetimes' differ considerably for different positions, and that among them there are special stability (or near-stability) positions wherein the 13 C n-spin has a rather large 'lifetime', exceeding that for other positions by two to four orders of magnitude. In figure 1 these relatively stable positions are shown in various colors, to highlight them against the other positions characterized by faster 13 C n-spin flipping rates.
Having calculated, along with the hfi characteristics, the coordinates of all C sites in the relaxed cluster, we determined the locations of the above-mentioned 'colored' stable positions in the cluster, which are shown in figure 2 where the corresponding positions are presented by the same colors as their 'lifetimes' shown in figure 1(b). The four most stable positions-7, 8, 269 and 505, shown in figure 1(b) in red-are just the expected, and previously considered [18][19][20], 'axial' positions. The other near-stable 'colored' positions in figure 1(b) having numbers in the 200s (plus positions 4, 5, and 6, shown in black) are located in the diamond lattice bilayer (see, e.g., [32]) being perpendicular to the NV axis and passing through the vacancy (see figure 2 for better illustration). Among them there are 18 quite stable positions, having 'lifetimes' of between 10 3 and 10 4 , which are shown in the enlarged view in the inset in figure 1(b). They can be classified as belonging to four families [10,18] of near-equivalent positions of 13 C in the cluster, which exhibit near-equal values of experimentallymeasurable hfi characteristics due to the axial symmetry of the NV center. We will refer here to these near-stable families as St1, St2, St3 and St4 below. They contain three, six, three and six members, respectively. More specifically (see the supplement), in our cluster, family St1 consists of the positions C222, C255 and C260; family St2 of positions C223, C225, C256, C263, C269 and C275; family St3 of positions C214, C267 and C277 and family St4 of positions C212, C216, C254, C264, C279 and C286. They are depicted in figure 2 in blue, green, purple and brown, respectively. The same colors are used for the characteristics of the members of these families in the table in the supplement. Spatially, the members of these near-stable families are located symmetrically with respect to the NV axis in the vacancy-containing diamond bilayer at near-equal distances from the axis, as shown in figure 2. Note that the first two families, St1 and St2, located well inside the simulated C 510 [NV] -H 252 cluster, were identified previously [18] in the smaller C 291 [NV] -H 172 cluster, where they were termed the K2 and Y families. Their members are the fourth and fifth neighbors of the vacancy, respectively. One can see from figure 2 that the members of the two additional near-stable families St3 and St4, which are both the seventh Positions exhibiting the longest lifetimes for the 13 C n-spin are depicted in different colors. Red bars show the most long-lived 13 C n-spins located at 'axial' positions, while the other less-stable positions are located basically in the diamond lattice bilayer being perpendicular to the NV axis and passing through the vacancy (see figure 2). The inset shows 18 moderately stable positions having lifetimes exceeding those of the rest by at least two orders of magnitude. They can be divided into four families, termed St1, St2, St3 and St4, containing three, six, three and six members with lifetimes shown by blue, green, purple and brown bars, respectively. neighbors of the vacancy, are located not far from the edge of the cluster, so that their hfi characteristics can be influenced to some extent by the H-terminated cluster surface. Note also that the members of the St1 family which are nearest to the vacancy near-stable positions are situated in the lower sublayer of the vacancycontaining diamond bilayer (β-layer in the terminology of [32]), while all the others are in the higher sublayer (α-layer). It should be emphasized that the 13 C n-spin located in one of the above 18 positions is approximately two orders more stable than the next, less stable positions, shown in figures 1 and 2 in orange and black, and also located in the vacancy-containing diamond bilayer. Additionally, there are nine less stable positions, shown in figures 1 and 2 in yellow, which are located in two higher diamond bilayers disposed above the vacancy. The reason for the exceptionally high stability of 13 C n-spins at these sites is associated with the local symmetry of the spin density distribution in the vicinity of these lattice positions, as will be discussed later. Now we will present in more detail the hfi and spatial characteristics of the above four families of non-axial near-stable NV− 13 C spin systems, St1−St4. These characteristics can be important not only in finding the spin systems experimentally-and in interpreting these findings in a diamond sample-but also in being able to control their creation with increased probability during diamond synthesis by chemical vapor deposition (CVD). The analysis showed that the calculated hfi matrices A KL of different members in the same family were approximately related to each other by the unitary transformation of rotation around the NV axis by angles of 2π/3. In turn, it is easy to verify that hfi parameters A ZZ and T nd are invariant with respect to rotation about the NV axis, and are therefore characteristic for each family along with the observables constructed from these parameters, e.g. the zero-field hfi splitting D = + ( ) / T A ZZ 0 nd 2 2 1 2 of the substates m S =±1 and the hfi-induced 13 C n-spin flipping rates Γ 0 (or 'lifetimes' τ 0 =1/Γ 0 ). These quantities have been calculated for all possible positions of the 13 C n-spin in the cluster (see the supplement) and, in particular, for the St1−St4 families. Simultaneously, we found the spatial characteristics for each member of the near-stable families, viz the Z coordinate of the family member, and its distance from the NV axis and from the N atom of the NV center. Again, the coordinates of the family members were approximately related by the unitary transformation of rotation around the NV axis by an angle of 2π/3. Clearly, due to the finite size of the cluster there was no exact symmetry and, hence, no exact mutual correspondence of the calculated hfi matrices through the abovementioned 2π/3 Z-axial rotations, nor of the spatial characteristics for various members within the near-stable families. Therefore, in  calculation showed that in some cases the rhombicity contribution was essential, and even exceeded the axial component. In particular, for members of the St1 family the rhombicity was only approximately twice as small as the respective axial parts. Additionally, we found all cosines between different axes of the NV-PACS and 13 C-PACS. In table 1 we give only the values (averaged over the members of the St1−St4 families) of the ( ) U d 33 =cosZz, which is the cosine between respective Z and z axes in both coordinate systems.
One can see from table 1 that among axial atoms the highest hfi has the 13 C atom in the C7 position. The other axial C8 position from the vacancy side, and the C469 and C505 positions from the N atom side, exhibited the values A ZZ =86.5, 99.5 and 58.6 kHz. For all of them, the A ZZ values were positive and, due to negligible values of T nd , gave a basic contribution to the respective zero-field hfi splitting, Δ 0 , which can be measured experimentally. Note that our previous calculations [18] for the smaller cluster C 291 [NV] -H 172 predicted a slightly smaller value of Δ 0 =187 kHz for the C7 position. Moreover, in [18]  in the vicinity of these positions because these elements can be presented as (see, e.g., [33]): ò ò g g r g g r is the e-spin density at point  r in the cluster,  R is the position of the 13 C n-spin and S is the total e-spin of the system. Due to the nonlocal character of the anisotropic contribution T, the dipolar integral (5) is over the whole space of the e-spin density distribution, but because of the factor --  | | r R 5 the dominant contribution is made by the e-spin density in the nearest vicinity of the 13 C n-spin. From equation (5) it follows that the dipolar terms vanish if the spin density r +   ( ) r R S n observed from the point of the 13 C nucleus is highly symmetric. As pointed out previously [23], an arbitrary NV− 13 C spin system possesses symmetry C S , which means that in the specific NV-PACS having XZ plane passing through this definite 13 C nucleus and containing the NV axis, the spin density will be symmetric with respect to the change Y to -Y, resulting in the zeroing out of the matrix element T ZY in (5). We verified numerically that in the transformed hfi matrices ¢ A KL for each position of 13 C in our cluster the element T ZY (and, additionally, elements T XY =T YX , T YZ =T ZY ) was zero. For example, for the C222, C255 and C260 positions (belonging to the St1 family), the value of the Z axis rotation angle θ needed to transform from A KL (in the initial NV-PACS) to ¢ A KL (in the aforementioned special NV-PACSs) was found to be 8.91°, 128.049°and 69.042°, respectively.
To understand why the other elements T ZX for the stability positions also vanish, we simulated the spin density distribution over the cluster using the GAUSSIAN'09 program suite [34]. Three examples of calculated Table 1. Calculated hfi and spatial characteristics for the stability positions of the 13 C n-spin in the studied cluster which exhibit the lowest flipping rates Г 0 . Here Δ 0 is the zero-field splitting of the spin sublevels m S =±1 due to hfi. A ZZ and T nd are the calculated hfi matrix elements in the NV-PACS, cosZz are cosines of the angle between the Z axis of the NV-PACS and the z axis of the respective 13 C-PACS wherein the hfi matrix is diagonal, R CN is the distance from the respective 13 C atom to the N atom of the NV center, R CZ is the Z coordinate of the 13 C atoms, and R CXY is the distance from the 13 C positions to the Z axis.   figure 3, with the semi-transparent red/blue lobes corresponding to positive/negative values of the spin density. One can see from figure 3(a) that, as is well known (see, e.g., [9,33,[35][36][37][38][39][40]), the spin density is mostly positive and localized on the three nearest-neighbor 13 C atoms of the vacancy. At smaller absolute isovalues (=0.001 au −3 and 0.0001 au −3 ) there are both positive and negative lobes extending far enough from the Z axis and localized basically near and above the diamond bilayer containing the vacancy of the NV center that is just in the area where the non-axial stability positions are located. Figure 4 shows more clearly the symmetry of the local distribution of spin density at absolute isovalue 0.0001 in the vicinity of the near-stable positions C222, C255 and C260 (belonging to the St1 family) and, additionally, in the vicinity of near-stable positions C223 and C223 of the St2 family. One can see that the lobes of negative spin density near, for example, the C222 position look like axially symmetric bubbles having the axis nearly coinciding with the X axis of the special NV-PACS for NV− 13 C222 where the T ZY element is equal to zero. (More careful analysis showed that the symmetry axis is along the bond from the C222 site to the neighboring C237 site of the diamond bilayer i.e. it consists of the tetrahedral angle 109.5°with the Z axis.) It is because of this axial symmetry of the local spin density distribution around this position that the T ZX element in the hfi matrix for the 13 C nucleus in the C222 position also becomes practically equal to zero, resulting finally in T nd = 0 and, hence, non-flipping 13 C n-spin in this position. Analogous local symmetry of the spin density distribution occurred for the other nearstable non-axial positions of 13 C in the cluster. It is important to note that a vanishing or small perpendicular hfi coupling strength

Position/family
2 which is essential for n-spin robustness, leads in general to a vanishing or slow n-spin manipulation mediated by the e-spin. This circumstance can be overcome when an artificial perpendicular coupling is introduced by a resonant radiofrequency (RF) field.
One can see from figures 3 and 4 that the DFT-calculated spin density distribution over the cluster is rather complicated, especially when considering small isovalue surfaces. Meanwhile, according to figure 3(a), the basic part of the spin density r  ( ) r S is localized near the three carbon atoms being the nearest neighbors of the vacancy of the NV center. It would be instructive to compare the above correct DFT analysis of hfi characteristics in the studied cluster with similar results obtained in the framework of the frequently used point-dipoles approximation where one neglects the NV e-spin distribution. Assuming that the NV spin is completely localized in the middle point between the three nearest neighbors of the NV center vacancy, and using the simulated coordinates of carbon atoms in the relaxed C 510 [NV] -H 252 cluster we calculated (see appendix), in the dipole−dipole approximation, the hfi characteristics for all possible positions of the 13 C n-spin in the cluster. We show that this approximation also predicts the presence of near-stable non-axial positions for 13 C n-spins located in the diamond bilayer containing the vacancy and being perpendicular to the NV axis. However, comparison with the data from the DFT simulation shows the importance of using an accurate, detailed spindensity distribution of e-spin in the diamond clusters hosting the NV center, and demonstrates this with examples of many cases where the simplified dipole−dipole description of hfi in NV− 13 C spin systems gives incorrect results.

Comparison with experiment
As mentioned above, a number of studies [17,[25][26][27][28] have reported the experimental detection of almost-stable NV− 13 C spin systems with different characteristics. In particular, the stable systems exhibiting Δ 0 ≈A ZZ of 201 kHz [17], 89 kHz [26,27] and 50 kHz [28] have been found which were close to the predicted values for the St2 family, axial C8 position and the St3 or St4 families, respectively. To confirm the above theoretical predictions in more detail an additional experiment was undertaken using the near-stable system NV− 13 C of [28], which was used there to enhance quantum metrology by repetitive quantum error correction. The experiment was performed on a single NV center in an engineered diamond with 0.1% 13 C abundancy and a near-stable NV− 13 C spin system exhibiting a lifetime of seconds at a low magnetic field B =340 gauss (B||OZ) was found. To determine the strength of the diagonal hfi component A ZZ , an electron nuclear double resonance experiment (see figure 5(a)) was performed. In this experiment a RF pulse with variable frequency was used to identify the nuclear Larmor frequency in the NV sublevels m s =−1 and m s =0. The results are shown in figure 5(c). As one can see, in the case of the m S =0 substate the signal peaked at ∼362.2 kHz, which is just the nuclear Larmor frequency γ n B of the substate m S =0 in the system NV− 13  To determine the off-diagonal component we used dynamical decoupling spectroscopy [12]. The sequence consisted of periodically-applied microwave pulses, which flip the e-spin of the NV center, and a free evolution between pulses which is adjusted to last half of the nuclear Larmor period. The resonant interaction based on the off-diagonal hfi component is shown in figure 5(d). To determine the strength of T nd we compared the results with numerical simulations and estimated a strength of (1.4±0.1) kHz. The error takes a misalignment of the external magnetic field of the order of 0.1 degrees into account. The simulation was performed in an interaction picture rotating with the e-spin Larmor frequency, and pure dephasing of the electronic part with an adjustable decay rate was included. The finite duration of the microwave pulses was taken into account and the experimental pulse sequence was numerically reproduced using the model Hamiltonian hfi s z 0 n C By inserting the known magnetic field strength as well as the parallel hfi component from above, and testing different values for the decay rate and the perpendicular hfi component, the optimal values describing our experimental results were found. The stability of this n-spin is verified in figure 5(e), which shows a repetitive single shot readout [41] of the n-spin. From such time traces we estimated an n-spin lifetime of 4 s at a magnetic field of 340 gauss. When we compare the experimentally-observed hfi with the results of our DFT simulations, we find good agreement with the stable n-spin groups St3 and St4.

Conclusion
Using DFT, we simulated the H-terminated cluster C 510 [NV] -H 252 hosting the NV center and found, for the first time, the hfi and spatial characteristics of a new class of near-stable NV− 13 C spin systems wherein the 13 C n-spin exhibits negligible hfi-induced flipping rates due to near-symmetric local spin density distribution. Spatially, these positions of stability for the 13 C n-spins are located in the diamond bilayer passing through the vacancy of the NV center and being perpendicular to the NV axis. Analysis of available publications showed that, apparently, some of the predicted non-axial near-stable NV− 13 C systems have already been observed experimentally. A special experiment performed on one of these systems confirmed the prediction made.
We hope that the data in table 1 will help experimentalists find, identify and use these robust NV− 13 C spin systems in emerging quantum technologies. It should also be pointed out that the predicted non-axial nearstable NV− 13 C spin systems can be created during CVD growth of (111) diamond samples by analogy with recent creation [44][45][46] of NV centers perfectly aligned in the (111) direction. Analogous stable n-spins coupled to other paramagnetic centers can also be presented in diamond and in SiC. and A zz ) for the cluster C 510 [NV] -H 252 with a correct account of spatial delocalization of the NV center e-spin over the cluster described by the complicated spin-density distribution r  ( ) r S (figures 3, 4), it would be instructive to compare them with analogous data predicted in the framework of a simplified (but intuitive and therefore popular) description based on the interaction of point magnetic dipoles Following a previous projective measurement [41], the n-spin is in one of the  ñ | eigenstates. The quantization axis is defined by an external, static magnetic field aligned with respect to the NV symmetry axis. After applying an RF pulse (R x ) with a certain duration and a variable frequency, another projective measurement of the n-spin follows and the corresponding n-spin flip probability is calculated. The projective measurement consists of a consecutively-applied block. Each block consists of a CNOT gate (selective, microwave-driven spin-flip of electron when carbon spin is in ñ | -state), which correlates the n-spin with the e-spin, and an optical readout of the e-spin, which also polarizes the spin. The n-spin state is then verified after the accumulation (N=10 000) of NV fluorescence. (b) Measurement sequence for determining the non-diagonal hfi component T nd . Firstly, a microwave pulse (R π/2 ) is applied to create e-spin coherence. Secondly, a periodical sequence of microwave pulses is applied. Each pulse (R π ) flips the e-spin. The time period of the pulse sequence is 2τ. In addition, the rotation axis of the pulses is varied to fulfill robustness against pulse error (XY8 sequence see [42,43]). : where, as previously, γ e and γ n are the gyromagnetic ratios for the NV center and 13 C, respectively,  S and  I are the e-and n-spin operators,  R is the NV→ 13 C vector (coordinates x y z) and  T is the dipole−dipole interaction tensor without averaging over e-spin distribution. As was pointed out in the main text (see also [9,[35][36][37][38][39][40]), the e-spin of the NV center is mostly localized on the three nearest-neighbor carbon atoms of the vacancy. Therefore, within the model of interacting point magnetic dipoles, we assume that the NV spin is completely localized in the middle point of the three carbon atoms (atoms C4, C5, C6 in the table in the supplement) and have the coordinates (0, 0, 2.29) Å. Using the simulated coordinates of carbon atoms in the relaxed C 510 [NV] -H 252 cluster and the values γ e =1.8574 * 10 −20 erg gauss −1 and γ n =7.0944 * 10 −24 erg gauss −1 we calculated the dipole−dipole matrices  T for all NV− 13 C systems in the cluster using the expression (see also the appendix in [48] where a nn =1.54 Å is the nearest-neighbor distance for diamond, and matrix elements are in MHz. Matrices (A2) are the hfi matrices in the dipole−dipole approximation (with elements in MHz). Figure A1 A similar bar graph showing diagonal elements T zz of the matrix (A2) for different 13 C positions in the cluster is presented in figure 2(a). Again for comparison, figures 2(b) and (c) show two (dipolar-like T zz and isotropic A iso ) DFT-simulated contributions to the diagonal element A zz of the hfi matrices.
One can see from these two figures that, for many remote positions of 13 C in the cluster, the approximation of the point dipoles works well and its predictions are close to the results of the DFT simulation. At the same time, for other positions its predictions differ substantially from those of the DFT simulation. Evidently, the reason is in the complicated spin-density r  ( ) r S distribution over the cluster, which is not taken into account within the framework of a simple dipole−dipole approximation. Moreover, the dipole−dipole approximation does not take into account the isotropic contribution to hfi which, according to figure A2(c), can be essential even for rather distant positions of the 13 C n-spin.  In figure A3 we compare the main result of our main article, figure 1, with the predictions of dipole−dipole calculations.
One can see from figure A3 that dipole−dipole approximation also predicts the presence of near-stable nonaxial positions for 13 C n-spins located in the diamond bilayer containing the vacancy and being perpendicular to the NV axis. In table A1 below we give the values of non-diagonal (T nd ) and diagonal (T zz and A iso ) elements of hfi matrices for 18 near-stable positions of 13 C n-spin (members of the St1−St4 families, see main text), calculated in the dipole−dipole approximation and using DFT. One can see that for the St1 family the dipolar approximation does not work well due to the close proximity of the NV center, whereas for the more distant St2 −St4 families it gives more-or-less similar results for dipolar contributions to the hfi matrices. Again, in all 18 cases the isotropic contributions were essential. It should also be noted that, according to figure A3(a), the dipole −dipole approximation prediction contradicts that of DFT in the number of near-stable non-axial positions, viz, along with the 18 near-stable positions belonging to the St1−St4 families and described both in the main text and in table A1, the dipole−dipole approximation predicts many additional near-stable positions which are absent in the DFT description.
In conclusion, the above analysis shows the importance of obtaining an accurate description of the complicated spin-density distribution of e-spin in the diamond clusters hosting the NV center, and demonstrates this with examples of many cases where the simplified dipole−dipole description of hfi in NV− 13 C spin systems gives incorrect results. Table A1. Comparison of non-diagonal and diagonal elements of the hfi matrices for stable positions of 13 C n-spin in the cluster, belonging to the St1−St4 families (see main text), calculated using DFT simulation and a dipole−dipole approximation.