Non-Markovian longitudinal spin relaxation model applied to point defect qubit systems

qubit systems Viktor Ivády1, 2, ∗ 1Wigner Research Centre for Physics PO Box 49, H-1525, Budapest, Hungary 2Department of Physics, Chemistry and Biology, Linköping University, SE-581 83 Linköping, Sweden (Dated: January 8, 2020) Abstract Controllable, partially isolated few level systems in semiconductors have recently gained multidisciplinary attention due to their widespread nanoscale sensing and quantum technology applications. Quantitative simulation of the dynamics and related applications of such systems is a challenging theoretical task that requires faithfully description not only the few level systems but also their local environments. Here, we develop a method that can describe relevant relaxation processes induced by a dilute bath of nuclear and electron spins. The method utilizes an extended Lindblad equation in the framework of cluster approximation of a central spin system. We demonstrate that the proposed method can accurately describe T1 time of an exemplary solid-state point defect qubit system, in particular NV center in diamond, at various magnetic fields and strain.


I. INTRODUCTION
Controllable solid-state spin systems have attracted considerable scientific and technological interest over the last decades. Point defect-based applications are among the most recent use of solid-state spins that allow full control over a set of electron and nuclear spins. NV center, substitution nitrogen-carbon vacancy complex point defect in negative charge state in diamond 1-5 is a magneto-optically active electron spin system that can be isolated to a large degree of freedom from the environmental disturbances. NV center's triplet electron spin can be initialized by pumping through optical excited triplet and meta stable singlet states. 3 The very same process gives rise to spin dependent optical decay that allows high fidelity read-out even at single NV center level. [6][7][8] In association with nuclear spins, NV center can implement few qubit nodes to realize high fidelity gates. 9 Coherence time may exceed a millisecond 10 and the qubit nodes can operate even above 600 K. 11 These attributes made NV center interesting for a broad range of quantum technology applications, especially in the field of quantum sensing [12][13][14][15] and quantum information processing 2,16,17 . Beside NV center, there have been several akin point defect qubit systems demonstrated in various wide band gap semiconductors. [18][19][20] Environmental spins, such as point defect and nuclear spins, play crucial role in spin relaxation and decoherence processes that are often the major limiting factors in quantum technology applications. Due to the complexity of some environmental spins' inner energy level structure, decay processes often depend on external control parameters, such as magnetic, electric, and microwave fields. In case of strong qubit-environment couplings, pumped point defect qubit systems serve as efficient spin polarization sources that can be utilized in hyperpolarization applications 21,22 either for enhancing the sensitivity of magnetic resonance experiments or for cooling environmental spins to reduce local magnetic field fluctuations.
Deeper understanding and numerical description of decoherence, spin relaxation, and polarization transfer over a wide range of environmental conditions are essential for advanced future applications. Lindblad master equation that describes Markovian decay processes is frequently applied when dynamical properties are considered. On the other hand, this approach relies on experimental decay rates and neglects the complexity of environmental interactions that may cause loss of quantitative accuracy and predictive power. To overcome these limitations numerous theoretical studies have been recently reported in this subject.
Temperature dependence of spin-phonon coupling induced spin relaxation of NV center was recently studied by analytic 4,40,41 and ab initio 42 approaches. Theoretical studies on spin bath induced spin relaxation processes focused on strong environmental coupling regions where dynamical nuclear polarization can be achieved. 21,[43][44][45][46] Much less attention have been paid, however, to the calculation of spin bath assisted relaxation processes and related decay time T 1 of point defect qubits at general control parameter settings where spin flip-flops are suppressed to a large degree of freedom. In a very recent study, CCE method was generalized to describe spin flip-flops of a NV center interacting with a bath of 13 C nuclear spins. 47 Rest of the article is organized as follows. Section II describes formulation of the theoretical approach and details of the implementation. Section III provides numerical results on the spin relaxation of NV center in diamond. Finally, section IV summarizes and concludes our findings.

