# Dipolar spin relaxation of divacancy qubits in silicon carbide

Oscar Bulancea Lindvall1 Nguyen T. Son1 Igor A. Abrikosov1 Viktor Ivády1,2,3,∗

###### Abstract

Divacancy spins in silicon carbide implement qubits with outstanding
characteristics and capabilities in an industrial semiconductor host. On the
other hand, there are still numerous open questions about the physics of
divacancy point defects, for instance, spin relaxation has not been thoroughly
studied yet. Here, we carry out a theoretical study on environmental spin
induced spin relaxation processes of divacancy qubits in 4H-SiC. We reveal all
the relevant magnetic field values where the longitudinal spin relaxation time
T${}_{\\_}1$ drops resonantly due to the coupling to either nuclear spins or
electron spins. We quantitatively analyze the dependence of the T${}_{\\_}1$
time on the concentration of point defect spins and the applied magnetic field
in the most relevant cases and provide an analytical expression. We
demonstrate that dipolar spin relaxation plays a significant role both in as-
grown and ion implanted samples and it often limits the coherence time in
4H-SiC.

Department of Physics, Chemistry and Biology, Linköping University, SE-581 83
Linköping, Sweden

Wigner Research Centre for Physics, Hungarian Academy of Sciences, PO Box 49,
H-1525, Budapest, Hungary

Max-Planck-Institut für Physik komplexer Systeme, Nöthnitzer Street 38,
D-01187 Dresden, Germany

Correspondence should be sent to<EMAIL_ADDRESS>
August 27, 2024

## Introduction

Through the example of the nitrogen-vacancy center in diamond[1, 2, 3, 4] (NV
center), point defects in wide-band gap semiconductors have demonstrated their
potential for quantum-enhanced technologies. In particular, the NV center
based devices are about to revolutionize sensing at the nanometer scale[5, 6,
7, 8, 9, 10, 11, 12], while NV centers coupled to adjacent nuclear spins can
serve as nodes for quantum internet[13, 14, 15] and quantum computations[16,
17, 18]. While the physics of the NV center is understood to a very large
degree, other point defect qubits with similar or even superior capabilities,
such as the neutral silicon and carbon vacancy pair (divacancy)[19, 20, 21]
and the negatively charged silicon vacancy (silicon vacancy)[22, 23] in
silicon carbide (SiC) are still subjects of active fundamental studies. For
example, recent reports on the neutral divacancy have demonstrated 64 ms
coherence time[24], implementation of spin-to-photon interface[25], nuclear
spin operations for quantum memory applications[26, 27] and dynamic nuclear
polarization[26, 28, 29], and a room temperature spin contrast as high as
30%[30]. The fact that SiC is a technologically mature semiconductor with
controllable p- and n-type doping and established nanofabrication techniques
further enhances its importance and paves the way toward affordable integrated
quantum devices[31, 32, 33].

Besides coherent properties of the qubits, longitudinal spin relaxation, with
the corresponding decay time T${}_{\\_}1$, is of great importance as it sets
the fundamental limit for several applications, e.g. for dynamic decoupling
techniques and sensing.[34, 35, 36, 37] Relaxation of the NV center’s electron
spin at different magnetic fields, temperatures, and environmental spin defect
concentration has been studied.[38, 39] Furthermore, optical signatures of
strong environmental couplings were revealed, enabling novel microwave free
spectroscopy applications.[40, 41, 42, 43, 44] In contrast, much less is known
about the longitudinal spin relaxation of divacancy qubits in SiC. The
T${}_{\\_}1$ time was reported for specific values of the magnetic field at
various temperatures [45, 46, 47], however, no systematic studies on the
environmental couplings and resulting spin relaxation processes have been
carried out yet. Even less is known about environmental resonances that may
give rise to T${}_{\\_}1$ spectroscopy and dynamic nuclear polarization in
SiC.

In this article, we study the longitudinal spin relaxation of divacancy qubits
in 4H-SiC due to various environmental spins, such as spin-1/2, spin-1, and
spin-3/2 point defects and 13C and 29Si nuclear spins. First, we identify the
most relevant resonance magnetic field values where efficient coupling to the
environmental spins leads to enhanced spin relaxation, and second, we
quantitatively study the spin relaxation rate, T${}_{\\_}1^{-1}$, over a wide
range of magnetic field values, with two orders of magnitude variation of the
defect concentration, and various isotope abundances. We provide a simple
analytical formula with theoretically fitted constants that can be used to
obtain ensemble averaged T${}_{\\_}1$ time in a sample of given defect
concentration and isotope abundance.

## Background

###  0.1 Theory of longitudinal spin relaxation

For a given quantisation axis, specified for example by the direction of the
external magnetic field, the magnetization of an integer or half-integer spin
can be described by the longitudinal and the transversal components. The
longitudinal, field parallel, magnetization is time independent for isolated
spins, while the transversal, perpendicular to magnetic field, magnetization
component rotates around the quantization axis with the Larmor frequency
$\omega=-\gamma B$, where $\gamma$ is the gyromagnetic ratio and $B$ is the
external magnetic field.

Due to the interaction of the spin with its environment, the above dynamics is
altered in two ways, the longitudinal magnetization becomes time dependent and
the Larmor frequency may be perturbed, i.e. the spin accumulates a phase shift
compared to the isolated case during its precession around the magnetic field.
While the total magnetization of a single spin is always preserved, ensemble
measurement and time averaged quantities exhibit the relaxation of both the
longitudinal and transversal magnetization due to the interaction with the
environmental. For example, the longitudinal magnetization of an ensemble
vanishes with time and approaches the thermal state with zero magnetization,
where the spins are pointing in any direction with equal probability. In a
large enough ensemble, the longitudinal magnetization relaxes exponentially to
the thermal state with relaxation time T${}_{\\_}1$. Similarly, the
accumulation of varying phase shifts either in an ensemble or over time gives
rise to the reduction of the net transversal magnetization, often referred to
as dephasing or decoherence.

In order to quantitatively describe spin relaxation and decoherence of spin a
ensemble in an open quantum system, one may use the density operator that can
be expressed as

$\hat{\varrho}=\sum_{\\_}n^{N}\left|\phi_{\\_}n\right\rangle\\!\\!\left\langle\phi_{\\_}n\right|\text{,}$
(1)

where $\phi_{\\_}n$ the wavefunction of spin $n$ in an ensemble of $N$ spins.
The density matrix of spin-1/2 objects can be written in a given quantization
basis as a $2\times 2$ matrix

$\varrho=\begin{pmatrix}\varrho_{\\_}{\uparrow\uparrow}&\varrho_{\\_}{\uparrow\downarrow}\\\
\varrho_{\\_}{\downarrow\uparrow}&\varrho_{\\_}{\downarrow\downarrow}\end{pmatrix}\text{,}$
(2)

where $\varrho_{\\_}{ij}=\left\langle
i\right|\hat{\varrho}\left|j\right\rangle=\sum_{\\_}n\left\langle
i\middle|\phi_{\\_}n\right\rangle\\!\\!\left\langle\phi_{\\_}n\middle|j\right\rangle$,
where $\left|i\right\rangle$ and $\left|j\right\rangle$ are the eigenstates of
the spin $z$ operator $S_{\\_}{z}$ and
$i,j=\left\\{\uparrow,\downarrow\right\\}$. The diagonal elements of the
matrix specify the net longitudinal magnetization of the spin ensemble as
$\bar{\mu}_{\\_}{\parallel}=g_{\\_}e\mu_{\\_}B\left(\varrho_{\\_}{\uparrow\uparrow}-\varrho_{\\_}{\downarrow\downarrow}\right)/2$,
while the off-diagonal elements of the density matrix specifies the coherence
properties of the ensemble [48, 49, 50]. In case of a spin-1 ensemble the
density matrix can be expressed by a $3\times 3$ matrix whose diagonal
elements $\varrho_{\\_}{k}=\varrho_{\\_}{kk}$, where
$k=\left\\{-1,0,+1\right\\}$, define the longitudinal magnetization.

In this paper we focus only on the longitudinal spin relaxation of point
defect spins. One may find further information on the theory and calculation
of the decoherence in Refs. [[48, 49, 50]]. In semiconductors, there are two
main contributions to the longitudinal spin relaxation dynamics of point
defect qubits. Spin-orbit coupling and phonon scattering assisted processes
are dominating the temperature dependence and are the most efficient
relaxation processes at high temperature in general. Due to this reason the
phenomena of longitudinal spin relaxation is frequently referred to as spin-
lattice relaxation in the magnetic resonance community.

In general, there are several contributions to the spin relaxation rate[51,
52] ($1/T_{\\_}1$) that can be written as

$\frac{1}{T_{\\_}1}=\frac{1}{T_{\\_}{1}^{dd}}+A_{\\_}{\text{dir,0}}+A_{\\_}{\text{dir,1}}\mathcal{T}+A_{\\_}{\text{Ram}}\mathcal{T}^{n}+\frac{A_{\\_}{\text{Orb}}}{e^{\Delta/k\mathcal{T}}-1}+A_{\\_}{\text{loc}}\frac{e^{\Delta_{\\_}{\text{loc}}/k\mathcal{T}}}{\left(e^{\Delta_{\\_}{\text{loc}}/k\mathcal{T}}-1\right)^{2}}+A_{\\_}{\text{therm}}\\!\left(\mathcal{T}\right)\text{,}$
(3)

