In the past years, there has been rapid progress in the development of quantum magnetic sensing technologies based on nitrogen-vacancy (NV) color centers in diamonds (e.g., see (Schwartz, 2019; Barry et al., 2020; Pezzagna and Meijer, 2021) for reviews). The use of single NV centers makes it possible to implement a magnetometer that provides nanometer-scale spatial resolution and extraordinary sensitivity, allowing the detection and imaging of multiple individual ¹³C nuclear spins (see, e.g., Zhao et al., 2012; Kolkowitz et al., 2012; Taminiau et al., 2012; Dréau et al., 2012; Mȕller et al., 2014; Zopes et al., 2018; Sasaki et al., 2018; Abobeih et al., 2019; Cujia et al., 2022; Vorobyov et al., 2022; van de Stolpe et al., 2023) inside the diamond. Due to their large coherence times, resulting from high isolation from the environment, these ¹³C nuclear spins are widely used as a quantum memory in emerging quantum technologies (see, e.g., Awschalom et al., 2018). In addition to single ¹³C nuclear spins, it has also been proposed (Reiserer et al., 2016; Chen et al., 2017) to use for the purpose pairs of coupled ¹³C nuclear spins or ¹³C–¹³C dimers (Zhao et al., 2012), since, in singlet state, such spin systems demonstrate exceptionally long coherence times (Stevanato et al., 2015; Levitt, 2019). Experimentally, individual dimers were observed and analyzed earlier in diamond with NV centers (Zhao et al., 2011; Shi et al., 2014; Ma and Liu, 2016). Recently, they have been studied in detail in (Abobeih et al., 2018; Bradley et l., 2019; Yang et al., 2020; Bartling et al., 2022; van de Stolpe et al., 2023) where a single NV center was used not only to detect and to characterize multiple ¹³C nuclear spins and few individual ¹³C–¹³C dimers in the specific spin environment of the studied NV center with sub-Hertz spectral resolution but also to initialize, control and readout the states of these dimers (Abobeih et al., 2018; Yang et al., 2020). In particular, in (Bradley et al., 2019), it was shown that the inhomogeneous dephasing time for the studied ¹³C–¹³C dimers was about 1 min at room temperature, the longest reported for individually controlled qubits. It should be noted also, that the parameters of the ¹⁴NV center itself can now be measured at Hz-level precision (Xie et al., 2021).