II. METHODOLOGY
In this section we discuss cluster approximation of a many particle system in the framework of an extended Lindblad formalism to simulate spin relaxation processes. Hereinafter, we use the following terminology. We name subsystems of a closed or open system as spins. Spins are either elementary building blocks of the system or complex, many level systems. In the latter case spins can be defined based on the difference of internal and external coupling strength. We assume that inter-spin couplings are weaker than intra-spin couplings. Furthermore, we name processes that change the diagonal elements of spins' reduced density matrices as spin flip-flop processes.  Master equation of the open system S can be written aṡ where the Hamiltonian H 0 can be written as where h 0 is the Hamiltonian of the central spin, h i is the Hamiltonian of the coupled spin s i , and h 0i describes the coupling of the central spin and the spin bath spin s i . The last term on the right hand side of Eq. (1) accounts for environmental effects, such as temperature dependent effects and spin relaxation due to spins that are not included in S, through the Lindbladian E. The size of the problem, i.e the dimension of the Hilbert space, exponentially increases with n that makes exact solution unfeasible for large n.
To model the dynamics of S we divide it into a cluster C N of overlapping cluster systems, where N is the order of the cluster approximation. In first order cluster approximation cluster C 1 consists of n + 1 cluster systems c 0 and c i , where i ∈ {1, ..., n}. Expect from c 0 that includes only s 0 , all other cluster systems include the central spin and one coupled spin s i , see Fig. 1 illustration. Hamiltonians of the cluster systems can be written as, We may rationalize the above clustering by considering each c i cluster system as an implement to measure spin flip-flops induced by the coupled spin s i . Cluster system c 0 serves as a reference system where the central spin evolves freely without interacting with other spins.
Master equations of the cluster systems can be written aṡ where the dimensions of the density matrices are given by dim As the central spin is included in all cluster systems, there are altogether n + 1 definitions for the reduced density matrix of the central spin, i.e.
In the following subsections, we introduce couplings between the cluster systems to approximate the dynamics of the many spin system S. First, we extend Eq. (3) and Eq. (4) to account for an effective intra-spin bath field that time dependently shifts energy eigenvalues and preserves diagonal elements of the reduced density matrix of the central spin. Second, we extend Eqs. (5) by additional time dependent Lindbladian terms to account for interactions that induce spin flip-flops and cause variation of the diagonal elements of the central spin's density matrix.

Mean intra-spin bath field
The interaction Hamiltonian h 0i may include terms that do not induce spin flip-flops of the central spin but rather shift the energy levels. As such interactions alter the energy level structure of the system, they may affect the dynamics of the central spin too. According to Eq. (4) and Eq. (5), cluster system c i describes energy shifts solely due to spin bath spin s i , as the Hamiltonian h c i does not depend on other spin bath spin degrees of freedom. Energy shifts due to other spins, however, may be taken into account by introducing an effective field acting on the central spin.
This field is of course different in all cluster systems.
In order to account for the effective field of environmental spins included in other cluster systems, we extend h 0 as where β 0 and β i describe effective fields acting on s 0 in cluster system c 0 and c i , respectively.
To define β i , we first calculate the internal field α i in each cluster system c i obtained from the polarization of the environmental spin s i through a semi-classical formula and where I d i is the identity matrix of d i dimension. The extended Hamiltonians of the cluster systems can be written as, We note that the effective internal field defined by Eq. (10) and Eq. (11) act solely on the central spin. Note furthermore that the total effective field in each cluster systems is equal to β 0 as Tr i (β i + α i ) = β 0 . When nuclear spin bath is considered, the effective field β 0 of the polarized nuclear spin bath may be referred as the Overhauser field. Finally, note that the internal effective field can be utilized to account for dephasing effects in a semi classical approximation. Study of such processes is out of the scope of the present article. First of all, we extend master equations of cluster systems c 0 and c i by adding time dependent