where $1/T_{\\_}{1}^{dd}$ is the temperature independent spin-spin or dipolar
relaxation rate, the main subject of the paper, $A_{\\_}{\text{dir,0}}$ and
$A_{\\_}{\text{dir,1}}$ are the coefficient of the temperature independent and
the temperature dependent direct one phonon relaxation process, $\mathcal{T}$
is the temperature, $A_{\\_}{\text{Ram}}$ is the coefficient of the Raman
process, $n=\left\\{5,6,7\right\\}$ is the exponent of the temperature
dependence of the Raman process[51, 53, 54, 52], $A_{\\_}{\text{Orb}}$ is the
coefficient of the Orbach process, $\Delta$ is the energy separation between
the ground state and the excited state involved in the Orbach process, $k$ is
the Boltzmann constant, $A_{\\_}{\text{loc}}$ is the coefficient of the
contribution from local-vibrational mode scattering,
$\Delta_{\\_}{\text{loc}}$ is the energy of the local vibrational mode, and
$A_{\\_}{\text{therm}}$ describes thermally activated processes. In wide band
gap semiconductors, such as diamond and SiC the first four terms on the right
hand side of Eq. (3) are the most relevant ones.[55, 38, 39, 56, 57, 58]
Temperature dependence of the spin relaxation rate of divacancy qubits in SiC
has been reported very recently.[58] We fit Eq. (3) to the data in Ref. [[58]]
to quantify the strength of the various relaxation processes. The fitted
parameters for the dipolar spin relaxation, direct process, and Raman process
can be found in Table 1. Orbach, local, and temperature activated relaxation
processes do not play a significant role in the reported temperature range for
divacancy qubits. We note that the temperature dependence in the 200-300 K
interval is best fitted by $\mathcal{T}^{7}$ Raman processes[52]. Indeed, when
$n$ is also varied in the least square fit, we obtain 6.953 for the exponent
in the sample implanted with 1014 cm-2 N${}_{\\_}2$ ion dose. Note also that
the temperature independent spin relaxation contribution plays crutial role,
and its value depends only slightly on the implantation dose. From the fit we
cannot determine the weights of the dipolar spin relaxation and the
temperature independent direct one phonon relaxation processes. The latter
contribution is however sample independent. Further, discussion on this
observation in the light of our results can be found in the Discussion
section.

Table 1: Parameters of Eq. (3) for divacancy spin in 4H-SiC implanted with different dose of nitrogen molecule ion in Ref. [[58]]. $A_{\\_}{\text{Orb}}$, $A_{\\_}{\text{loc}}$, and $A_{\\_}{\text{therm}}$ are set to zero as the corresponding processes do not play a role.  N${}_{\\_}2$ ion dose (cm-2) |  $\left(\frac{1}{T_{\\_}{1}^{dd}}+A_{\\_}{\text{dir,0}}\right)$ (ms-1) |  $A_{\\_}{\text{dir}}$ (ms-1K-1) |  $A_{\\_}{\text{Ram}}$ (ms-1K-7) | $n$
---|---|---|---|---
1014 | 1.765 | 0.005866 | 2.715$\times 10^{-17}$ | 7
1013 | 1.575 | 0.004980 | 1.916$\times 10^{-17}$ | 7

Recent advances in first principles theory of spin-phonon relaxation processes
have made parameter free calculations of temperature dependent spin relaxation
rate possible, see Refs. [[59, 60, 61, 62, 63]] for further details. The first
principles simulation of dipolar spin relaxation of point defect spins has
also recently been addressed in Ref. [[64]]. This relaxation process may be
critical for shallow defects in nanometer scale sensing applications and for
integrated devices, where the dipolar spin relaxation may enhance due to the
increased doping level. Coupling to defects with various local spin
Hamiltonians gives rise to a rich phenomena with numerous resonances where the
spin relaxation rate can be enhanced by several orders of magnitude.[39, 64,
44] In this paper we utilize this latter method to quantitatively study the
temperature independent longitudinal spin relaxation rate of divacancy spin
qubits in various spin environments. In the following few paragraphs we
shortly review the theory and basic features of this method.

Decoherence[48, 49] and dipolar spin relaxation are inherently many-body
phenomena as they originate from the coupling of a central spin to its
environment of countless other spins. As the size of the Hilbert space
exponentially increases with the number of environmental spins, exact solution
of a sufficiently large model is unfeasible. Therefore, approximations must be
used. Furthermore, for point defect qubits the time scale of the longitudinal
relaxation is often $10^{6}$-$10^{8}$ longer that the characteristic time
scale of the coherent oscillations, thus application of state-of-the-art
quantum dynamics simulation approaches, such as tDMRG, are impractical.

Figure 1:  Illustration of the approximations used in the calculation of
dipolar spin relaxation (a) Many-spin system. (b) Cluster approximation and
central spin approximation of the many-spin system. (c) Coupled Lindbladian
master equation to retain relevant properties of the modelled system.

In order to address this computationally highly demanding problem, various
cluster expansion[48, 49] and cluster approximation models have been
developed. In the calculation of the coherence function of a many-spin system,
the cluster correlation expansion (CCE) method has gained considerable
popularity in the recent years due to its computational efficiency and
accuracy.[48, 49] On the other hand, to quantitatively include spin relaxation
effects one needs to go even beyond the standard CCE method.[50] For a
qualitatively accurate description of spin relaxation a modified methodology
is needed, which was developed in Ref. [[64]]. In this theory, a cluster
approximation is utilized to reduce the Hilbert space of the many-spin system.
The central spin and a few environmental spins are included in each cluster
system, see Fig. 1(b). In contrast to the CCE method[48, 49], the sub-clusters
of smaller order are not considered here. The ensemble of these small cluster
systems are meant to model the dynamics of the many-spin system. Due to the
truncation of the Hilbert space, correlation among the environmental spins are
neglected. This approximation is valid either when the coherence time of the
bath spins is short or when the coupling between spins are weak, such that the
intra-bath correlation cannot have an effect on the dynamics of the system in
reasonable time scales. The former is valid, for example, when the qubit is
coupled to a bath of electron spin defects with short coherence time, while
the latter is valid for a weakly coupled nuclear spin bath. Note, that by
increasing the size of the cluster systems, i.e. including an increasing
number of environmental spins, the exact behaviour may be asymptotically
recovered.

As can be seen, the cluster approximation naturally fits to the problem of
dipolar spin relaxation, however, the dynamics of a system of independent
cluster systems cannot properly describe the diagonal elements of the central
spin’s density matrix in the many spin system. The most obvious reason is that
the total magnetization of the cluster of independent smaller systems is not
preserved. Furthermore, since the central spin is included in all the cluster
systems, multiples of its magnetization give rise to a false thermal state. In
order to overcome these limitations, the dynamics of the cluster systems are
bound together in Ref. [[64]] through a coupled master equation, see Fig.
1(c), which can be written in the simplest case as

$\varrho_{\\_}{i}=-\frac{i}{\hbar}\left[\mathcal{H}_{\\_}i,\varrho_{\\_}i\right]+\mathcal{L}\\!\left(\varrho_{\\_}i\right)\text{,}$
(4)

where $\varrho_{\\_}{i}$ is the density matrix of cluster system $i$,
$\mathcal{H}_{\\_}i=H_{\\_}d+H_{\\_}i+H_{\\_}{di}$ is the Hamiltonian of the
cluster system including the Hamiltonian of the defect spin ($H_{\\_}0$) and
one coupled environmental spin ($H_{\\_}i$) and their interaction
($H_{\\_}{di}$), and $\mathcal{L}\\!\left(\varrho_{\\_}i\right)$ is an
extended Lindbladian that realizes the coupling between the cluster systems.
The last term on the right hand side of Eq. (4) can be expressed as

$\mathcal{L}\\!\left(\varrho_{\\_}i\right)=\sum_{\\_}{j=1,j\neq
i}^{N}\sum_{\\_}l\frac{\Delta
a_{\\_}{jl}}{\text{Tr}\\!\left(C_{\\_}{l}^{\dagger}C_{\\_}{l}\varrho_{\\_}i\right)}\left(C_{\\_}l\varrho_{\\_}iC_{\\_}l^{\dagger}-\frac{1}{2}\left(\varrho_{\\_}iC_{\\_}l^{\dagger}C_{\\_}l+C_{\\_}l^{\dagger}C_{\\_}l\varrho_{\\_}i\right)\right)\text{,}$
(5)

where $C_{\\_}l$ form a complete set of jump operators that can realize any
population transfer between the diagonal elements of the reduced density
matrix of the central spin and $\Delta a_{\\_}{jl}$ are the time dependent
rates that measure spin flip-flops caused by environmental spin $j$ in the
corresponding cluster system.