In the works (Zhao et al., 2011; Shi et al., 2014; Ma and Liu, 2016; Abobeih et al., 2018; Bradley et al., 2019; Yang et al., 2020; Bartling et al., 2022; van de Stolpe et al., 2023), the experimental results have been interpreted on the basis of theoretical analysis of a quantum system NV-¹³C–¹³C with the account of the Zeeman interaction of its spins with an external magnetic field, the hyperfine interaction of the electron spin of the NV center with the ¹³C nuclear spins forming the dimer, and of the direct dipole-dipole interaction of ¹³C nuclear spins. The magnetic field lifts the degeneracy of the NV center states with projections m_S = ±1 and initiates the precession of the ¹³C nuclear spins of the dimer, the frequencies of which depend on the projection m_S of the electronic spin S of the NV center. When the NV center is in the m_S = 0 state, both spins of the dimer do not hyperfine-coupled to the NV center and have the same Larmor frequencies (the case of equivalent nuclear spins in the Nuclear Magnetic Resonance (NMR) terminology (Wasylishen, 2007). Taking into account the interaction of nuclear spins with each other, the dimer can be considered (Wasylishen, 2007) as four-state system involving three triplet states ↑ ↑ , ↑ ↓ + ↓ ↑ / 2 , ↓ ↓ and one singlet state ↑ ↓ − ↓ ↑ / 2 , where the arrows mean the projections of ¹³C nuclear spins. In this case only two transitions ↓ ↓ ↔ ↑ ↓ + ↓ ↑ / 2 , ↑ ↓ + ↓ ↑ / 2 ↔ ↑ ↑ between triplet states are possible, giving two lines in the NMR spectrum, known as the Pake doublet (Pake, 1948). The two lines separation is determined mainly by the direct dipolar coupling of nuclear spins, the magnitude of which depends on the distance between them and on the orientation of the dimer relative to the external magnetic field. Typical values of the direct dipolar coupling for neighbor spin pairs range from a few hundred Hz to several kHz (Wasylishen, 2007). In the case of ¹³C nuclear spins in diamond this coupling is ∼2 kHz (Shi et al., 2014; Bartling et al., 2022) when the dimer is oriented parallel to the magnetic field, and approximately 3 times less (∼690 Hz) if it is oriented at a tetrahedral angle of 109.5° relative to the external magnetic field (Shi et al., 2014). For the dimer consisting of non-nearest neighbors, the magnitude of the interaction drops to ∼200 Hz (Bartling et al., 2022).

Different situation takes place, when the NV center is transferred to the m_S = ±1 states, for example, by microwaves having resonant frequency for respective transition. In these cases both nuclear spins of the dimer experience a hyperfine interaction with the electron spin of the center, which can either differ significantly (Shi et al., 2014) or be quite similar (Bartling et al., 2022) in characteristics depending on the position of the dimer with respect to the NV center. As a result, the ¹³C nuclear spins belonging to the dimer cease to be equivalent and four-level nuclear-spin subsystems of a dimer with the NV center in m_S = ±1 states are not purely triplet and singlet. In such cases, NMR transitions between all nuclear-spin states of a NV-¹³C–¹³C spin system become possible, giving the corresponding four lines in the NMR spectra of the NV center in each of states m_S = ±1. In addition, in a NV-¹³C–¹³C system, microwaves can initiate transitions between different nuclear-spin sublevels of the states of the NV center with projections, for example, m_S = 0 and m_S = −1. The frequencies and amplitudes of all these transitions are determined by the strength and direction of the magnetic field acting on the system, the hyperfine interactions of the NV center with the nuclear spins of the dimer and the interaction of these nuclear spins with each other. The accuracy of characterization of dimers achieved in the above works, in particular, the measurement of their resonant frequencies, requires a detailed analysis and consideration of all interactions actually presenting in such spin systems.

The most obvious additional spin-spin interaction, which was not taken into account earlier, is indirect nuclear spin–spin interaction (J-coupling) that arise due to second-order hyperfine interactions with electrons from chemical bonds connecting the studied nuclei (see, e.g., (Ramsey, 1953; Vandersypen and Chuang, 2004; Wasylishen, 2009)). Usually this interaction is expected to be weaker in comparison with the above mentioned ones, but it still can have a significant impact on the evolution of spin systems. As far as the J-coupling is mediated by electrons, it is sensitive to subtle changes in geometry, conformation and electronic structure of molecules and solids, thus being a valuable source of important structural and dynamical information usually revealed by NMR spectroscopy. Typically, the strength of J-couplings is about tens or hundreds Hz for proton and carbon nuclei separated by a single covalent bond, and quickly decreases as the number of inter-mediate bonds grows.

Generally, a second-rank tensor ⁿ J _KL(K, L = X, Y, Z) is required to fully describe J-coupling between two nuclei (Harris et al., 2009). However, until recently, most conventional high-resolution NMR experiments were focused on measuring only isotropic scalar constant ⁿ J _iso = Sp(ⁿ J _KL)/3 because the anisotropic parts of the J-tensor were averaged out to zero by fast molecular motion in solution-state NMR or fast magic-angle spinning (MAS) in solid-state experiments (Frydman, 2001; Vaara et al., 2002; Harris et al., 2009; Reif et al., 2021). In the case of crystalline solids, the constituent atoms are located in a certain order determined by the crystal structure so that many important NMR interactions are orientation-dependent, and the information about anisotropic NMR interaction tensors becomes essential (Vaara et al., 2002; Harris et al., 2009). In particular, both the symmetric ¹ J _iso =(¹ J _ZZ+¹ J _XX+¹ J _YY)/3 (=70 Hz) and the asymmetric Δ¹ J = ¹ J _ZZ -(¹ J _XX +¹ J _YY )/2 (=90 Hz) parts of the J-coupling tensor ¹ J _KL for the nearest-neighbor (N–N) ²⁹Si nuclear spins have been determined in single-crystal silicon (Christensen and Price, 2017). This was achieved by measuring the NMR lineshapes, which are sensitive to the value of Δ¹ J, at four different crystallographic orientations relative to the applied magnetic field.

In diamond with ¹³C nuclear spins, a similar experiment was performed many years ago (Lemann et al., 1994), however, at that time, any effects associated with the presence of indirect interactions of ¹³C spins could not be studied due to the low detection sensitivity. Currently, owing to the development of new highly sensitive NV-based methods of NMR spectroscopy (Boss et al., 2017; Schmitt et al., 2017; Bucher et al., 2020) used in the above-mentioned papers (Abobeih et al., 2018; Bradley et l., 2019; Yang et al., 2020; Bartling et al., 2022; van de Stolpe et al., 2023), direct observation and study of J-coupling effects have become possible. In this respect, predicting their characteristics for studied spin systems is essential. Up to the recent time, there were no theoretical works on the quantum-chemical simulation of J-coupling of ¹³C nuclear spins in diamond. The first such attempt was made in our article (Nizovtsev et al., 2022), where we used a rather simple density functional theory (DFT) level to simulate full tensors J _KL describing the J-couplings of nuclear spins ¹³C in small H-terminated diamond clusters, as well as in a cluster hosting the NV-color center. Here, we extended findings using an additionally more advanced DFT level and presented our results in conjunction with those obtained in (Nizovtsev et al., 2022).

Basically, the theoretical foundations of J-coupling are well established (Ramsey, 1953; Wray, 1979; Helgaker et al., 1999; Krivdin, 2004; Helgaker et al., 2008; Krivdin, 2021), and there has been considerable progress in calculating the J-coupling characteristics for many simple molecules (e.g., see (Wray, 1979; Kamienska-Trela, 1995; Helgaker et al., 1999; Krivdin, 2004; Krivdin and Contreras, 2007; Antušek et al., 2008; Krivdin, 2018), including ¹³C–¹³C pairs (Wray, 1979; Kamienska-Trela, 1995; Jaszunsi et al., 2003; Krivdin, 2004; Peralta et al., 2004; Krivdin, 2018). However, previously, most software packages were aimed specifically at the calculation of the scalar J-coupling constants. Only recently has it become possible to get information about full J-coupling tensors. Here, we have used for the purpose version 5.0.3 of the ORCA package (Neese, 2012; Neese, 2022). The diamond crystal was modeled by H-terminated carbon clusters. We optimized the cluster geometry with the package at the B3LYP/UKS/def2/J/RIJCOSX level of theory and then simulated the n-bond J-coupling tensors ⁿ J _KL(Ci, Cj) for all possible ¹³Ci–¹³Cj pairs in the clusters using two levels of theory. The first one is the B3LYP/UKS/TZVPP level of theory which was used in (Nizovtsev et al., 2022) according to (Makhyoun et al., 2019) where it was shown that the functional B3LYP in combination with TZVPP doubly polarized triple-zeta basis set is quite successful in comparative calculations of J-coupling parameters performed using different levels of theory. Additionally, here we have also used the PBE0/UKS/pcJ-2 level of theory with standard DFT hybrid level PBE0 employing a special pcJ-2 segmented contracted basis set (Jensen, 2008), which was specially developed for spin–spin coupling constants calculation. The successful application of the ORCA software with the second level of theory was demonstrated in (Grimme et al., 2017).

In order to test the ORCA opportunities, we first calculated the J-tensors for all possible pairs ¹³C–¹³C in the diamond-like adamantane molecule C₁₀H₁₆ (see Figure 1A), for which the isotropic J-coupling constants ¹ J _iso for N-N nuclear spins ¹³C was experimentally measured to be 31.4 ± 0.5 Hz (Gay et al., 1991). Having obtained for them the values of ∼29.9 Hz and ∼30.1 Hz (illustrated in Figure 2A) when using basis sets 1 and 2, respectively, which both were quite close to the experimental value, we performed similar calculations for the H-terminated carbon cluster C₃₅H₃₆ (Figure 1B), as well as for the cluster C₃₃[NV⁻]H₃₆ hosting the NV color center (Figure 1C). It should be noted that the choice of these small clusters was due to the fact that, as is known (Krivdin, 2004; Helgaker et al., 2008; Krivdin, 2021), the calculations of the J-coupling characteristics are very computationally demanding for even modest-sized molecules.

The package returns matrices describing the diamagnetic, paramagnetic, Fermi-contact, spin-dipolar, and spin-dipolar/Fermi contact cross-term contributions to the total ⁿ J _KL tensor in the coordinate systems indicated in Figure 1. Using them and taking into account the known coordinates of carbon atoms belonging to some definite ¹³Ci–¹³Cj pair in the cluster, one can find respective J-coupling matrices in the other coordinate system. In particular, for neighboring nuclear spins ¹³C, separated by a single bond in diamond (∼1.54 Å), the total ¹ J _KL matrix becomes diagonal with J _XX ≈ J _YY in the coordinate system in which the Z-axis is directed along this bond (Christensen and Price, 2017). In this case, it is conventional to describe an axial J-coupling tensor in terms of two parameters: the scalar constant ¹ J _iso and the asymmetric part Δ¹ J. Since the magnitude of the J-coupling decreases rapidly with bond order, we will mainly consider here such N–N nuclear spins.

Using quantum chemistry software ORCA at two theory levels we simulated full tensors describing the indirect interaction (J-coupling) of ¹³C nuclear spins in H-terminated diamond cluster C₃₅H₃₆ and in the cluster C₃₃[NV⁻]H₃₆ hosting the NV center. We found that the PBE0/UKS/pcJ-2 level of theory provides, in the case of adamantine molecule, better agreement with the only available experimental data on the isotropic constant ¹ J _iso in comparison with previously used (Nizovtsev et al., 2022) more simple and less time-consuming B3LYP/UKS/TZVPP level of theory. At the same time, it is shown that the calculations of full J-coupling tensors J _KL performed using these two different levels of theory give fairly close results, which, in particular, is important in itself for developing methods for calculating the total J-coupling tensors.

Here, we have focused mainly on the calculation of the full tensors J _KL for nuclear spins ¹³C in diamond separated by one bond for which the J-coupling is strongest. For such ¹³C–¹³C dimers, the calculated J-coupling tensors are almost diagonal in the coordinate system in which the Z-axis is directed along the bond connecting the nuclei. We found that in addition to the usually considered isotropic scalar ⁿ J-coupling constant, the anisotropic contributions to the ⁿ J-coupling tensor are essential. It is also shown that the NV center affects the characteristics of the J-coupling of ¹³C nuclear spins, especially if they are located near the vacancy of the NV center.

The obtained data on the J-coupling of ¹³C nuclear spins, being supplemented by the previously simulated data (Nizovtsev et al., 2018) on the hyperfine interaction of the NV center with the ¹³C nuclear spins forming some specific dimer, as well as on the spatial locations of these nuclear ¹³C spins relative to the NV center, forms a complete set of data allowing to simulate numerically energy levels and eigenstates of studied spin systems by analyzing their spin-Hamiltonians and to predict frequencies and strengths of transitions between states of the system. Such analysis of multi-spin systems NV–¹³C–¹³C with the account of both direct dipole–dipole interactions of ¹³C nuclear spins and of their J-coupling has been done recently in (Nizovtsev et al., 2023). Simulating numerically the specific spin system we demonstrated that for accurate determination of the Pake doublet splitting in a state with the spin projection of the NV center m_S = 0, along with the main dipole-dipole interaction of ¹³C nuclear spins in the dimer, it is necessary to take into account their anisotropic J-coupling. The account of the J-coupling can also be important in cases where the gradient of hyperfine interactions for the nuclear spins of the studied dimer is small. In addition, we also performed simulations of their EPR and NMR spectra and discuss effective ways to convert the studied dimer to the desired long-lived singlet state using microwave and radiofrequency pulses having predicted characteristics.

Theoretical simulation of J-coupling tensors is also important in conjunction with studies aimed at the creation of nanoscale NV-based quantum sensors for the detection of molecules/radicals adsorbed on the surface of nanostructured diamond (see, e.g., Glenn et al., 2018) and the determination of their chemical structure. The data obtained can also be useful for studies of NMR in the zero-to ultralow-field (ZULF) regime (Theis et al., 2013; Blanchard and Budker, 2016; Jiang et al., 2018; DeVience et al., 2021), where the internal spin interactions are dominated in their natural environment.