Extended Lindbladian
Lindbladian terms L c 0 and L c i , aṡ and˙ where the Lindbladians are defined in the form of and whereḃ 0l andḃ il ≥ 0 are time dependent rates and C 0l and C il are Lindblad operators. We consider C 0l and C il operators that describe solely spin flip-flop transitions of the central spin. Therefore, C 0l and C il operators can be written as where C l Lindblad operators of d 0 dimension are identical for all cluster systems. Altogether d 0 (d 0 − 1) number of independent C l operators can be defined. We define these operators as where |m and |n are states of an orto-normal basis that spans the Hilbert space of s 0 . Hereinafter, we use l index as a shorthand notation of mn indices. Note that {C l } includes operators that drive spin flip-flops both forward and backward, i.e. C k and C T k ∈ {C l }. This condition is required by the irreversible effect of the extended Lindbladians in Eq. (16) and (17) and the positivity ofḃ 0l andḃ il rates. Furthermore, we note that Eq. (16) and (17) require that i.e. |m n| spin flip-flop processes are only possible when the population in the initial |n state is non-zero. To explicitly handle the exception when Tr Furthermore, we draw attention to a specific property of the definitions while where we explicitly use (mn) = l indices and dt is an infinitesimal time period. The right hand side of Eq. (21) and Eq. (22) describe how the diagonal elements of the density matrix c 0 change due to L c 0 over dt propagation.
We utilize L c 0 and L c i Lindbladians to carry out such spin flip-flops of the central spin in c 0 and c i that happen in cluster system c j due to coupling to s j for j = i. This effectively makes the diagonal elements of the reduced density matrix of the central spin to be identical in all cluster systems during the time evaluation, ie.
for any i at any time t, where diag is the vector of diagonal elements of . We utilize time dependent ratesḃ 0l andḃ il in Eqs. (14)- (15) to achieve this goal.
Before definingḃ 0l andḃ il , we need to quantify internal flip-flop rates in each cluster systems.
To do so, we defineȧ 0l andȧ il positive rates in such a way that the following equality are satisfied, and Note that the parentheses on the left hand side of Eq. (24) and Eq. (25)   To measure differences of the spin flip-flop rates between cluster systems c 0 and c i during the time evolution, we calculate The role of cluster system c 0 that includes only the central spin is apparent from Eq. (26). As s 0 in c 0 interacts with no environmental spin directly, a 0l measure flip-flop rates intrinsic to the central spin. Thus ∆ȧ il quantifies spin flip-flop rates of the central spin induced solely by environmental spin s i in cluster system c i .
Finally, let us define the time dependent ratesḃ 0l andḃ il entering Eq. (16) and Eq. (17)  is the additivity of the time dependent rates ∆a 0l and ∆a il . In Appendix A, we demonstrate that additivity is a good approximation for a non-entangled or partially entangled central spin system over an infinitesimal dt time evolution. Note that the first order cluster approximation, that neglects entanglement within the spin bath, and self-consistent solution of the equations ensure that the additivity holds at any time t during the time evolution of C 1 .
It is clear from the above discussion that the main approximation in the description of the central spin-spin bath coupling is the assumption of non-entangled spin bath. This maybe a good approximation when the coherence time of the spin bath specimens is shorter than the inverse coupling strength between the central spin and the spin bath specimens, i.e.
As we will see in the numerical calculations, this approximation is either satisfied or not satisfied depending on the type of environmental spins. Note, however, that in the latter case the approximation can be systematically improved by including spin bath interactions and inter spin bath correlations in higher order cluster approximations, see section II B.
It it apparent from the equations that the method does not require additional approximation of the Hamiltonian beyond the central spin approximation. Furthermore, restrictions due to central spin approximation can be remedied by using higher order cluster approximation, see section II B.
Note however that the approximation depends somewhat on the choice of the basis states used to span the Hilbert space of s 0 . In addition, the formalism ensures that conservable quantities are conserved in the cluster model. Let assume that the h 0 and h i are defined so that extensive quantity E is preserved in all cluster systems c 0 and c i . Property shown in Eqs. we assume that clustering is based on the indices of the environmental spins.
The Hamiltonian of higher order cluster system c I can be defined as where h ij describes interactions of environmental spins. Lindblad equation of the density matrix c I of cluster system c I can be written aṡ We note that definitions for the reference system c 0 are the same in every order of the cluster approximation. Furthermore, the definitions and approximations introduced in section II A 1 and II A 2 are irrespective to the order of the cluster approximation. In general the corresponding definitions can be obtained by substituting index i with index I.
Note that the approximations introduced in section II A 1 and II A 2 are systematically improvable by increasing the clustering order. Larger cluster systems describe coupling and entanglement that are completely neglected in first order cluster approximation. Ultimately, for N = n we return to the exact case.
In step (vi) of the propagation cycle, we quantify in each cluster system the spin flip-flops occurred during the short propagation calculated in the previous step. To do so, we restrict Eq. (24) and Eq. (25) to infinitesimal time evolution and obtain and where da 0l =ȧ 0l dt and da il =ȧ il dt.
Before discussing the next step of the cycle, we discuss through a few examples how to obtain and Eq. (35) as The first (second) summation on the right hand side adds up transition amplitudes of flip-flop processes that transform population to (from) state |j . From the solution of the above systems of linear equations one can obtain da 0l and da il .
It is important to notice that, da 0l and da il are not always uniquely defined. Altogether Having all da 0l and da il transition amplitudes defined, in step (vii) of the propagation cycle we and where In step (viii), we determine variation of the density matrices due to the extended Lindbladian defined in Eq. (14) and Eq. (15) as and In step (ix) of the propagation cycle, the cluster density matrices at t + dt are determined as and Finally, repetition of this procedure by substituting c 0 (t) and c i (t) by c 0 (t + dt) and c i (t + dt) allows one to simulate the dynamics of the many spin system S.