The Lindbladian coupling is defined so that the diagonal elements of the
reduced density matrix of the central spin are identical in all cluster
systems and the global magnetization is preserved at any time during the
simulation, see Ref. [[64]] for further details. This way the decaying
dynamics of the coupled cluster system approaches a proper thermal state.
Relaxation of the central spin originates from the coupling to the spin bath
and the decoherence of the bath.[64] Decoherence properties of the bath thus
have a quantitative effect on the relaxation properties. Consequently, the
method is most accurate when the cluster approximation, and resultant
reduction of the Hilbert space, is suitable for the considered system. As we
discussed before, cluster approximation expectedly holds for electron spin
environment of short coherence time and nuclear spin environments of weak
couplings. When a spin bath includes more than one strongly coupled nuclear
spins, inclusion of intra-spin-bath correlation is required, which can be
achieved by increasing the number of bath spins within the cluster systems.

###  0.2 Silicon carbide, the divacancy qubits, and relevant spin defects

Silicon carbide is a polytypic material with the 4H polytype being the most
stable[65] and most commonly used. The primitive cell of 4H-SiC consists of
four Si-C double layers in hexagonal stacking order. Due to the stacking of
the double layers, there are two non-equivalent positions for single site
point defects in 4H-SiC. The next nearest neighbor sites of the so-called $h$
and $k$ defect configurations show hexagonal-like and cubic-like local
arrangements, respectively. In a perfect 4H-SiC lattice, the divacancy, as an
adjacent pair of vacancies, has four distinguishable configurations $hh$,
$kk$, $hk$, and $kh$, where the first and second letter indicate the carbon
and silicon vacancy sites, respectively. These configurations in the above
specified order were assigned[66] to the PL1, PL2, PL3, and PL4
photoluminescence (PL) signals[45]. Here, we note that there are additional
divacancy related PL lines observed in experiment that were named as
PL5-PL6[45, 46] as well as PL5’-PL6’ and PL7-PL10[67]. Recently, some of these
lines were tentatively assigned to stacking fault-divacancy complexes.[68] In
this article, we consider the four regular divacancy configurations (PL1-4)
only.

There are several spin defects in 4H-SiC that can influence the spin dynamics
of a point defect qubit. Considering the total spin momentum and the kind of
the spins, we identify four groups of spin defects that we consider separately
in this article: 1) Spin-1/2 nuclear spins due to paramagnetic 13C (1.07%) and
29Si (4.68%) isotopes, which are naturally abundant in SiC. 2) Dopants and
intrinsic defects with a spin-1/2 electronic ground state are the dominant
electron spin defects in SiC. Besides these two common spin defects, we
consider two, more exotic categories, namely, 3) nearby divacancies as
environmental spin defects, and 4) spin-3/2 silicon vacancies as environmental
spin defects. While these defects have generally low concentrations (
$<10^{14}$ cm-3) in high purity samples[69], due to irradiation and ion
implantation, often used to create divacancy and silicon vacancy qubits, their
concentration may be considerably enhanced either throughout the sample or
locally. For example, after irradiation and annealing, the silicon vacancy
concentration can be in the $10^{16}$ cm-3 range or even more in heavily
irradiated samples.[70]

In the second group of paramagnetic defects, there are several spin-1/2
dopants and defects that can be present in SiC with varying concentrations
depending on the growth process and after-growth treatments. In as-grown high
purity SiC epilayers, the dominating impurities and defects with non-zero
electron spin are: a) nitrogen substitutional defects with typical
concentration of $10^{14}$-$10^{15}$ cm-3 for commercial layers, $10^{13}$
cm-3 for ultra-pure epilayers[69], $10^{14}$-$10^{15}$ cm-3 for high-purity
semi-insulating (HPSI) 4H-SiC.[71] It is worth mentioning that the nitrogen
nucleus also possesses a spin of $I=1$ (99.63% natural abundance) that couples
through an isotropic hyperfine interaction of 50.97 MHz to the electron
spin.[72] b) Shallow boron defect with $10^{14}$ cm-3 concentration for
commercial layers, $10^{13}$ cm-3 concentration for pure layers, and
$10^{14}$-$10^{15}$ cm-3 for HPSI 4H-SiC.[71] While boron also possesses a
none zero nuclear spin, the spin density of the defect is localized on other
neighboring atoms with low paramagnetic isotope abundance[72], thus
hereinafter we neglect the hyperfine interaction for the shallow boron
acceptor. c) Carbon vacancy with $10^{13}$ cm-3 in pure epilayers and
$10^{15}$ cm-3 in HPSI 4H-SiC.[69] d) Carbon antisite-vacancy (CAV) defects
can be present in HPSI materials with concentration $10^{14}$-$10^{15}$
cm-3.[71] e) In p-type samples the concentration of aluminum is typically in
the range of 1-5$\times 10^{18}$ cm-3. The hyperfine splitting due to
paramagnetic aluminum isotope is not resolved in 4H-SiC.[72]

## Results

###  0.3 Resonant spin bath couplings

Figure 2:  a) Magnetic field dependence of the energy levels of $hh$ (dark
blue) and $hk$ (green) divacancies coupled to a single spin-1/2 point defect.
b) Population of the $m_{\\_}S=\pm 1$ spin states of the $hh$ (dark blue),
$kk$ (light blue), $hk$ (green), and $kh$ (amber) divacancies induced by a
bath of spin-1/2 environmental defects at the concentration of 1018 cm-3. The
depicted curves show the magnetic field dependence of the population
enhancement of the $m_{\\_}S=\pm 1$ states of the divacancies. Each curve
exhibits two resonances corresponding to
$B_{\\_}{\text{GSLAC}}\\!\left(D,E\right)$ (group I) and
$B_{\\_}{\text{GSLAC}}\\!\left(D,E\right)/2$ (group II), where $D$ and $E$ are
the ZFS parameters of the corresponding divacancy qubit. In a) and b) the
magnetic field is parallel to the quantization axis of the divacancy qubits in
all cases. At $t=0$, the central divacancy is initialized in the $m_{\\_}S=0$
state, while the spin bath in an high temperature thermal state.

First we study resonance effects due to the most common spin defects in SiC,
i.e. nuclear spin and spin-1/2 point defects. Spin mixing curves due to 29Si
and 13C nuclear spins exhibit a single resonance place that corresponds to
ground state level anticrossing (GSLAC) of the divacancy electron spin states.
At the magnetic field value of the GSLAC,
$B_{\\_}{\text{GSLAC}}\\!\left(D,E\right)=\left(D-E\right)/g_{\\_}e\mu_{\\_}B$
the Zeeman shift of the $m_{\\_}S=-1$ divacancy spin states compensates the
zero-field-splitting. The gapless electron spin states can be efficiently
mixed either with each other through the precession of the electron spins due
to transverse magnetic field, or with nuclear spin states due to the hyperfine
interaction. The latter effect gives rise to strong relaxation effect at
$B_{\\_}{\text{GSLAC}}$. Note that for different divacancy configurations the
GSLAC resonance appears at different magnetic field values due to the
different $D$ and $E$ values of the centers. Note also that the different
divacancy configurations possess different neighboring shells, which, however,
on average has minor effects on the overall strength of the mixing.

Figure 3:  Population of the $m_{\\_}S=\pm 1$ spin states for the case of
spin-1/2 defect with strong hyperfine interaction. Solid light blue line shows
the population variation of the $hh$ divacancy spin states due to
substitutional nitrogen shallow donor. Dotted dark blue depicts population
variation due to spin-1/2 electron spin defects without the hyperfine
coupling. Hyperfine interaction leads to a splitting of the center peak into
three peaks at $B_{\\_}{\text{GSLAC}}/2$ as well as to an appearance of two
satellite peaks next to the GSLAC peaks.

Spin-1/2 point defects, like boron, aluminum, nitrogen dopants, carbon
vacancies, and CAV defects give rise to two distinct level crossings at
$B_{\\_}{\text{GSLAC}}$ and $B_{\\_}{\text{GSLAC}}/2$, see Fig. 2(a), that
result in resonant spin relaxation, see group I and II peaks in Fig. 2(b). The
GSLAC resonance is due to the precession of the divacancy spins induced by the
transverse field of the spin defects. At $B_{\\_}{\text{GSLAC}}/2$, the
$\left|m_{\\_}S,m_{\\_}{\mathcal{S},i}\right\rangle=\left|0,+1/2\right\rangle$
and $\left|-1,-1/2\right\rangle$ states cross. Due to the dipole-dipole
interaction of the neighboring paramagnetic defects, a LAC is formed at
$B_{\\_}{\text{GSLAC}}/2$ and efficient spin flip-flop interaction takes
place. While divacancies can polarize nuclear spin at $B_{\\_}{\text{GSLAC}}$,
they can polarize spin-1/2 electron spins only at $B_{\\_}{\text{GSLAC}}/2$.
Note that when a strong hyperfine coupling is not expected at the spin defect
site, a featureless resonance peak is observed as seen in Fig. 2(b). On the
other hand, the nitrogen dopant, exhibits a strong hyperfine interaction with
its abundant 15N nuclear spin, resulting in a fine structure, see the dashed
line in Fig. 2 for the $hh$ divacancy. The hyperfine interaction at the spin
defect site leads to a splitting of the resonance peaks in group II at
$B_{\\_}{\text{GSLAC}}/2$ according to the nuclear spin states. In addition,
two satellite peaks appear on the GSLAC resonance peak, see Fig. 3. The curves
are similar to those reported for the NV center in diamond interacting with
substitutional nitrogen point defects [41, 44]. However, in hexagonal SiC
polytypes there is only one allowed symmetry distortion and the hyperfine
interaction is dominantly isotropic[72], which leads to a less complicated
spin structures at $B_{\\_}{\text{GSLAC}}/2$.

