# Fileset [MDR_JJAP_Inoue.pdf](https://mdr.nims.go.jp/filesets/251c4ab5-d66a-474d-8a2c-ce27df8e42d3/download) ## Creator [Jun-ichi Inoue](https://orcid.org/0000-0001-6743-4258) ## Rights © 2024 The Japan Society of Applied Physics
This is an author-created, un-copyedited version of an article accepted for publication/published in Japanese Journal of Applied Physics. IOP Publishing Ltd is not responsible for any errors or omissions in this version of the manuscript or any version derived from it. The Version of Record is available online at https://doi.org/10.35848/1347-4065/ad1719.[In Copyright](http://rightsstatements.org/vocab/InC/1.0/) ## Other metadata [Dynamic frequency shift in NV⁻ center in diamond induced by anisotropic in-plane vacuum fluctuation](https://mdr.nims.go.jp/datasets/0c85fb10-eda5-451f-8d0b-a5b025c2c88f) ## Fulltext Dynamic Frequency Shift in NV− Center inDiamondJun-ichi Inoue ∗National institute for Materials Science, Tsukuba, Ibaraki 305-0044, JapanApril 16, 2024AbstractSensor applications of negatively charged nitrogen vacancy (NV− center indiamond are now in practical use, yet for finer sensitivity, a comprehensive under-standing of various kinds of sources that cause detrimental relaxation and dampingis still required. During the course of theoretical study regarding this, we foundthat Gaussian white noise with the zero-mean has a substantial effect, which man-ifests itself in a period of free induction decay (FID) oscillation. This effect isexperimentally detectable through comparison with zero-field splitting fixed by,e.g., optically detected magnetic resonance (ODMR). The result is corroboratedby a different analytical framework, Lindblad master equation. Our finding in FIDoscillation period, or an equivalent energy shift, is concluded to fall into a class ofdynamic frequency shift.1 IntroductionSince usefulness of color centers in diamond was recognized, extensive studies re-garding these have been continued. Among others, a nitrogen vacancy center withan excess electron, NV−, is the most promising one. It is now used as sensors andother applications, to name a few, for local magnetic fields [1, 2, 3, 4], electron [5]and viscous [6] flows in graphene, temperature [7, 8, 9], and even vortexes in cupratehigh-temperature superconductors [10]. One of the reasons behind this ability lies inits long relaxation times even at room temperature. Aiming at finer sensitivity, furthereffort mainly in terms of material science has been made to obtain clean and high-puritysamples, which are expected to be ideal for longer relaxation times. However, as shownin the counter-intuitive example that ultra long relaxation time was reported in a sys-tem with phosphorus intentionally-doped [11], underlying physics regarding relaxationphenomena is not crystal clear, and a comprehensive understanding of various sourcesthat cause relaxation and damping is highly required. As such, NV−1 is still an idealand unique test bed to study, in particular, non-equilibrium open quantum physics, and∗E-mail: inoue.junichi@nims.go.jp1along this line, the inhomogeneous relaxation time, T∗2 , of NV− has been studied asa typical example of relaxation times. While primary interests lie in what sources asnoise dominate relaxation phenomena and how, in the course of the study, we theoreti-cally find that anisotropic noise with the zero-mean can bring about an energy shift innon-trivial manner, which is the main result of this work.The relaxation time T∗2 characterizes free induction decay (FID) and Ramsay fringe.In spite of many advanced theoretical and computational methods to discuss the relax-ation time [12, 13], its core can be captured by simply considering an effective modelof a two-level system, or a spin 1/2, in a fictitious magnetic field along z-directionwhose amplitude fluctuates, and in this case, pure dephasing part of T∗2 is focused. Themodeling is verified by the following consideration [14].An NV− center, where direction toward N- from V-sites is taken as z axis, locallyhas C3v point group symmetry. An electronic ground state of NV− is known to be S = 1triplet states, whose behavior is described by an S = 1 spin Hamiltonian respecting thesymmetry. The most general Hamiltonian for the triplet state up to the second order ofspin operators can be read as [14] (Zeeman term against an externally applied magneticfield is omitted in order to avoid unnecessary confusion in the later discussion)H = (d0 + ξ)S2z + −Πx(SxSy + SySx) − Πy(S2x − S2y), (1)where ®S = (Sx, Sy, Sz) are S = 1 operators, and Sz is diagonal with its eigenvaluesm = {+1, 0,−1}. In the Hamiltonian, d0 ≈ 2.87GHz is zero-field splitting, and the otherterms with constants, ξ and Πx,y , represent contribution from local lattice deformationand electric field induced by intrinsic charge density fluctuation. Diagonalize theHamiltonian and pick up the lower two energy states, then it is equivalent to a two-levelsystem in a fictitious static magnetic field.As a Hamiltonian describes unitary time evolution of a given system, some “noise” isnecessary to describe pure dephasing. For the purpose, the constants in the Hamiltonian,ξ and Πx,y , are allowed to fluctuate around their certain mean values, by regarding theparameters in Eq, (1) as stochastic variables [14]. Since NV− centers in most caseslie at room temperature, thermal lattice fluctuation and accompanied charge fluctuationare naturally considered as its origin. However introduction of such a fluctuation inthe Hamiltonian turns itself into time dependent one. In this course, the eigenvalueproblem delivered by the static Hamiltonian might lose its firm ground. The ideaof the authors of Ref. [14] is that, by considering a snapshot of the time dependentHamiltonian, the stochastic effect is incorporated as fluctuation of energy differencebetween the two states in the effective two-level system. The system obtained is thusequivalent to a spin 1/2 in fictitious fluctuating magnetic field, Bz(t). The analysis ofT∗2 on the basis of this model unveiled qualitatively different behavior of T∗2 against anexternally applied magnetic field, depending on which parameter’s fluctuation in theHamiltonian dominates, and showed agreement with experiments [14].The essentials of Ref. [14] are to assume that the influence of stochastic variables inthe Hamiltonian can be read as fluctuation of the energy levels derived by diagonalizingthe snapshot Hamiltonian. Here, being fueled by its success and reviewing the work,we notice that a margin still remains in this procedure: not only diagonal componentsin the effective two-level system, or energy levels, but also off-diagonal components2in diagonalized snapshot Hamiltonian can also be considered as stochastic variables,which have the zero-means. This corresponds to encompassing vacuum fluctuationof fictitious in-plane magnetic fields in the effective two-level system. To examineconsequences of this extension, below we consider a spin 1/2 system in a fictitiousmagnetic field whose direction fluctuates as well as amplitude, and show that the in-plane vacuum fluctuation brings about an energy shift, or an equivalent frequency shiftas an experimentally detectable effect. This frequency shift falls into a class of dynamicfrequency shift defined below.2 Theory2.1 stochastic modelConsider dynamics of a spin 1/2 exposed to a fictitious fluctuating magnetic field,®B(t) = {Bµ(t)} (µ = x, y, z), which is assumed to involve Gaussian white noise with themean values (Bx(t), By(t), Bz(t)) = (0, 0, d0) and fluctuation amplitude γµ. A connectionof γz to the parameters in the original Hamiltonian, Eq. (1), can be postulated asin Ref.[14], while ones of γx and γy can not, even though their origins could betraced back to the parameters in the original Hamiltonian. This is because these twoparameters describe fluctuation around the zero off-diagonal components of the matrixalready diagonalized, where dependence of γx,y on the parameters in the Hamiltonianis unseen. Therefore γx and γy are phenomenologically treated, and the values of γµ’sare reasonably assumed to be different from each other, partly because of an intrinsicsymmetric property. Here we would like to emphasize, in advance, that it is unequalityof γx and γy which plays an essential role in this work.A two-component state vector, |ψ(t)⟩, obeys (ℏ = 1)i∂∂t|ψ(t)⟩ = 12®B(t) · ®σ |ψ(t)⟩, (2)with usual Pauli matrix, ®σ. Introducing Bloch vector, defined by ®r = (x, y, z) ≡⟨ψ(t)| ®σ/2|ψ(t)⟩, the dynamics is equivalently described by Bloch equationddt®r = ®B(t) × ®r = ©«0 −Bz(t) By(t)Bz(t) 0 −Bx(t)−By(t) Bx(t) 0ª®¬ ©«xyzª®¬ , (3)where the matrix in Eq. (3) is simply denoted by M(t).Time profiles of FID of a two-level system are described by x(t) and y(t) in thisusage. With a time-ordering operator T , the formal solution of this equation can bewritten as®r(t) = T exp[∫ t0M(t ′)dt ′]®r0, (4)for a certain initial state ®r0 ≡ (x0, y0, z0) at t = 0. To deal with the time-orderingoperator, we follow the method given in Ref. [15]. Divide the full-time evolution from3t ′ = 0 to t ′ = t into a series of time slices with width ∆, which is finally taken by thelimit of ∆→ 0. Since the noises are assumed to be white, those belonging to differenttime slices are uncorrelated, and the matrix is regarded as “constant” during the thintime slice ∆. Taking ensemble average, denoted by overbar, the solution can then bewritten as®r(t) =[exp(M∆)] t/∆®r0. (5)A series expansion with respect to ∆ gives exp[M∆] = Mxyz + O(∆3), whereMxyz =©«1 −(γy2 +γz2)∆ −d0∆ 0d0∆ 1 −(γz2 +γx2)∆ 00 0 1 −(γx2 +γy2)∆ª®®®¬ . (6)In the derivation, the mean values, Bx(t) = By(t) = 0 and Bz(t) = d0, and the secondorder correlationsBµ(t)Bν,µ(t) = 0, B2µ(t) =γµ∆, (7)are substituted. The last condition is consistent with white noise in time slice withwidth ∆. Taking the limit ∆→ 0, the solution turns to®r(t) =[Mxyz] t/∆ ®r0∆→0−→ exp[−Nxyz t]®r0, (8)whereNxyz =©«12(γy + γz)d0 0−d012 (γx + γz) 00 0 12(γx + γy)ª®¬ . (9)The final form of z(t) is immediately obtained asz(t) = exp(−γx + γy2t)z0, (10)where z0 = 0 is usually prepared after π/2 pulse injection. As for x(t) and y(t), let Ωbe a real constant asΩ =√16d20 − (γx − γy)2, (11)whose reality comes from a reasonable assumption that, for NV− center, d0 > γx,y,z .Then the final solutions are represented as(x(t)y(t))= exp[−Nxyt] (x0y0), (12)4where the matrix exponential isexp[−Nxyt]= exp(−γall4t) [N (1)xy cos(Ω4t)+ N (2)xy sin(Ω4t)], (13)withN (1)xy =(1 00 1), (14)N (2)xy =1Ω(γx − γy −4d04d0 −γx + γy), (15)γall = γx + γy + 2γz . (16)Note that both matrices, N (1,2)xy , are regular as their determinants are identically unity.From these equations, FID signal shows sinusoidal oscillation with exponential decay,whose relaxation time characterized by γall corresponds to T∗2 . The in-plane vacuumfluctuation γx,y , as expected, additively contributes to the relaxation constant, whoseamount is half of that of longitudinal fluctuation γz .In the above equations, our main result is included: a period of the oscillation, τ, isfound to be also influenced by the in-plane vacuum fluctuation, in the form |γx − γy |, asτ =2πΩ/4 = 42π√16d20 − (γx − γy)2. (17)Since the off-diagonal components in Eq. (2) are stochastic variables with the zero-means, substantial effects in energy are, at first sight, not expected. However, Eq. (17)indicates that the in-plane fluctuation causes an energy, or frequency shift. The pointof interest is that a condition of non-zero γx and γy is not sufficient, and unequality ofthe two is required to obtain the finite contribution.We claim that the effect caused by the anisotropic in-plane vacuum fluctuation isexperimentally detectable. When an oscillation period of FID shows a deviation fromthe one obtained by zero-field splitting, the origin of it can be ascribed to the vacuumfluctuation effect discussed here. In previous studies, the oscillation period in FID isconsidered to be solely determined by zero- field splitting d0, whose value is ∼ 2.87GHz, measured by e.g. ODMR at zero external magnetic field. Since ODMR isfree from the frequency shift induced by the in-plane vacuum fluctuation, careful FIDmeasurements can unveil the effects as a difference of the period from the conventionalvalue 2π/d0.The system addressed above is omnipresent [16, 17], and indeed often referred to asAbragam model in literature [18]. However most of the previous studies conventionallyassumed γx = γy . Once this equality is imposed, the effect due to the in-plane vacuumfluctuation in frequency shift disappears, as in Eq. (17). In this respect, our finding isthought to stem from a property given by the condition of γx , γy , detail of which isdescribed in latter section.One might think that the substantial effect from the in-plane vacuum fluctuation asin the form of frequency shift is an artifact of naive mathematical handling of equations.5Therefore in order to corroborate the results, we demonstrate below that the identicalsolutions are obtained within the framework of master equation with Lindblad form,which is one of the most-established methods for open quantum physics [12, 13].2.2 Equivalent Lindblad equationConsider a phenomenological Lindblad master equation for a density operator ρ of ageneral two-level systemdρdt= −i[H, ρ] +∑µ=x,y,z(LµρL†µ − 12L†µLµρ −12ρL†µLµ), (18)where Hamiltonian, describing unitary time evolution, is H = (d0/2)σz , while Lindbladoperators, characterizing relaxation processes, are given by {Lµ} = (√γµ/2)σµ withphenomenological relaxation constants γµ. The equivalence of the treatment in thelast subsection and that with Eq.(18) is demonstrated by agreement of their explicitsolutions.From the Lindblad equation Eq.(18), analytical solutions of Sµ=x,y,z(t) ≡ Tr[ρ(t)σµ]for an initial condition (Sx(t = 0), Sy(t = 0), Sz(t = 0)) = (sin θ cos ϕ, sin θ sin ϕ, cos θ) ≡(x0, y0, z0) on the Bloch sphere are straightforwardly obtained asSx(t) = sin θ exp[−γx + γy + 2γz4t]×cos ϕ cosh©«√(γx − γy)2 − 16d204tª®®¬+(γx − γy) cos ϕ − 4d0 sin ϕ√(γx − γy)2 − 16d20sinh©«√(γx − γy)2 − 16d204tª®®¬ , (19)Sy(t) = sin θ exp[−γx + γy + 2γz4t]×sin ϕ cosh©«√(γx − γy)2 − 16d204tª®®¬−(γx − γy) sin ϕ − 4d0 cos ϕ√(γx − γy)2 − 16d20sinh©«√(γx − γy)2 − 16d204tª®®¬ , (20)Sz(t) = cos θ exp[−γx + γy2t]. (21)In particular, when |γx−γy | < 4d0, using a real constant as before,Ω =√16d20 − (γx − γy)2,6above results are summarized as(Sx(t)Sy(t))= exp[−γall4t]× ©«cos(Ωt4)+γx−γyΩsin(Ωt4)− 4d0Ωsin(Ωt4)4d0Ωsin(Ωt4)cos(Ωt4)− γx−γyΩsin(Ωt4)ª®¬(x0y0), (22)andSz(t) = exp[−γx + γy2t]z0. (23)These are identical with Eqs.(12) and (10), respectively. By this agreement, we can saythat the frequency shift caused by the anisotropic in-plane vacuum fluctuation γx andγy is now on firm ground.3 Discussion3.1 role of anisotropy in in-plane vacuum fluctuationThe peculiarity of γx , γy can be ascribed to an algebraic property of the Hamiltonianand the Lindblad operators in Eq.(18). For the two operators, the following commutationrelation [[H, ρ],∑µ=x,y,zγµ(LµρL†µ − 12L†µLµρ −12ρL†µLµ)]∝ d0(γx − γy), (24)holds [19], meaning that any non-zero component of the matrix commutator is propor-tional to d0(γx − γy). This commutation relation indicates that only when γx = γy ,the evolution induced by the Lindblad operators is invariant under a rotation aroundthe z axis, latter of which is induced by the Hamiltonian. We can expect that this non-commutativity is the origin of the frequency shift appeared in FID oscillation period,and if it is the case, the shift should be referred to as dynamic frequency shift, by thedefinition of the shift [20, 21]. This expectation is indeed supported by the followinganalysis.The abstract algebraic property in Eq. (24) turns more explicit once the masterequation is given in a matrix form. From Eq. (18), a set of differential equations foreach component of the density matrix, ρi j ≡ ⟨ϕi |ρ|ϕ j⟩ ({i, j} = {1, 2}), is written asddt©«ρ11(t)ρ12(t)ρ21(t)ρ22(t)ª®®®¬ =(MHxyz + MLxyz) ©«ρ11(t)ρ12(t)ρ21(t)ρ22(t)ª®®®¬ , (25)7whereMHxyz =©«0 0 0 00 −id0 0 00 0 +id0 00 0 0 0ª®®®¬ , (26)MLxyz =©«−γx+γy4 0 0 γx+γy40 −γx+γy+2γz4γx−γy4 00 γx−γy4 −γx+γy+2γz4 0γx+γy4 0 0 −γx+γy4ª®®®®¬. (27)In these equations, MHxyz characterizes unitary evolution, while MLxyz relaxation. Thesematrices do not commute when d0(γx − γy) , 0, as[MHxyz, MLxyz]= d0(γx − γy)©«0 0 0 00 0 −i/2, 00 i/2, 0, 00 0 0 0ª®®®¬ . (28)This is a representation of Eq.(24) in case of the present Hamiltonian and Lindbladoperators.The origin of this non-commutativity comes from the (2, 3) and (3, 2) componentsof MLxyz , both of which are proportional to γx − γy . This statement turns more obviousby comparing with a different master equation. Take the most typical textbook exam-ple [22], a Lindblad master equation with H = (d0/2)σz and the Lindblad operator{Lµ=±,z} = (√γ+/2)σ+ , (√γ−/2)σ− , (√γz/2)σz}, where σ± ≡ (1/2)(σx ± σy). Thisequation describes photon emission and absorption with phase relaxation. In the samemanner as Eq. (25), two matrices in a set of differential equations in the present case,MH±z + ML±z , are MH±z = MHxyz andML±z =©«−γ−4 0 0 γ+40 −γ++γ−+4γz8 0 00 0 −γ++γ−+4γz8 0γ−4 0 0 −γ+4ª®®®¬ , (29)which now commute: [MH±z, ML±z] = 0. From this comutativity, an FID oscil-lation period derived from the master equation with H and {L±,z} is in turn an-ticipated to be free from the dynamic frequency shift. To see this, solutions of(Sx(t), Sy(t), Sz(t)) for the master equation with initial condition (Sx(0), Sy(0), Sz(0)) =(sin θ cos ϕ, sin θ sin ϕ, cos θ), should be given. They are readily obtained asSx(t) = exp[−γ+ + γ− + 4γz8t]cos(d0t + ϕ) sin θ, (30)Sy(t) = exp[−γ+ + γ− + 4γz8t]sin(d0t + ϕ) sin θ, (31)Sz(t) = exp[−γ+ + γ−4t] (cos θ − γ+ − γ−γ+ + γ−)+γ+ − γ−(γ+ + γ−), (32)8and indeed, the period of the sinusoidal behavior is solely fixed by d0. Althoughthe two Lindblad equations with {Lx,y,z} and with {L±,z} seem to be equivalent fromthe view points of Lie algebra [23], the relaxation processes described by each arequalitatively different. Accordingly, the Lindblad equation with {Lx,y,z} consideredin the last section is not a cosmetic change of the Lindblad equation with {L±,z}, andthe difference between these two are crucial. From these consideration, we concludethat the frequency shift in FID oscillation period due to γx , γy originates from thenon-commutativity, Eq.(24) and that the shift is dynamic frequency shift [20, 21].The Lindblad form of Eq.(18), which was phenomenologically introduced, remindsus similarity with a relaxation process of depolarization, but it is not exactly thecase. Usually, the process is described by Lindblad master equation with d0 = 0and γx = γy = γz in Eq.(18) [25, 24]. Possible generalization to address the case ofd0 , 0 and unequality of three relaxation constants was discussed in Refs. [19, 26]. Inthese studies, microscopic derivation was performed of the Lindblad master equation todescribe generalized depolarization, starting from a model Hamiltonian consisting of atwo-level system and three independent boson bathes [26]. However in these studies,γx = γy was finally obtained. This result still survives even when coupling constantsof the two-level system and each of three boson bathes with Ohmic spectral functionsare different, and thereby we call Eq.(18) phenomenological. Since there is a crucialdifference between γx = γy and γx , γy cases, as shown above, construction of amicroscopic model system which yields the Lindblad form with γx , γy through aconventional manner, is in need as a future work from the view point of open quantumphysics.3.2 overdamped fluctuating fieldIn discussions so far, a condition 4d0 > |γx −γy | is throughout assumed. The followingtwo subsections mention the results when 0 < 4d0 ≤ |γx − γy |, as a complement.Since formulation in Sec.2.1 maintains generality to some extents, it is applicableto any case equivalent to a two-level system in a fluctuating magnetic field, like su-perconducting charge qubit[27], and ion trap [28], to name a few. Let the scheme beapplied to a certain system where 0 < 4d0 < |γx − γy | is somehow satisfied. Then thesolutions, x̃(t) and ỹ(t), are now given as(x̃(t)ỹ(t))= exp[−Ñxyt] ( x̃0ỹ0), (33)whereexp[−Ñxyt] = exp(−γall + Γ4t)Ñ (1)xy + exp(−γall − Γ4t)Ñ (2)xy , (34)Ñ (1)xy =12Γ(Γ − (γx − γy) 4d0−4d0 Γ + (γx − γy)), (35)Ñ (2)xy =12Γ(Γ + (γx − γy) −4d04d0 Γ − (γx − γy)), (36)9and a real constant Γ =√(γx − γy)2 − 16d20 . These solutions show monotonic ex-ponential decays without oscillation, which correspond to an overdamped case in aclassical harmonic oscillator. Marked contrast to the case in Sec.2.1 can be seen in amatrix property: determinants of both Ñ (1),(2)xy are null, while those of N (1),(2)xy in Eq.(13)are identically unity. Thus after applying π/2 pulse and set such an initial condition(x0, y0, 0) as [Γ − (γx − γy)]x0 + 4dy0 = 0, (37)or [Γ + (γx − γy)]x0 − 4dy0 = 0, (38)the motion of the Bloch vector, (x̃(t), ỹ(t), z̃(t) = 0), is determined by a solely Ñ (2)xy termwhen Eq.(37) is satisfied since Ñ (1)xy (x0, y0)T = 0, or by solely Ñ (1)xy term when Eq.(38)is satisfied, since Ñ (2)xy (x0, y0)T = 0, where the superscript T denotes transpose of thecolumn vector. In the former case, the longest relaxation time in this system, apart fromwhether this should be termed T∗2 or not, is expected [15].3.3 possible exceptional pointsIn addition to the algebraic difference between Lindblad equations with {Lx,y,z} and{L±,z}, there is another qualitative difference in matrices MHxyz +MLxyz and MH±z +ML±z ,under the condition d0 , 0. While both matrices have rank three in common, due toconstraint of Trρ = 1, a striking contrast is that MHxyz + MLxyz can have an exceptionalpoint, but MH±z + ML±z not. Indeed, four eigenvalues of MHxyz + MLxyz are{0, −γx + γy2, −14(γall ± iΩ)}, (39)while, those of MH±z + ML±z are{0, −γ+ + γ−4, −γ+ + γ− + 4γz8± id0}. (40)For the former, when Ω = 0, or 4d0 = |γx − γy |, the eigenvalues degenerate and, atthe same time, corresponding eigenvectors coalease. This corresponds to a criticaldamping condition for a simple classical harmonic oscillator. On the other hand, in thelater case it will not happen. The role of exceptional points in Lindblad type masterequation recently draws considerable attention [29, 30]. Studies along this line shouldbe of interest.4 ConclusionWe theoretically studied dynamic frequency shift which is expected to be observedin NV− center in diamond through FID measurement. Our analytical method was10primarily based on an extension of the procedure in Ref.[14], and novelty over it was anintroduction of anisotropic in-plane vacuum fluctuation. The result was corroboratedby a master equation with phenomenological Lindblad form, which corresponds tofull-fledged depolarization process. The origin for the shift to manifest itself lies in thenon-commutative property, whose explicit presentation is given in terms of the masterequation. The findings in this work would stimulate experimental researchers in NVcenter community and also pose such a fundamental question as a microscopic modelconstruction for the interesting Lindblad form.The author thanks T. Teraji, Y. Masuyama, C. Shinei and T. Yamaguchi for theirinputs on NV− from an experimental side, and M. Arai for his valuable comments.References[1] J. R. Maze, P. L. Stanwix, J. S. Hodges, S. Hong, J. M. Taylor, P. Cappellaro, L.Jiang, M. V. G. Dutt, E. Togan, A. S. Zibrov, A. Yacoby, R. L. Walsworth, and M.D. Lukin, Nature 455, 644 (2008).[2] G. Balasubramanian, I. Y. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin,C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, T. Hanke, A. Leitenstorfer, R.Bratschitsch, F. Jelezko, and J. Wrachtrup, Nature 455, 648 (2008).[3] C. L. Degen, Appl. Phys. Lett. 92, 243111 (2008).[4] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R.Hemmer, A.Yacoby, R. Walsworth, and M. D. Lukin, Nat. Phys. 4, 810 (2008).[5] J.-P. Tetienne, N. Dontschuk, D. A. Broadway, A. Stacey, D. A. Simpson, and L.C. L. Hollenberg, Sci. Adv. 3, e1602429 (2017).[6] M. J. H. Ku, T. X. Zhou, Q. Li, Y. J. Shin, J. K. Shi, C. Burch, L. E. Anderson, A. T.Pierce, Y. Xie, A. Hamo, U. Vool, H. Zhang, F. Casola, T. Taniguchi, K. Watanabe,M. M. Fogler, P. Kim, A. Yacoby, and R. L. Walsworth, Nature 583, 537 (2020).[7] V. M. Acosta, E. Bauch, M. P. Ledbetter, A. Waxman, L.-S. Bouchard, and D.Budker, Phys. Rev. Lett. 104, 070801 (2010).[8] D. M. Toyli, F. Charles, D. J. Christle, V. V. Dobrovitski, and D. D. Awschalom,Proc. Natl. Acad. Sci. U.S.A. 110, 8417 (2013).[9] P. Neumann, I. Jakobi, F. Dolde, C. Burk, R. Reuter, G. Waldherr, J. Honert, T.Wolf, A. Brunner, and J. H. Shim, Nano Lett. 13, 2738 (2013).[10] S. Nishimura, T. Kobayashi, D. Sasaki, T. Tsuji, T. Iwasaki, M. Hatano, K. Sasaki,and K. Kobayashi, Appl. Phys. Lett. 123, 112603 (2023).[11] E. D. Herbschleb, H. Kato, Y. Maruyama, T. Danjo, T. Makino, S. Yamasaki, I.Ohki, K. Hayashi, H. Morishita, M. Fujiwara and N. Mizuochi, Nat. Commun. 10,3766 (2019).11[12] H.-P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (OxfordUniversity Press, Oxford, 2002).[13] A. Rivas and S. F. Huelga, Open Quantum System (Springer, New York, 2014).[14] P. Jamonneau, M. Lesik, J. P. Tetienne, I. Alvizu, L. Mayer, A. Dréau, S. Kosen,J.-F. Roch, S. Pezzagna, J. Meijer, T. Teraji, Y. Kubo, P. Bertet, J. R. Maze, and V.Jacques, Phys. Rev. B 93, 024305 (2016) and summplement.[15] M. Berry, Ann. N. Y. 755, 303 (1995).[16] de Sousa, in Electron Spin Resonance and Related Phenomena in Low-Dimensional Structures (Springer, New York, 2009).[17] C. P. Slichter, The Principle of Magnetic Resonances (Springer, New York, 1989).[18] A. Yoshimori and J. Korringa Phys. Rev. 128, 1054 (1962).[19] J. Denis and J. Martin, Phys. Rev. Reserch, 4, 013178 (2022),[20] P. S. Hubbard, Phys. Rev. 128, 650 (1962).[21] L. G. Werbelow, in Encyclopedia of Magnetic Resonance (John Wiley & Sons,Germany, 2011).[22] H. Carmichael, An open systems approach to quantum optics (Springer, New York,1993).[23] F. Iachello, Lie Algebras and Applications (Springer, New York 2015).[24] A. Rivas and A. Luis, Phys. Rev. A 88, 052120 (2013).[25] M. A. Nielsen and I. L. Chuang, Quantum Computation and Qunatum Information(Cambridge University Press, 2000).[26] M. Arsenijević, J. Jeknić-Dugić, and M. Dugić, J. Phys. 47, 339 (2017).[27] Y. Nakamura, Y. A. Pashkin, and J. S. Tsai, Nature 398, 786 (1999).[28] F. S.Kaler, H. HÃďffner, M. Riebe, S. Gulde, G. P. T. Lancaster, T. Deuschle, C.Becher, C. F. Roos, J. Eschner and R. Blatt, Nature 422, 408 (2003).[29] M. S. Sarandy and D. A. Lidar, Phys. Rev. A 71, 012331 (2005).[30] J. Larson, E. Sjoqvist and P. Ohberg, Conical Intersections in Physics (Springer,New York, 2020)12