III. NUMERICAL SIMULATIONS
In this section, we computationally demonstrate through the example of NV center in diamond that the above described method can account for spin relaxation processes in different spin environments at various external fields. First, we provide spin Hamiltonians for the considered systems, then we describe the details of our ab initio density funtional theory calculations used to parameterize NV center-nuclear spin bath interactions. In the subsequent sections we study different spin bath induced relaxation processes of an NV center's spin polarization.
A. Background and methodology

Spin Hamiltonian
We study NV center-spin bath coupled systems. In particular, we consider P1 center (neutral substitutional nitrogen atom with spin-1/2 ground state), NV center, and 13 C nuclear spin reservoirs interacting with the central NV center (s 0 ). For simplicity, we ignore NV centers' nitrogen nuclear spin that gives rise only to a fine structure at the ground state level anti crossing (GSLAC) 43 where S z is the spin z operator defined in a coordinate system with z axis parallel to the NV axis, D = 2.870 GHz 6 is the zero field splitting, g e is the electron g-factor, µ B is the Bohr magneton, and d and d ⊥ account for parallel and perpendicular strain coupling 48 , respectively, where variables with tilde symbol are defined in a coordinate system with z axis parallel to the C 3v axis of the Jahn-Teller distorted configuration of P1 center, J i is the 14 N nuclear spin operator, P = 5.01 is the quadrupole splitting for which we use the value of the NV center 49 , A is the hyperfine tensor whose diagonal elements, A zz = 114 MHz and A xx = A yy = 81 MHz, are determined by our first principles electronic structure calculations, see below, g 14N = 0.4038 is the nuclear g-factor of the spin-1 14 N nucleus, µ N is the nuclear magneton, where S i is the spin operator defined in a coordinate system with z axis parallel to the symmetry axis of the environmental NV center, and where I i and g 13C = 1.4048 are the nuclear spin operator and the nuclear g-factor of the spin-1/2 13 C nucleus, respectively.
Coupling tensors between the central NV center' spin and electron spin bath specimens, such as P1 centers and environmental NV centers, are obtained by neglecting spacial distribution of the spin densities through the dipole-dipole interaction Hamiltonian, where S 0 and S i are spin operators of the central spin and the coupled spin, respectively, µ 0 is the vacuum permeability, and r 0i is a vector pointing from the central spin to the coupled spin. When nuclear spin bath is consider, NV center-nuclear spin couplings are described by the hyperfine interaction Hamiltonian, where I i and A i are the nuclear spin operator and the hyperfine coupling tensor in cluster system j, respectively.
The following cluster Hamiltonians are used to model different spin environments,