Figure 4:  Magnetic field dependence of the population of the divacancy
$m_{\\_}S=\pm 1$ spin states due to spin flip-flops induced by V1 and V2
silicon vacancy centers. a) and b) depict the magnetic field dependence of the
spin mixing of $hh$ and $kk$ and $hk$ and $kh$ divacancy centers,
respectively. In addition to the primary peaks I.a and II.a, observed for
spin-1/2 defect, satellite peaks I.b, I.c, II.b, and II.c can be observed to
the left of their respective primary peaks. The satellite peaks correspond to
double and triple quantum jumps of environmental silicon vacancies. In a) and
b) the magnetic field is parallel to the quantization axis of the divacancy
qubits in all cases. At $t=0$, the central divacancy is initialized in the
$m_{\\_}S=0$ state, while the spin bath in an high temperature thermal state.
Figure 5:  Identification of resonance peaks for the $hh$ divacancy
interacting with spin-3/2 silicon vacancy centers. a) and c) show the energy
levels structure of $hh$ divacancy-V1 silicon vacancy and $hh$ divacancy-V2
silicon vacancy two spin systems, respectively. b) depicts the derived
population curve for $hh$ divacancy. At each resonance peak, a pair of numbers
in parenthesis provide the jump of the spin quantum numbers of the divacancy
and silicon vacancy centers at the assigned level crossings.

More exotic couplings are possible when the divacancies interact with other
divacancies or with spin-3/2 silicon vacancy centers. This may be the case in
highly irradiated samples or when the sample is bombarded with positive ions
creating hundreds of defects in a small volume[73]. In Fig. 4, we depict the
magnetic field dependence of the population of the divacancy $m_{\\_}S=\pm 1$
spin states due to spin flip-flops induced by V1 and V2 silicon vacancy
centers[74]. The curves resemble in their main character the mixing curves of
the spin-1/2 defects. It is indeed expected as the quartet spin of the silicon
vacancy contains a $m_{\\_}{\mathcal{S}}=\left\\{+1/2,-1/2\right\\}$ manifold
that can interact with the divacancy spin states as the spin-1/2 defect spins.
Accordingly, we observe two main resonance I.a and II.a at magnetic field
values $B_{\\_}{\text{GSLAC}}$ and $B_{\\_}{\text{GSLAC}}/2$. In addition to
these dominant peaks several satellite peaks can be observed on the left side
of both primary resonances. Majority of the peaks can be connected to level
crossings, see Fig. 5. To allow spin mixing at these crossings the quantum
number of the silicon vacancy spin must jump either by two or by three, see
Fig. 5. Note that dipole-dipole interaction does not include such Hamiltonian
terms. However, higher order Zeeman terms are non-zero for quartet spin states
in C${}_{\\_}{3v}$ symmetry, see in Eq. (14), which in combination with the
dipole-dipole interaction makes a weak mixing of the divacancy and silicon
vacancy spin states possible. Peaks II.d require the triple quantum jumps,
i.e. ($\Delta m_{\\_}S,\Delta m_{\\_}\mathcal{S}$) = (1,3), while peaks I.b,
II.b, and II.c require double quantum jumps from the silicon vacancy spin,
i.e. ($\Delta m_{\\_}S,\Delta m_{\\_}\mathcal{S}$) = (1,2).

Figure 6:  Population of the $m_{\\_}S=\pm 1$ spin states of divacancies
interact with other divacancy spins. a) and b) show the mixing curves on a
logarithmic scale of the $hh$ and $kk$ and the $hk$ and $kh$ divacancy
centers, respectively. In all cases the spin bath contains all the possible
divacancy configurations other than the central divacancy under consideration.
In a) and b) the magnetic field is parallel to the quantization axis of the
divacancy qubits in all cases. At $t=0$, the central divacancy is initialized
in the $m_{\\_}S=0$ state, while the spin bath in an high temperature thermal
state.  Figure 7:  Identification of resonance peaks for the $hh$ divacancy
interacting with other divacancy centers. a) and c) show the energy levels
structure of $hh$ divacancy-$hk$ divacancy and $hh$ divacancy-$kh$ divacancy
two spin systems, respectively. At zero magnetic field the spin states split
into three well separated manifolds. b) depicts the derived spin population
curve for the $hh$ divacancy.

In order to investigate the interaction between divacancy centers we studied
each divacancy configuration interacting with a bath of other possible
divacancies. The obtained spin polarization curves exhibit a complicated
pattern which has not been seen so far, see Fig. 6. Besides the GSLAC
resonance, C${}_{\\_}{3v}$ symmetric configurations exhibit two other groups
of resonance peaks, peaks III and IV. In order to understand the origin of
these resonances, we shortly discuss the energy level structure of divacancy
pairs. At zero magnetic field, the states of all divacancy pairs, irrespective
to their local quantization axis, split into three manifolds, 1) the single
degenerate
$\left|m_{\\_}S,m_{\\_}{\mathcal{S}}\right\rangle=\left|0,0\right\rangle$
manifold with the lowest energy, 2) the $\left\\{\left|0,\pm
1\right\rangle,\left|\pm 1,0\right\rangle\right\\}$ manifold, and C) the
$\left\\{\left|\pm 1,\pm 1\right\rangle,\left|\mp 1,\pm
1\right\rangle\right\\}$ manifold with the highest energy, see Fig. 7. The
GSLAC resonance is due to the crossing between the $\left|0,0\right\rangle$
state and states of the $\left\\{\left|0,\pm 1\right\rangle,\left|\pm
1,0\right\rangle\right\\}$ manifold. Group III peaks appear at crossings of
manifolds 2) and 3), see Fig. 7. Note that spin mixing at such level crossings
is only possible when divacancies with two different quantization axes are
interacting. In this case the Zeeman interaction partially mixes the spin
states of one of the divacancies. Combined with dipole-dipole interaction this
makes spin flip-flop interaction possible at the crossings. Accordingly, peaks
III.a and III.c and peak III.b in Fig. 6(a) originate from interactions with
$hk$ and $kh$ basal plane oriented divacancy centers, respectively, see Fig.
7. Group III peaks include four peaks for basal plane oriented divacancy
configurations (Fig. 6(b)), each of them corresponds to interaction with
divacancies, $hh$, $kk$, $hk$, and $kh$, whose quantization axis is 109 degree
aligned to the quantization axis of the considered basal plane divacancy.
Group IV. peaks correspond to crossings between the states of the
$\left\\{\left|0,\pm 1\right\rangle,\left|\pm 1,0\right\rangle\right\\}$
manifold, see Fig. 7. As the splitting of the states in the middle manifold is
small compared with the $D$ zero-field-splitting parameter, the resonances
peaks appear at small magnetic field values in most cases. The positions of
the crossings largely depended on the differences of the $D$ splittings and
the values of the $E$ splitting of the coupled divacancy centers. The ZFS
parameters of the $hk$ divacancy deviates the most from the ZFS parameters of
the other configurations. Consequently, the crossing region stretches over a
magnetic field interval of 175 G. Furthermore, accidental degeneracies can be
observed in the energy level structure (not shown) that give rise to broad
resonance peaks as can be seen in Fig. 6(b) for the $hk$ divacancy center.

As the different spin defects give rise to different signatures, study of spin
mixing can be utilized to identify the dominant spin-spin interactions in a
sample. In case of single defects, spin mixing, detectable through the spin
dependent luminescence of the centers, can provide valuable information on the
local environment of the centers.

###  0.4 Longitudinal spin relaxation

Figure 8:  Magnetic field dependence of the longitudinal spin relaxation rate
T${}_{\\_}1^{-1}$ due to a) nuclear spin of natural abundance and b) spin-1/2
point defects of various concentrations. a) Relaxation due to 29Si (solid dark
blue line) and 13C (solid light blue line) nuclear spins are considered
separately. When 13C spin bath is considered the nearest neighbor sites are
excluded. For 29Si we also depict the largest (dashed dark blue line) and the
smallest (dotted dark blue line) relaxation rate found in our ensemble of 100
configurations. The shaded are between these curves covers more than four
orders of magnitude.

In this section, we carry out quantitative study on the relaxation dynamics of
the divacancy qubits in different spin environments and determine the magnetic
field and concentration dependence the longitudinal spin relaxation rate
T${}_{\\_}1^{-1}$. Here, we focus only on the most relevant spin defects in
SiC, namely 29Si and 13C nuclear spins and spin-1/2 point defects. Figs. 8 a)
and b) depict the calculated spin relaxation rate for $hh$ divacancy qubit
embedded in a nuclear spin environment and spin-1/2 point defect environment
of different concentrations, respectively. Here, we note that when 13C spin
bath is considered the nearest neighbor sites are excluded. Comparing the 13C
nuclear spin environment induced spin relaxation rate with the one obtained
for the NV center[64], one can notice that the paramagnetic carbon nuclei of
natural abundance give rise to comparable but somewhat smaller relaxation rate
for the $hh$ divacancy in SiC. This is due to the fact that the lattice
constant of SiC is larger than in diamond, thus the hyperfine interaction
decreases faster within the neighboring shells. In addition, carbon atoms
occupy only half of the sites in the vicinity of the defect. For example, the
second neighborhood shell, that contains twelve sites with hyperfine
interaction $\approx$10 MHz, is occupied by silicon atoms. This very last
argument suggests that 29Si paramagnetic silicon nuclei of 4.6% natural
abundance have strong effect on the spin relaxation time. Indeed, according to
our simulations 29Si nuclear spins limit the ensemble averaged T${}_{\\_}1$
time to 100-500 ms farther away from the GSLAC resonance. We note that the
spin relaxation rate of individual divacancy centers can vary over several
orders of magnitude due to their distinct local environment. As can be seen in
Fig. 8, the largest relaxation rate found in our ensemble is an order of
magnitude larger than the average, while the smallest one is more than three
orders of magnitude smaller than the average.

