Spin decontamination for magnetic dipolar coupling calculations: Application to high-spin molecules and solid-state spin qubits

An accurate description of the two-electron density, crucial for magnetic coupling in spin systems, provides in general a major challenge for density functional theory calculations. It affects, e.g., the calculated zero-field splitting (ZFS) energies of spin qubits in semiconductors that frequently deviate significantly from experiment. In the present work, (i) we propose an efficient and robust strategy to correct for spin contamination in both extended periodic and finite-size systems, (ii) verify its accuracy using model high-spin molecules, and finally, (iii) apply the methodology to calculate accurate ZFS of spin qubits (NV− centers, divacancies) in diamond and silicon carbide. The approach is shown to reduce the dependence on the used exchange-correlation functional to a minimum.

The phenomenological spin Hamiltonian [1] is conventionally adopted to characterize the ZFS in terms of the parameters D = D zz − 1 2 (D xx + D yy ) and E = 1 2 (D yy − D xx ), where D ii are the eigenvalues of the symmetric and traceless 3 × 3 tensor D, and axially symmetric D tensors lead to vanishing E .
Within the DFT framework, the spin-dipolar coupling driven D tensor (referred to as spin-spin ZFS) is derived from the spatial distribution of Kohn-Sham orbitals obtained with standard spin-polarized self-consistent field calculations [32][33][34][35][36]. For the S 1 high-spin state of a system, the spin-spin D tensor is given by where the elements of the traceless 3 × 3 matrix d are constituted by a magnetic dipolar interaction between all the pairs of occupied Kohn-Sham states m and n (with a, b = x, y, z), 1. (a) Schematic illustration of the electron-electron magnetic dipole interaction in a spin-triplet system. The overall coupling d originates from the interaction between the two single occupied Kohn-Sham orbitals (denoted as I), as well the contributions arising from the double occupied orbitals (II). (b) A set of biradicals of different length (general formula C n H 2n ; adopted from Ref. [4] [37], PBE functional [38], and triple-zeta basis set with one polarization function, def2-TZVP [39]) with the direct (red crosses) and UNO-corrected (black crosses) approaches are presented.
In Eq. (3), α is the fine-structure constant, the charge densities n mn (r) = ψ m (r)ψ n (r) reflect the spatial distribution of the orbitals, and χ m,n originates from the matrix elements of spin operators: χ m,n = 1 when the orbital belongs to the same spin channel, and χ m,n = −1 otherwise. Throughout this Rapid Communication, we will also make use of the magnetic dipolar coupling parameter d = d zz − 1 2 (d xx + d yy ) analogous to the D value. Note that the expression in the square brackets in Eq. (3) is an approximation of the two-particle spin density which is otherwise not directly available from DFT.
Due to spin polarization, the Kohn-Sham states of a paramagnetic system have different energies and different spatial distributions in the spin-up and spin-down channels. Therefore, the spin-spin ZFS is not entirely determined by the coupling between the half-filled states [denoted as d I,I in the energy level diagram in Fig. 1(a), where a spin-triplet S = 1 case is shown as a prototype example]. Instead, it also contains a nonvanishing collective contribution from the two-particle spin densities of single-/double-(d I,II ) and double-/doubleoccupied (d II,II ) states. Note that in the absence of spin polarization, the contributions from the spin-up and spin-down electrons of the same double-occupied orbital would cancel out due to the spin-dependent factor χ m,n .
The spin contamination manifests itself in unphysical spindipolar interaction which induces spurious "on-site" terms in the two-particle spin density [4] entering d I,II . This can be nicely demonstrated by adopting the family of S = 1 biradicals used in Ref. [4] as a reference system [see Fig. 1 Here, the unpaired electrons are located at the opposite termini, so that the spin-spin ZFS must decrease according to the point-dipole approximation as the molecule gets longer and longer. This is, however, not the case as the unphysical on-site contributions [d I,II T in Fig. 1(c), where the subscript The considered biradicals represent a lucky case where the physical and spin contamination driven parts of the magnetic dipolar coupling tensor are well separated. Indeed, for the longer molecules, the coupling is completely described by d I,I T , while d I,II T is being entirely caused by spin contamination. In general, there is, however, no evidence that the spin contamination can be attributed only to selected terms in Eq. (3) or that it can be treated as an additive error. Therefore, it must be excluded in a more elaborate way. Here, we establish a systematic correction scheme based on the simple idea that the spin-contamination error in d (m S =S) of the m S = S high-spin state of a general S 1 system is one-to-one reflected by its m S = S − 1 state. In Fig. 1(c), this is illustrated for the prototypical S = 1 case. Its m S = 0 low-spin configuration is a single-Slater-determinant state and thus should be referred to as the broken-symmetry (BS) state. As a consequence, the average d (m S =S−1) among all the possible m S = S − 1 configurations restores the physical part while canceling out the spin-contamination error completely [see also the discussion in the Supplemental Material (SM) [40]], Each of these m S = S − 1 configurations can be obtained by changing the occupation of one of the half-filled Kohn-Sham orbitals from spin up to spin down and subsequently performing the self-consistent field calculation. Note that the S-dependent denominator in Eq. (4) is exactly the prefactor which relates the D values and the energetic distance between the m S = S and m S = S − 1 spin sublevels [40].
For an S = 1 system, this follows directly from Fig. 1(c). While showing no net spin, its m S = 0 BS state exhibits antiferromagnetic coupling between the spin-up and spindown electrons occupying two half-filled orbitals, which leads to nonvanishing d I,I BS [blue triangles in Fig. 1(c)]. Thus, in a "perfect" spin-triplet system, i.e., a system without spin contamination, d BS = −d T , is expected. Consequently, when the spin-spin ZFS of a triplet system is defined asD SS = (d T − d BS )/2, the physical d I,I T part is conserved, while the spin-contamination errors of the triplet and BS states cancel.
As an illustration, we apply the proposed strategy to the considered set of biradicals. Here and in the following, the spin-spin ZFS calculations [35] are based on the projector augmented-wave (PAW) formalism [41] in combination with norm-conserving pseudopotentials [42], as implemented in a modified version of the gauge-including PAW (GIPAW) module of the QUANTUM ESPRESSO software [43,44]. We used the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [38], and a plane-wave (PW) basis set with 700 eV kinetic energy cutoff. We find that the corrected ZFS exhibits perfect agreement with the UNObased all-electron DFT results [ Fig. 1(d)]. Most importantly, in contrast to the UNO approach, the established scheme is applicable for high-spin systems with periodic boundary conditions. Subsequently, we apply the correction to the NV − and VV 0 defects in 4H-SiC, and the NV − center in diamond. These spin-triplet defect centers are arguably among the most prominent solid-state spin qubits. Their ZFS is believed to be almost entirely caused by magnetic dipolar coupling and has been successfully treated with the plane-wave pseudopotential DFT [8,10,31]. Some recent publications [35,45,46], however, argued that the reported agreement between the calculated and the measured ZFS profits from error cancellation due to incomplete treatment of the pseudized electron density. Indeed, as shown in Ref. [35], the resulting ZFS depends significantly on the technical parameters adopted for the pseudopotential generation, unless the "true" character of the wave function within the atomic core region is reconstructed. Such a reconstruction can be achieved by the PAW method [41]. As pointed out by Ivády et al. [46], the fully PAW reconstructed ZFS of the spin-triplet centers in the SiC polytypes is, however, almost 30% larger than the values reported experimentally.
In diamond, there is only one magnetically distinct NV − center. It can be aligned with one of the 111 diagonals and always exhibits C 3v symmetry. In contrast, the hexagonal 4H polytype of SiC features a periodic sequence of quasicubic (k) and hexagonal (h) Si-C double layers along the crystal c axis, offering inequivalent crystallographic positions and configurations for the defect pairs. Those NV − and VV 0 centers that are aligned with the hexagonal c axis (the so-called axial pairs) have C 3v symmetry [see Fig. 2(a)]. For all other orientations, the symmetry is lowered towards C 1h . Irrespective of the host material and the crystallographic position, the NV − and VV 0 defects share almost the same m S = 1 electronic structure with two unpaired electrons either populating a double degenerate e orbital (axial pairs) or slightly split a and a orbitals (basal pairs). As shown in Fig. 2(c) and Table I, even for the axially symmetric NV − and VV 0 centers, the ZFS depends on the lattice site sensitively. Notably, our results for SiC reproduce the previously reported 30% discrepancy between the measured D values and D SS , if directly calculated via Eq. (2). It should be also mentioned that the coupling of the two half-filled Kohn-Sham states (i.e., the d I,I T term) explains in this case less than 70% of the measured value, rendering a spin-restricted approach (where the same wave functions are used for spin up and spin down) to be insufficient. The d I,II T and d II,II T terms are thus significant for NV − and must be retained in a systematic manner. In the following, we show that this can be accomplished by using the proposed methodology.
In contrast to the biradicals considered above (Fig. 1), the unpaired electrons of the axial NV − or VV 0 center are equally distributed among the three carbon dangling bonds surrounding the vacancy [cf. Fig. 2(a)]. Thus, in a m S = 0 configuration, the spin-up electron can localize at one of these dangling bonds, while the spin-down electron is being TABLE I. The spin-spin ZFS for the NV − center in diamond, and NV − , VV 0 , V − Si centers in 4H -SiC calculated without (second column, D SS ) and with (third column,D SS ) a spin-contamination correction in comparison with experiment. Experimental values taken from Ref. [31] (diamond), Refs. [8,10] (NV − in SiC), Ref. [13] (VV 0 ), and from Ref. [21] (V − Si ). The calculations performed using 512-atom (diamond) and 432-atom (4H -SiC) supercells with shifted 2 × 2 × 2 k-point grids. shared by the other two adjacent carbon atoms. As illustrated in Fig. 2(b), this provides three symmetry-reduced BS states of the axial defect. Thus, in order to apply Eq. (4) to the ZFS of NV − and VV 0 , we define the d (m S =0) tensor as the C 3v symmetrized average of the three resulting BS tensors. Figure 2(c) demonstrates that the spin decontamination applied in this way removes the overestimation and brings the calculated ZFS (D SS ) in almost perfect agreement with the measured D values. The improvement is especially decisive for the spin qubits in SiC. Note that a part of the remaining discrepancies (<100 MHz) can be attributed to second-order contributions to the ZFS due to spin-orbit coupling [47].
The same procedure can be applied in the case of the basal pairs in SiC. From Table I it can be seen that these spin qubits suffer from the spin contamination to the same extent as their axially symmetric counterparts. Averaging of their three no-longer equivalent BS configurations provides improved estimates for D, and also reasonable values for the rhombicity parameter E . Finally, we show that the proposed approach can be successfully used for higher-spin qubits, i.e., the V − Si centers in 4H-SiC that possess a S = 3/2 ground state. Notably, the D value of V − Si at the h site is found to manifest an enormous spin-contamination error, which is eliminated by our correction scheme (cf. Table I Fig. 1) calculated with different exchange-correlation functionals: Perdew-Zunger (PZ) [48], PBE, Perdew-Wang (PW91) [49], Becke three-parameter Lee-Yang-Parr (B3LYP) [50], and Heyd-Scuseria-Ernzerhof (HSE) [51]. The horizontal dashed line indicating theD SS value obtained with the PBE functional is provided as a guide for the eye.
accuracy provided by our approach, the reported results can be considered as a solid contribution that supports the assignment of the spectroscopic fingerprints observed in 4H-SiC as NV − , VV 0 , and V − Si centers. To summarize, in the present Rapid Communication we propose an efficient scheme to eliminate spin contamination in magnetic dipolar coupling. It allows us to identify the spin contamination of the two-electron density as a highly relevant source of a discrepancy between measured and calculated ZFS of high-spin defects in semiconductors. This discrepancy can be lifted by calculating the difference of m S = S and m S = S − 1 magnetic dipolar coupling tensors. This scheme is robust with respect to the exchange-correlation functional (see Fig. 3). In particular, it compensates the increased extent of spin contamination associated with the introduction of a fraction of an exact Hartree-Fock exchange in hybrid functionals [5,52,53]. This opens up the possibility to use either hybrid functionals or (semi)local [local density approximation (LDA) and generalized gradient approximation (GGA)] functionals while reaching the same level of accuracy for the spin-spin ZFS. The proposed correction will be also beneficial for an accurate description of magnetic dipolar interactions between distant spin centers, e.g., for modeling a qubit coupled to a spin bath [54][55][56] or a solid-state spin probe coupled to surface electron spins [57].
Numerical calculations were performed using grants of computer time from the Paderborn Center for Parallel Computing (PC 2 ) and the HLRS Stuttgart. The Deutsche Forschungsgemeinschaft (DFG) is acknowledged for financial support via the priority program SPP 1601 and the Transregional Collaborative Research Center TRR 142 (Project No. 231447078).