First principles electronic structure calculations
Hyperfine coupling tensors are key quantities when 13 C nuclear spin bath is considered. We use first principles Density Functional Theory (DFT) electronic structure and subsequent hyperfine tensor calculations to obtain relevant coupling tensors of NV center and P1 center spin systems in diamond. In our DFT calculations we use 1728 atom supercell, HSE06 hybrid functional 50 , PAW core potentials 51 , and plane wave basis set of 420 eV as implemented in VASP 52 .
It is possible to calculate hyperfine interaction with high accuracy 53 for atomic sites in close vicinity of the NV center, however, for farther sites the hyperfine interaction suffers from considerable finite size effects in supercell methods 54,55 . To overcome this issue we utilize real space grid combined with PAW method to calculate hyperfine tensors. Fermi contact term, dipole-dipole interaction within the PAW sphere, and core polarization correction are calculated within the PAW formalism 53 from the convergent spin density. Dipolar hyperfine contribution from spin density localized outside the PAW sphere is calculated by using a uniform real space grid. This procedure allows us to obtain hyperfine coupling tensors including no effects from periodic replicas of the spin density due to the periodic boundary condition. Additionally, we can calculate hyperfine tensors for atomic sites outside the boundaries of the supercell by neglecting Fermi contact interaction in that region.