Comparing Figs. 8(a) and (b), one can see that, spin-1/2 point defect induced
relaxation starts to dominate, except the resonance place at
$B_{\\_}{\text{GSLAC}}/2$, only at defect concentrations as high as $3\times
10^{18}$ cm-3. Therefore, we conclude that in samples of natural isotope
abundance, 29Si nuclei are the main non-thermal source of spin relaxation. In
29Si isotope depleted samples, depending on their net concentration, spin-1/2
point defects may dominantly contribute to spin relaxation at low temperature.
Fig. 8(b) shows the magnetic field dependence of the relaxation rate at
various concretions of spin-1/2 point defects ranging from 1016 cm-3 to 1018
cm-3. Note that the lowest concentration considered in this study is still
considered to be high in SiC, however, below this level 13C nuclear spins
start to dominate the relaxation dynamics.

Figure 9:  Concentration dependence of the fitting parameters of Eq. (6)
fitted to the simulated spin-1/2 defect concentration dependent spin
relaxation rate curves in Fig 8(b).

In order to make our numerical results more accessible, we parameterize all
the curves depicted in Fig. 8(b). A Lorentzian curve is adequate to describe
the effects of resonant spin state couplings in the few spin limit. Due to the
large number of spins included in the bath and the vast variation of spin bath
configurations in an ensemble, the actual spin relaxation curve at the
resonance places can be considered as a superposition of numerous Lorentzian
curves of varying width. The distribution of the width of the Lorentzian
curves is dictated by the distribution of the coupling strength between the
defect and the spins in the environment. For simplicity, we fit a single
Lorentzian at each resonance peak that can describe the decay of the
relaxation rate properly farther away from the resonance place, but may be
inaccurate right at the place of the level crossing. The relaxation rate is
thus parameterized as

$T_{\\_}1^{-1}\\!\left(B,C\right)=T_{\\_}1^{-1}\\!\left(0,C\right)+\frac{2A_{\\_}1}{\pi}\frac{\omega_{\\_}1}{4\left(B-B_{\\_}{\text{GSLAC}/2}\right)^{2}+\omega_{\\_}1^{2}}+\frac{2A_{\\_}2}{\pi}\frac{\omega_{\\_}2}{4\left(B-B_{\\_}{\text{GSLAC}}\right)^{2}+\omega_{\\_}2^{2}}\text{,}$
(6)

where $B$ is the magnetic field, $C$ is the concentration of spin-1/2 defect,
$B_{\\_}{\text{GSLAC}/2}\approx 240$ G and $B_{\\_}{\text{GSLAC}}\approx 480$
G. $T_{\\_}1^{-1}\\!\left(0,C\right)$, the relaxation rate at $B=0$,
$A_{\\_}1$, $A_{\\_}2$, $\omega_{\\_}1$, and $\omega_{\\_}2$ are fitting
parameters which are depicted in Fig. 9 as a function of the spin-1/2 defects
concentration. In order to be able to extrapolate to concentrations that are
not considered in our numerical simulations, we fitted polynomial curves to
the concentration dependence of the parameters of Eq. (6), as

$X\\!\left(C\right)=\alpha+\beta C^{n}\text{,}$ (7)

where $X$ is a parameter of the double-Lorentzian curve. The simple polynomial
form given in Eq. (7) fits well onto the points depicted in Fig. 9. Table 2
provides the parameters of these polynomial fits. To obtain the relaxation
rate for a particular defect concentration given in cm-3 unit, one can
determine the corresponding Lorentzian parameters using Eq. (7) and Table 2
and insert them into Eq. (6). This will provide the ensemble averaged
relaxation rate in Hz unit.

Table 2: Parameters of the polynomial fits describing concentration dependence parameters of the double-Lorentzian curve fitted to the numerical results. Parameters of the Lorentzians | $\alpha$ | $\beta$ | n
---|---|---|---
A${}_{\\_}1$ | 0 | $1.75\times 10^{-25}$ | 1.65
A${}_{\\_}2$ | 0 | $1.23\times 10^{-29}$ | 1.9
$\omega_{\\_}1$ | 2.87 | $1.6\times 10^{-15}$ | 0.875
$\omega_{\\_}2$ | 6.89 | $1.23\times 10^{-17}$ | 1
T${}_{\\_}1^{-1}\\!\left(B=0\right)$ | 0 | $1.6\times 10^{-35}$ | 2

Considering the concentration dependence of the spin-1/2 point defects induced
relaxation rate at zero magnetic field, we obtain a simple quadratic
concentration dependence of the relaxation rate,

$T_{\\_}1^{-1}\left(B=0,C\right)=\beta C^{2}$ (8)

where $\beta=1.6\times 10^{-35}$ Hz/cm-6. Note that the quadratic dependence
given by Eq. (8) differs from the earlier experimental reports[39] observing
linear dependence on the concentration. We attribute the difference to the
uncertainties in the concentration and relaxation rate measurements and to the
smaller concentration interval considered in the experiment.

## Discussion

As we have seen in Eq. (3), there are two main contributions to the
temperature independent part of the spin relaxation, namely the dipolar spin
relaxation, due to the surrounding spin bath, and the temperature independent
contribution of the direct one phonon scattering process. Separating the
contributions of these processes is sometimes cumbersome and requires low
temperature measurements[61]. As can be seen from the results of Ref. [[58]]
and from our fit provided in Table 1, direct relaxation plays an important
role for divacancy defects, thus its temperature independent contribution may
not be ruled out. On the other hand, this process is sample independent, in
contrast to the dipolar spin relaxation, as it depends only on the intrinsic
properties of the SiC. Therefore, the lowest ever reported relaxation rate of
divacancy spins in SiC sets an upper bound for the related constant
$A_{\\_}{\text{dir,0}}$. So far the longest relaxation time of divacancy
defects is reported in 6H-SiC and found to be in the 10 ms range at 20 K.[46]
Note that 4H and 6H-SiC polytypes share common properties and we expect that
the $A_{\\_}{\text{dir,0}}$ parameter is in the same order of magnitude in
4H-SiC. This assumption is further supported by the results of Ref. [[24]],
where the reported 64 ms coherence time indicates a relaxation time in the
order of 10 ms in 4H-SiC. This means that the $A_{\\_}{\text{dir,0}}<0.1$ kHz,
which is an order of magnitude smaller than the fitted values in Table 1.
Consequently, the temperature independent relaxation rate of divacancy qubits
in the nitrogen ion irradiated sample in Ref. [[58]] is dominantly due to
dipolar spin relaxation processes and it is found in the range of 1.8 ms-1 and
1.6 ms-1 at 180 G external magnetic field for 1014 and 1013 cm-2 implantation
doses, respectively.

Our results depicted in Fig. 8 clearly show that the hyperfine interaction
cannot account for the 1.8 ms-1 and 1.6 ms-1 dipolar relaxation rate at 180 G,
it is three orders of magnitude smaller. On the other hand, paramagnetic
defects created by nitrogen ion irradiation can account for the observed
values. Utilizing the analytical fit provided in Eqs. (6) and (7) and Table 2,
we obtain $C=4.1\times 10^{18}$ cm-3 and $3.9\times 10^{18}$ cm-3 defect
concentration for the 1014 cm-2 and 1013 cm-2 N${}_{\\_}2$ ion irradiated
samples. Such a high defect contraction may be expected due to the large
irradiation dose applied in the experiment. We note, however, that the
concentration only slightly depends on the implantation dose. A possible
explanation of the implantation dose independence could be given based on the
fact that the T${}_{\\_}1$ time based concentration measurement provides
information only on the local environment of the defects and not on the
concentration of the entire sample. Since positive ions create a large damage
locally in an implantation layer, it is expected that the defect concentration
is highly inhomogeneous. Furthermore, different ions create different local
damage, thus the observed $C\approx 4\times 10^{18}$ cm-3 defect
concentrations may characterize the method of implantation (30 keV
N${}_{\\_}2$ ions) rather than the dose and the sample’s overall defect
concentration. Finally, we note also that the selection of the external
magnetic field value is quite unfortunate, since at this magnetic field the
divacancy defects can couple relatively strongly to both spin-1/2 defects as
well as the spin-3/2 silicon vacancy, see Fig. 2 and Fig. 5.

## Summary

In this article we have reported on a systematic study of the spin relaxation
dynamics of divacancy qubits in 4H-SiC. First, we examined the magnetic field
dependence of the spin mixing induced by various spins in the local
environment of divacancy qubits. We demonstrated that neighboring divacancy
centers and the spin-3/2 silicon vacancy centers give rise to rich relaxation
patterns with multiple of relaxation peaks. The different magnetic field
dependent spin relaxation patterns reported in this article may allow one to
identify and study the local environments of a single or an ensemble of
divacancy qubits by optical means. Related applications, such as dynamic
nuclear polarization, microwave free sensing, and spin relaxation spectroscopy
are yet to be explored in SiC.

Moreover, we have simulated the magnetic field and the concentration
dependences of the spin relaxation time T${}_{\\_}1$ for the most relevant
spin defects in SiC. We have shown that in high purity samples the dominant
non-thermal contribution to the spin relaxation comes from the 29Si nuclear
spin bath, that maximizes the ensemble averaged T${}_{\\_}1$ time at 100 ms at
low temperature far away from the GSLAC resonance. For configurations of
adjacent nuclear spins the relaxation time may reduce to 40 ms at zero
magnetic field. Note that the 29Si spin bath limited spin relaxation time is
comparable or even shorter than the coherence time of 64 ms reported in Ref.
[[24]] for a decoherence protected subspace of divacancy qubits. In such
cases, the longitudinal spin relaxation may be the major limiting factor for
the lifetime of the coherence protected subspace.

We also demonstrated that dipolar spin relaxation of divacancies due to
paramagnetic point defects may be significant in ion implanted samples. The
analytical formula provided in this paper can be utilized either to estimate
T${}_{\\_}1$ in a given sample of known spin defect concentration or to
analyze the measured T${}_{\\_}1$ time to estimate the local spin defect
concentration of divacnacy qubits. By using this latter approach, we showed
that the local concentration of paramagnetic point defect in N${}_{\\_}2$ ion
implanted samples can be as high as $4\times 10^{18}$ cm-3 that maximizes the
coherence time in about 0.5 ms. This value is a factor of two smaller than the
1.3 ms coherence time achievable in 4H-SiC of natural isotope abundance.[49]

## Methodology

In order to faithfully simulate spin relaxation processes, we utilize the
method developed in Ref. [[64]]. In the simulations we consider a central
divacancy defect interacting with its environment, i.e. a bath of spins
consisting of one type of spin defects from categories 1)-4), in a central
spin arrangement. Considering only one type of spins, typically the closest
spin defects play the major role in thermalizing the central spin. We found
that inclusion of the 32 closest defects already provides convergent results
in terms of spin relaxation. The many-spin system is then divided into a
cluster of subsystems which includes the central divacancy and one defect spin
each, see Fig. 1 in Ref. [[64]]. The applied theoretical method[64] introduces
an effective interaction between the subsystems in such a way that the total
spin of the many-spin system is preserved throughout the simulation. This
ensures reliable simulation of spin relaxation and calculation of
T${}_{\\_}1$, as demonstrated for the case of the NV center in Ref. [[64]]. In
this article, we study ensemble averaged quantities. Therefore, in all cases
an ensemble of randomly generated local spin environments is considered. In
our qualitative (quantitative) study of divacancy ensembles we consider 100
(200) random spin configurations corresponding to a given defect concentration
on average. The time step is set to 25 ps (20 ps), while the simulation time
is set to 1 $\mu$s (600 $\mu$s) for the qualitative (quantitative) analysis of
the relaxation processes. When the hyperfine interaction induced spin
relaxation is investigated, 100 ps time step and 600 $\mu$s simulation time
are used, and ensemble averaging is carried out over 100 random configurations
of a bath containing 128 nuclear spins. Increasing the number of spins is
motivated by the spatial extension of the spin density of the divacancies,
which is comparable with the distances of the closest spins alike to the case
of point defect spins.

In our real time simulations we investigated the time dependent variation of
the initial population of the divacancy spin states. At $t=0$, the central
divacancy is initialized in the $m_{\\_}S=0$ state, while the spin bath in an
high temperature thermal state.

We simulate the dynamics of the coupled spin system without any additional
approximation on the local Hamiltonian of the divacancy, the spin defects, and
the environmental coupling of the divacancy. Intra-spin-bath interactions, on
the other hand, are neglected. The spin Hamiltonian thus can be written in the
form of

$H=H_{\\_}{\text{div}}+H_{\\_}{\text{bath}}+H_{\\_}{\text{coup}}\text{,}$ (9)

where the Hamiltonian $H_{\\_}{\text{div}}$ of the divacancy consists of three
terms

$H_{\\_}{\text{div}}=g_{\\_}e\mu_{\\_}{\text{B}}B_{\\_}zS_{\\_}z+D\left(S_{\\_}z^{2}-\frac{2}{3}\right)+\frac{E}{2}\left(S_{\\_}+^{2}+S_{\\_}-^{2}\right)\text{,}$
(10)

where $g_{\\_}e$ is the electron g-factor, $\mu_{\\_}e$ is the Bohr magneton,
$B_{\\_}z$ is the external magnetic field parallel to the quantization axis of
the central divacancy defect set by the
V${}_{\\_}\text{C}$-V${}_{\\_}\text{Si}$ defect axis, $S_{\\_}z$ is the
electron spin $z$ operator, $S_{\\_}+$ and $S_{\\_}-$ are the electron spin
ladder opearators, and $D$ and $E$ are the zero-field splitting parameters.
Due symmetry constraints the $E$ splitting is zero for the C${}_{\\_}{3v}$
symmetric $hh$ and $kk$ configurations, while it is non-zero for the
C${}_{\\_}{1h}$ symmetric $kh$ and $hk$ configurations. The experimental zero-
field-splitting parameters used in the simulations are reported in Ref.
[[Falk2013]]. The bath Hamiltonian and the coupling Hamiltonian depend on the
spin defect considered. For nuclear spins ($I=1/2$) the bath Hamiltonian
includes the nuclear Zeeman term only

$H_{\\_}{\text{bath}}^{\text{nuc}}=-\sum_{\\_}ig_{\\_}{\text{n},i}\mu_{\\_}{\text{N}}B_{\\_}zI_{\\_}{z,i}\text{,}$
(11)

where $\mu_{\\_}{\text{N}}$ is the nuclear magneton, $g_{\\_}{\text{n},i}$ is
the nuclear g-factor of nucleus $i$ being either a 29Si nucleus or a 13C
nucleus, and $I_{\\_}{z,i}$ is the nuclear spin $z$ operator of nucleus $i$.
For a nuclear spin bath, the coupling Hamiltonian includes the hyperfine term

$H_{\\_}{\text{coup}}^{\text{nuc}}=\sum_{\\_}iSA_{\\_}iI_{\\_}i\text{,}$ (12)

where the hyperfine tensors $A_{\\_}i$ are obtained from _ab initio_ supercell
calculations combined with real space integration methods[64] in order to
eliminate finite size effects.

The bath Hamiltonian for spin-1/2 point defects can be written as

$H_{\\_}{\text{bath}}^{\text{spin-1/2}}=\sum_{\\_}ig_{\\_}e\mu_{\\_}{\text{B}}B_{\\_}z\mathcal{S}_{\\_}{z,i}+\mathcal{S}_{\\_}i\mathcal{A}_{\\_}i\mathcal{I}_{\\_}i\text{,}$
(13)

where the first term is the Zeeman energy of the spin-1/2 bath and the second
term is the leading hyperfine interaction for a considered paramagnetic
defect. Here, $\mathcal{S}$ and $\mathcal{I}$ are the electron spin operator
vector of the environmental defect spin and the nuclear spin operator vector
of a nucleus strongly coupled to the defect’s electron spin through the
hyperfine tensor $\mathcal{A}$, respectively. Based on the discussions in
introduction, we consider only the nuclear spin of the 14N of the shallow
nitrogen donor in our study.

The bath Hamiltonian of environmental divacancy spins contains the Zeeman term
and the zero-field-splitting terms as specified in Eq. (10). It is worth
mentioning that the quantization axis of the environmental divacancies may not
be parallel to the quantization axis of the central divacancy and the external
magnetic field. In the random spin bath configurations, we consider all
possible arrangements.

According to Ref. [[75]] the bath spin Hamiltonian for the quartet silicon
vacancy can be written as

$\displaystyle
H_{\\_}{\text{bath}}^{\text{spin-3/2}}=\sum_{\\_}j\left(D_{\\_}{\text{VS},j}\left(\mathcal{S}^{2}_{\\_}{zj}-\frac{5}{4}\right)+\left[g_{\\_}e\mathcal{S}_{\\_}{zj}+g_{\\_}{3\parallel}\frac{\mathcal{S}_{\\_}{+j}^{3}-\mathcal{S}_{\\_}{-j}^{3}}{4i}\right]\mu_{\\_}{\text{B}}B_{\\_}\parallel\right)+$
$\displaystyle\sum_{\\_}{j}\left(g_{\\_}e\mu_{\\_}{\text{B}}\mathcal{S}_{\\_}{\perp}B_{\\_}{\perp}+g_{\\_}{3\perp}\mu_{\\_}{\text{B}}\frac{\left\\{\mathcal{S}^{2}_{\\_}{+j},\mathcal{S}_{\\_}{zj}\right\\}B_{\\_}+-\left\\{\mathcal{S}_{\\_}{-j}^{2},\mathcal{S}_{\\_}{zj}^{2}\right\\}B_{\\_}-}{2i}\right)$
(14)