B. T 1 time of NV center
In the following computational example, we study NV center's longitudinal spin relaxation in different spin environments over a wide range of external magnetic fields and strain. In the simulations, we neglect spin-orbit and phonon assisted decay processes. The former effect is negligible in the ground state of the NV center, while the latter approximation is valid at low temperatures (below ≈50 K).
First, we investigate spin relaxation due to P1 center spin environment. In the simulations, we considered an ensemble of 50 randomly generated configuration of 31 P1 centers. The concentration of the P1 center spin bath is set to 50 ppm, which corresponds sample S2 in Ref. [56]. Except Spin relaxation rate (1/T 1 ) as a function of the extarnal magnetic field is depicted in Fig. 3(b).
Note that the distribution of the relaxation rates over the considered ensemble is highly asymmetric, meaning that there is a low but non-zero probability of finding centers with very large relaxation rates. Such a distribution cannot be faithfully characterized by the usual statistical quantities, such as mean and standard deviation, therefore in Fig. 3(b), we provide additional quantities to properly describe the relaxation rate distribution. We depict relaxation rate of configurations with the lowest (Min) and the highest (Max) relaxation rates, as well as, the mean (Mean) and the mode (Mode) of the distribution. Note that the mean and mode relaxation rates are different due to the asymmetric distribution, implying that the average T 1 time can, in general, be shorter than the most probably T 1 time of individual centers in an ensemble.
The simulations reveal two magnetic field values, 51 mT and 102 mT, where enhanced spin relaxation takes place. By looking at the energy level structure of NV-P1 coupled system depicted in Fig. 3(a), the relaxation peaks can be assigned to level crossings. Reduction of energy gaps enhances spin flip-flop rates, which is captured by the simulations. We note that at 0 mT, no enhanced relaxation can be observed despite the crossing of the levels at this field. This is due to the fact that the spin-1 NV center exhibits a large zero-field splitting, while electron spin sublevels of the P1 center are degenerate at zero magnetic field. Therefore, couplings are efficiently suppressed.
The relaxation peak at 51 mT exhibits a fine structure not fully resolvable in Fig. 3(b). In Fig. 3(c), we depict the relaxation rate of a representative spin bath configuration, including either parallel or 109 • aligned P1 centers. In both cases a five-peak fine structure can be seen with different spacings due to the different orientation of the hyperfine principal axis in the two cases.
A different fine structure is obtained at 102 mT, see Fig. 3(d). The peaks at 51 mT are due to center. Related PL signatures were recently reported in Ref. [60].
The main approximation of the methodology proposed in this article is the neglect of entanglement between environmental spins. For P1 center spin bath we obtained T 1 > 1 ms for most of the magnetic field values considered in the simulations. As the T P1 2 time of the P1 centers at 50 ppm is expectedly much shorter than 1 ms, the relation T P1 2 << T NV 1 is satisfied. This validates the approximation of non-entangled spin bath. Next, we investigate the magnetic field and strain dependence of the spin relaxation rate of a central NV center interacting with a number of environmental NV centers. Settings for the simulations are the same as for P1 center environment, except the defect concentration, which is set to 12 ppm, in accordance with sample S2 in Ref. [56], and the initialization of the environmental NV centers, where we set 90% polarization in the m S = 0 state.
The energy level structure and the corresponding theoretical spin relaxation rate of a central NV center in 109 • aligned NV center environment are depicted in Fig. 4(a) and (b), respectively. We obtain highly anisotropic distribution for the relaxation rates characterized by the minimal, maximal, mean, and mode values in Fig. 4(b). Three relaxation peaks can be found in the investigated magnetic field interval at 0 mT, 59 mT, and 102 mT. Related PL features at 59 mT were reported in experiment. 58,60 The peaks correspond to crossings between the energy levels of the coupled two NV center system depicted in Fig. 4(a). Since the central NV center and the environmental NV centers exhibit the same zero-field splitting, spin states can be mixed at zero magnetic field that gives rise to a relaxation peak, in contrast to the P1 center environment.
Magnetic field oriented NV center environment gives rise to a distinct relaxation pattern, see The degenerate states can be mixed by the dipole-dipole coupling that gives rise to a constant very high relaxation rate. We note that this high relaxation rate can be substantially reduced in experiment due to two effects. (i) The relaxation rate is linearly depends on the polarization difference between the central NV center and the environmental NV centers. In our simulations we set 10% difference which may be higher than in sample upon measurement. (ii) In the simulation the states are degenerate due to the identical level structures of the two centers, however, magnetic field and strain inhomogeneities may make the centers distinguishable.
In order to investigate these effects in NV center-NV bath systems, we actuate parallel and perpendicular strain on the central spin and environmental NV centers of parallel and 109 • alignment.
In Fig. 5(a), polarization dependence of the 109 • aligned NV center bath induced spin relaxation rate is depicted. We find that the relaxation rate is approximately a factor of three larger in the case of the polarized, 90% in m S = 0 state with 109 • aligned quantization axis, NV center bath.
Due to the optical pumping, polarization of the NV bath is expected, however, at magnetic field of enhanced coupling to other spin species, e.g. at B = 59 mT, low polarization is more probable.
Strain dependence of the relaxation rate at specific magnetic fields is depicted in Fig. 5 that the coupling of the NV center to perpendicular strain is an order of magnitude larger than the coupling to parallel strain, thus the range of considered strain field is larger in the former case.
Note, furthermore that, similar but reduced effects can be found at B = 59 mT, not shown. At B = 102 mT we see, however, distinct behavior. Relaxation rate insensitive to a large extent to parallel strain applied on the central NV center and to both parallel and perpendicular strain applied on environmental NV center. On the other hand perpendicular strain applied on the central NV center mixed the spin states efficiently 48 that gives rise to a prominent increase of the relaxation rate. Relaxation rate distribution of NV centers in parallel aligned NV center environment is characterized in Fig. 5(d). It is apparent from the figure that the strain shift reduces the relaxation rate substantially. This effect, however, vary considerably with the spin bath configurations.
When the central spin-environment couplings are weak, even a small strain shifts can induce large reductions in the rates.
Similarly to the P1 center environment, the condition T NV 2 << T NV 1 is expectedly satisfied in the modeled sample. Therefore, the approximations of the applied method hold. Next, we investigate NV center-13 C spin bath systems. The settings for the simulations are similar as for the P1 center environment, except for the concentration of the spin defects, for which we used the natural abundance of 13 C. Due to the fairly simple level structure of the NV center-13 C nuclear spin system, see Fig. 6(a), the relaxation rate curves shown in Fig. 6(b) exhibit only a single peak at 102 mT that correspond to the GSLAC 43 .
For simplicity, here we use the same, first order cluster approximation as before, i.e. cluster systems include the central spin and only one nuclear spin. As the nuclear spins have very long coherence time, the relation T 13C is solely due to 13 C spins, may not be satisfied. In this case an overestimation of the relaxation rates is expected. Therefore, the results presented in Fig. 6(b) may be considered as an upper bound for the 13 C spin bath induced relaxation.
Finally, we combine our theoretical results in order to compare with experimental measure- (c) Combined theoretical relaxation rate is compared with the experimental spin relaxation rate reported in Ref. [56] for sample S2 at low temperature.
ments reported for sample S2 in Ref. [56]. The total spin relaxation rate can be given as The theoretical relaxation rate curves with uncertainties deduced from experimental uncertainties in the defect concentrations are depicted in Fig. 7(a) and (b). To determine the uncertainties, we use linear concentration dependence for the relaxation rates 56 . As there is no available data on the strain and magnetic field inhomogeneity as well as the polarization variation of parallel and 109 • aligned NV centers in the sample modeled here, we make the following assumptions. We assume 1) O(0.1 MHz) parallel strain and magnetic field inhomogeneity, 2) O(1 MHz) perpendicular strain, and 3) 1% variance in the polarization of parallel NV centers. The resultant curves are plotted in Fig. 7(b). It is apparent from the results that environmental NV centers have dominant effect on the spin relaxation rate.
When compared with experiment, we find that the theoretical curve follows within error bars the measurements over a wide range of the magnetic field considered in the experiment. Higher discrepancy can be seen at B = 59 mT, where the theoretical curve overestimates the experimental relaxation rate. This can be attributed to the neglect of depolarization of the environmental NV centers. As we have seen in Fig. 5(a), depolarization of the bath reduces relaxation rate.
Depolarization of parallel and 109 • aligned NV centers is expectedly mutual when they couple at B = 59 mT. Inclusion of this effect can lower the theoretical relaxation rate to the level of experimental measurements.
The numerical results demonstrate that the proposed theoretical method can account for the reported magnetic field dependent spin relaxation patterns induced by P1 centers, NV centers, and 13 C nuclear spins. This is due to the non-approximate description of the pair interactions between the central spin and the environmental spins. Furthermore, numerical simulations validate the approximations introduced by first order cluster approximation in the case of P1 center and NV center spin environments. This makes it possible to obtain quantitative results comparable with experiment.

IV. SUMMARY
In summary this paper describes a non-Markovian spin bath model for calculating spin relaxation effects in central spin approximation. To this end an extended Lindbladian formalism was introduced to account for spin flip-flops in a many spin system. Validity of the approximation is In the numerical simulations NV center's spin relaxation rate (1/T 1 ) was investigated. P1 center, NV center, and 13 C spin baths are considered at various magnetic fields and strain. The method captures all the known characteristics of the relaxation rate of specific spin bath systems. By taking all the relevant relaxation effects into account, the theoretical spin relaxation rate curve is quantitatively comparable with the measured one over a wide range of magnetic fields. for which one gets