where the higher order g-factors $g_{\\_}{3\parallel}$ and $g_{\\_}{3\perp}$
are measured in Ref. [[75]], the zero-field-splitting constant
$D_{\\_}{\text{VS},j}$ of the environmental silicon vacancy centers is set
either to 2.6 MHz or to 35 MHz for the two possible silicon vacancy
configurations[74], $B_{\\_}\parallel$ and $B_{\\_}\perp$ are the component of
the magnetic field parallel and perpendicular to the C${}_{\\_}{3v}$ axis of
the silicon vacancy, and $B_{\\_}\pm=B_{\\_}{\perp x}\pm iB_{\\_}{\perp y}$.
For C${}_{\\_}{1h}$ symmetric divacancy configurations, the high symmetry axis
and the quantization axis of the silicon vacancy is 109.5 degree aligned to
the quantization axis of the central defect and the external magnetic field.

For point defect spin environments, the coupling Hamiltonian includes dipole-
dipole interaction term:

$H_{\\_}{\text{coup}}^{\text{point}}=\sum_{\\_}i-\frac{\mu_{\\_}0}{4\pi}\frac{g_{\\_}e^{2}\mu_{\\_}{B}^{2}}{\left|r_{\\_}{0i}\right|^{3}}\left(3\left(Sr_{\\_}{0i}\right)\left(\mathcal{S}_{\\_}ir_{\\_}{0i}\right)-S\mathcal{S}_{\\_}i\right)\text{,}$
(15)

where $r_{\\_}{0i}$ is the vector pointing from the position of the central
divacancy to the position of spin defect $i$.

We utilize our methodology to identify magnetic field values where spin bath
couplings and derived relaxation processes may be resonantly enhanced due to
level crossings and cross relaxation effects. In order to capture all
resonances, we simulate high concentrations spin environments, that is the
concentration is set to 1018 cm-3 corresponding to doped or highly irradiated
samples. In our quantitative study we focus on the depopulation of the highly
populated $m_{\\_}S=0$ state due to spin precession and spin flip-flops
induced by a spin bath. As the luminescence of divacancies is sensitive to the
population of the $m_{\\_}S=0$ state, our results are comparable with magnetic
field dependent PL plots. The intensity of the divacancy luminescence can be
estimated from our results according to the following expression

$\mathcal{L}\approx\mathcal{L}_{\\_}0\left(1-Cp_{\\_}{m_{\\_}S=\pm
1}\right)\text{,}$ (16)

where $p_{\\_}{m_{\\_}S=\pm 1}=\varrho_{\\_}{+1}+\varrho_{\\_}{-1}$ is the
population of the $m_{\\_}S=\pm 1$ states, $C$ is the optical spin contrast,
and $\mathcal{L}_{\\_}0$ is the luminescence of the bright $m_{\\_}S=0$ state.

Furthermore, we study spin relaxation quantitatively, i.e. calculate the spin
relaxation rate $T_{\\_}1^{-1}$, for the most relevant cases. To simulate spin
relaxation we do not include additional jump and decoherence operators that
would require adjustable rate parameters in the Lindblaian formalism.

## Data and code availability statement

The data that support the findings of this study and the software package used
to carry out the numerical simulations are available from the corresponding
author upon reasonable request.

## Acknowledgments

We acknowledge support from the Knut and Alice Wallenberg Foundation through
WBSQD2 project (Grant No. 2018.0071). Support from the Swedish Government
Strategic Research Area SeRC and the Swedish Government Strategic Research
Area in Materials Science on Functional Materials at Linköping University
(Faculty Grant SFO-Mat-LiU No. 2009 00971) is gratefully acknowledged. VI
acknowledges the support from the MTA Premium Postdoctoral Research Program,
the Hungarian NKFIH grants No. KKP129866 of the National Excellence Program of
Quantum-coherent materials project, and the NKFIH through the National Quantum
Technology Program (Grant No. 2017-1.2.1-NKP-2017-00001) and the Quantum
Information National Laboratory sponsored by Ministry of Innovation and
Technology of Hungary. N. T. S. acknowledges the support from the Swedish
Research Council (Grant No. VR 2016-04068), the EU H2020 project QuanTELCO
(Grant No. 862721). The calculations were performed on resources provided by
the Swedish National Infrastructure for Computing (SNIC) at the National
Supercomputer Centre (NSC) partially funded by the Swedish Research Council
through grant agreement No. 2018-05973.

## Author contributions

V.I. and O.B.L. wrote the manuscript with input from all authors. Calculations
were performed by O.B.L and V.I. with input from N.T.S. and I.A.A. The
computational results were analyzed with contributions from all authors. The
research was initiated and supervised by V.I.

## Competing interests

The authors declare no competing interests.

## References

  * [1] du Preez, L.  Ph.D. thesis, University of Witwatersrand (1965).
  * [2] Wrachtrup, J. & Jelezko, F.  Processing quantum information in diamond.  _Journal of Physics-Condensed Matter_ 18, S807–S824 (2006).
  * [3] Maze, J. R. _et al._ Properties of nitrogen-vacancy centers in diamond: the group theoretic approach.  _New Journal of Physics_ 13, 025025 (2011).
  * [4] Doherty, M. W. _et al._ The nitrogen-vacancy colour centre in diamond.  _Physics Reports_ 528, 1 – 45 (2013).
  * [5] Balasubramanian, G. _et al._ Nanoscale imaging magnetometry with diamond spins under ambient conditions.  _Nature_ 455, 648–651 (2008).
  * [6] Taylor, J. _et al._ High-sensitivity diamond magnetometer with nanoscale resolution.  _Nat. Phys._ 4, 810–816 (2008).
  * [7] Kucsko, G. _et al._ Nanometre-scale thermometry in a living cell.  _Nature_ 500, 54–58 (2013).
  * [8] Plakhotnik, T., Doherty, M. W., Cole, J. H., Chapman, R. & Manson, N. B.  All-Optical Thermometry and Thermal Properties of the Optically Detected Spin Resonances of the NV– Center in Nanodiamond.  _Nano Letters_ 14, 4989–4996 (2014).
  * [9] Schirhagl, R., Chang, K., Loretz, M. & Degen, C. L.  Nitrogen-vacancy centers in diamond: Nanoscale sensors for physics and biology.  _Annual Review of Physical Chemistry_ 65, 83–105 (2014).
  * [10] Degen, C. L., Reinhard, F. & Cappellaro, P.  Quantum sensing.  _Rev. Mod. Phys._ 89, 035002 (2017).
  * [11] Glenn, D. R. _et al._ High-resolution magnetic resonance spectroscopy using a solid-state spin sensor.  _Nature_ 555, 351–354 (2018).
  * [12] Schmitt, S. _et al._ Submillihertz magnetic spectroscopy performed with a nanoscale quantum sensor.  _Science_ 356, 832–837 (2017).
  * [13] Bernien, H. _et al._ Heralded entanglement between solid-state qubits separated by three metres.  _Nature_ 497, 86–90 (2013).
  * [14] Pfaff, W. _et al._ Unconditional quantum teleportation between distant solid-state quantum bits.  _Science_ 345, 532–535 (2014).
  * [15] Wehner, S., Elkouss, D. & Hanson, R.  Quantum internet: A vision for the road ahead.  _Science_ 362 (2018).
  * [16] Weber, J. R. _et al._ Quantum computing with defects.  _PNAS_ 107, 8513–8518 (2010).
  * [17] Kurizki, G. _et al._ Quantum technologies with hybrid systems.  _Proceedings of the National Academy of Sciences_ 112, 3866–3873 (2015).
  * [18] Zhang, J., Hegde, S. S. & Suter, D.  Efficient implementation of a quantum algorithm in a single nitrogen-vacancy center of diamond.  _Phys. Rev. Lett._ 125, 030501 (2020).
  * [19] Koehl, W. F., Buckley, B. B., Heremans, F. J., Calusine, G. & Awschalom, D. D.  Room temperature coherent control of defect spin qubits in silicon carbide.  _Nature_ 479, 84 (2011).
  * [20] Christle, D. J. _et al._ Isolated electron spins in silicon carbide with millisecond coherence times.  _Nat. Mater._ 14, 160–163 (2015).
  * [21] Son, N. T. _et al._ Divacancy in 4h-sic.  _Phys. Rev. Lett._ 96, 055501 (2006).
  * [22] Soltamov, V. A., Soltamova, A. A., Baranov, P. G. & Proskuryakov, I. I.  Room temperature coherent spin alignment of silicon vacancies in 4h- and 6h-sic.  _Phys. Rev. Lett._ 108, 226402 (2012).
  * [23] Widmann, M. _et al._ Coherent control of single spins in silicon carbide at room temperature.  _Nat. Mater._ 14, 164–168 (2015).
  * [24] Miao, K. C. _et al._ Universal coherence protection in a solid-state spin qubit.  _Science_ 369, 1493–1497 (2020).
  * [25] Christle, D. J. _et al._ Isolated spin qubits in sic with a high-fidelity infrared spin-to-photon interface.  _Phys. Rev. X_ 7, 021046 (2017).
  * [26] Falk, A. L. _et al._ Optical polarization of nuclear spins in silicon carbide.  _Phys. Rev. Lett._ 114, 247603 (2015).
  * [27] Bourassa, A. _et al._ Entanglement and control of single nuclear spins in isotopically engineered silicon carbide.  _Nature Materials_ 19, 1319–1325 (2020).
  * [28] Ivády, V. _et al._ Theoretical model of dynamic spin polarization of nuclei coupled to paramagnetic point defects in diamond and silicon carbide.  _Phys. Rev. B_ 92, 115206 (2015).
  * [29] Ivády, V. _et al._ High-fidelity bidirectional nuclear qubit initialization in sic.  _Phys. Rev. Lett._ 117, 220503 (2016).
  * [30] Li, Q. _et al._ Room temperature coherent manipulation of single-spin qubits in silicon carbide with high readout contrast.  _arXiv:2005.07876_ (2020).
  * [31] Whiteley, S. J. _et al._ Spin–phonon interactions in silicon carbide addressed by Gaussian acoustics.  _Nature Physics_ 15, 490–495 (2019).
  * [32] Anderson, C. P. _et al._ Electrical and optical control of single spins integrated in scalable semiconductor devices.  _Science_ 366, 1225–1230 (2019).
  * [33] Son, N. T. _et al._ Developing silicon carbide for quantum spintronics.  _Applied Physics Letters_ 116, 190501 (2020).
  * [34] Pham, L. M. _et al._ Enhanced solid-state multispin metrology using dynamical decoupling.  _Physical Review B_ 86, 045214 (2012).
  * [35] Bar-Gill, N., Pham, L. M., Jarmola, A., Budker, D. & Walsworth, R. L.  Solid-state electronic spin coherence time approaching one second.  _Nature Communications_ 4, 1743 (2013).
  * [36] Romach, Y. _et al._ Spectroscopy of surface-induced noise using shallow spins in diamond.  _Phys. Rev. Lett._ 114, 017601 (2015).
  * [37] Barry, J. F. _et al._ Sensitivity optimization for NV-diamond magnetometry.  _Reviews of Modern Physics_ 92, 015004 (2020).
  * [38] Takahashi, S., Hanson, R., van Tol, J., Sherwin, M. S. & Awschalom, D. D.  Quenching Spin Decoherence in Diamond through Spin Bath Polarization.  _Physical Review Letters_ 101, 047601 (2008).
  * [39] Jarmola, A., Acosta, V. M., Jensen, K., Chemerisov, S. & Budker, D.  Temperature- and Magnetic-Field-Dependent Longitudinal Spin Relaxation in Nitrogen-Vacancy Ensembles in Diamond.  _Physical Review Letters_ 108, 197601 (2012).
  * [40] Wang, H.-J. _et al._ Optically detected cross-relaxation spectroscopy of electron spins in diamond.  _Nature Communications_ 5, 4135 (2014).
  * [41] Wickenbrock, A. _et al._ Microwave-free magnetometry with nitrogen-vacancy centers in diamond.  _Applied Physics Letters_ 109, 053505 (2016).
  * [42] Wood, J. D. A. _et al._ Wide-band nanoscale magnetic resonance spectroscopy using quantum relaxation of a single spin in diamond.  _Phys. Rev. B_ 94, 155402 (2016).
  * [43] Wood, J. D. A. _et al._ Microwave-free nuclear magnetic resonance at molecular scales.  _Nature Communications_ 8, 15950 (2017).
  * [44] Ivády, V. _et al._ Photoluminescence at the ground-state level anticrossing of the nitrogen-vacancy center in diamond: A comprehensive study.  _Physical Review B_ 103, 035307 (2021).
  * [45] Koehl, W. F., Buckley, B. B., Heremans, F. J., Calusine, G. & Awschalom, D. D.  Room temperature coherent control of defect spin qubits in silicon carbide.  _Nature_ 479, 84–87 (2011).
  * [46] Falk, A. L. _et al._ Polytype control of spin qubits in silicon carbide.  _Nature Communications_ 4, 1819 (2013).
  * [47] Yan, F.-F. _et al._ Room-temperature coherent control of implanted defect spins in silicon carbide.  _npj Quantum Information_ 6, 1–6 (2020).
  * [48] Yang, W. & Liu, R.-B.  Quantum many-body theory of qubit decoherence in a finite-size spin bath.  _Physical Review B_ 78, 085315 (2008).
  * [49] Seo, H. _et al._ Quantum decoherence dynamics of divacancy spins in silicon carbide.  _Nature Communications_ 7, 12935 (2016).
  * [50] Onizhuk, M. _et al._ Probing the Coherence of Solid-State Qubits at Avoided Crossings.  _PRX Quantum_ 2, 010311 (2021).
  * [51] Shrivastava, K. N.  Theory of Spin–Lattice Relaxation.  _physica status solidi (b)_ 117, 437–458 (1983).
  * [52] Norambuena, A. _et al._ Spin-lattice relaxation of individual solid-state spins.  _Physical Review B_ 97, 094304 (2018).
  * [53] Radczyk, T., Hoffmann, S. & Goslar, J.  Applications of the Transport Integrals in Solid-State Physics and in Electron Spin Relaxation.  _Acta Physica Polonica A_ 5, 469–477 (2003).
  * [54] Hoffmann, S. K. & Lijewski, S.  Raman electron spin–lattice relaxation with the Debye-type and with real phonon spectra in crystals.  _Journal of Magnetic Resonance_ 227, 51–56 (2013).
  * [55] Redman, D. A., Brown, S., Sands, R. H. & Rand, S. C.  Spin dynamics and electronic states of N-V centers in diamond by EPR and four-wave-mixing spectroscopy.  _Physical Review Letters_ 67, 3420–3423 (1991).
  * [56] Simin, D. _et al._ Locking of electron spin coherence above 20 ms in natural silicon carbide.  _Physical Review B_ 95, 161201 (2017).
  * [57] Soltamov, V. A. _et al._ Relaxation processes and high-field coherent spin manipulation in color center ensembles in 6H-SiC.  _Physical Review B_ 103, 195201 (2021).
  * [58] Lin, W.-X. _et al._ Temperature dependence of divacancy spin coherence in implanted silicon carbide.  _arXiv:2104.12089 [quant-ph]_ (2021).
  * [59] Restrepo, O. D. & Windl, W.  Full First-Principles Theory of Spin Relaxation in Group-IV Materials.  _Physical Review Letters_ 109, 166604 (2012).
  * [60] Gugler, J. _et al._ Ab initio calculation of the spin lattice relaxation time T 1 for nitrogen-vacancy centers in diamond.  _Physical Review B_ 98, 214442 (2018).
  * [61] Astner, T. _et al._ Solid-state electron spin lifetime limited by phononic vacuum modes.  _Nature Materials_ 17, 313–317 (2018).
  * [62] Park, J., Zhou, J.-J. & Bernardi, M.  Spin-phonon relaxation times in centrosymmetric materials from first principles.  _Physical Review B_ 101, 045202 (2020).
  * [63] Xu, J. _et al._ Spin-phonon relaxation from a universal ab initio density-matrix approach.  _Nature Communications_ 11, 2780 (2020).
  * [64] Ivády, V.  Longitudinal spin relaxation model applied to point-defect qubit systems.  _Phys. Rev. B_ 101, 155203 (2020).
  * [65] Kawanishi, S. & Mizoguchi, T.  Effect of van der Waals interactions on the stability of SiC polytypes.  _Journal of Applied Physics_ 119, 175101 (2016).
  * [66] Davidsson, J. _et al._ First principles predictions of magneto-optical data for semiconductor point defect identification: the case of divacancy defects in 4H–SiC.  _New Journal of Physics_ 20, 023035 (2018).
  * [67] Magnusson, B. _et al._ Excitation properties of the divacancy in $4h$-sic.  _Phys. Rev. B_ 98, 195202 (2018).
  * [68] Ivády, V. _et al._ Stabilization of point-defect spin qubits by quantum wells.  _Nature Communications_ 10, 5607 (2019).
  * [69] Nagy, R. _et al._ High-fidelity spin and optical control of single silicon-vacancy centres in silicon carbide.  _Nature Communications_ 10, 1954 (2019).
  * [70] Kasper, C. _et al._ Influence of Irradiation on Defect Spin Coherence in Silicon Carbide.  _Physical Review Applied_ 13, 044054 (2020).
  * [71] Son, N. T., Carlsson, P., ul Hassan, J., Magnusson, B. & Janzén, E.  Defects and carrier compensation in semi-insulating $4h\text{$-$}\mathrm{Si}\mathrm{C}$ substrates.  _Phys. Rev. B_ 75, 155204 (2007).
  * [72] Greulich-Weber, S.  EPR and ENDOR Investigations of Shallow Impurities in SiC Polytypes.  _physica status solidi (a)_ 162, 95–151 (1997).
  * [73] Kuate Defo, R. _et al._ Energetics and kinetics of vacancy defects in $4h$-sic.  _Phys. Rev. B_ 98, 104103 (2018).
  * [74] Ivády, V. _et al._ Identification of si-vacancy related room-temperature qubits in $4h$ silicon carbide.  _Phys. Rev. B_ 96, 161114 (2017).
  * [75] Simin, D. _et al._ All-optical dc nanotesla magnetometry using silicon vacancy fine structure in isotopically purified silicon carbide.  _Phys. Rev. X_ 6, 031014 (2016).