Interplay between charge and spin noise in the near-surface theory of decoherence and relaxation of $C_{3v}$ symmetry qutrit spin-1 centers

Decoherence and relaxation of solid-state defect qutrits near a crystal surface, where they are commonly used as quantum sensors, originates from charge and magnetic field noise. A complete theory requires a formalism for decoherence and relaxation that includes all Hamiltonian terms allowed by the defect's point-group symmetry. This formalism, presented here for the $C_{3v}$ symmetry of a spin-1 defect in a diamond, silicon cardide, or similar host, relies on a Lindblad dynamical equation and clarifies the relative contributions of charge and spin noise to relaxation and decoherence, along with their dependence on the defect spin's depth and resonant frequencies. The calculations agree with the experimental measurements of Sangtawesin $\textit{et al.}$, Phys. Rev. X $\textbf{9}$, 031052 (2019) and point to an unexpected importance of charge noise.


I. INTRODUCTION
Coupling of a spin-1 center in a solid, usually associated with a dopant or defect, to electric and magnetic fields provides a direct method of sensing nanoscale fields [1][2][3][4][5][6][7][8][9][10][11][12][13], of tuning the optical emission linewidth for optically-active defects [14][15][16], and of coupling to electric or magnetic excitations to realize hybrid quantum coherent systems [17][18][19][20][21][22][23][24][25][26][27][28][29].Conversely this coupling also makes the defect spin dynamics very susceptible to charge and magnetic noise, contributing to decoherence and relaxation of the spin qubit states [2,[30][31][32][33][34][35][36][37][38][39][40][41][42][43][44], and increasing the photoluminescence linewidth [14][15][16]45].Many approaches have been explored to diminish the effect of charge noise on defects, e.g., controlling the termination of the diamond surface [31], embedding diamonds in materials with a high dielectric constant [33], covering the diamond surface with an extra layer [34], and placing the spin center in the depletion region of a p-n diode [15,45].Nevertheless, surfaces present a useful laboratory for the study of noise sources, as the nature of these fluctuations can be specific to the surface type and also to the depth below the surface.For spins acting as quantum sensors for nanoscale fields, the surface noise limits how near to the surface a spin can be placed while still retaining experimentally resolvable coherent dynamics, and thus the spatial resolution achievable with the sensor.Thus a complete formalism for decoherence and relaxation will permit the surface properties to be optimized within practical parameters and will enable the optimal depth of a spin to be determined when it is acting as a quantum sensor for nanoscale fields.
In this work we provide a complete quantitative theory for the influence of the charge and magnetic noise on the dephasing and relaxation processes of the three states (|T + , |T 0 , and |T − , corresponding to the spin-1 projection along the symmetry axis) of shallow spin-1 (qutrit) solid-state spin centers with C 3v symmetry, embedded in hosts such as diamond and silicon carbide (Fig. 1).Our work includes all electric (d ⊥ E z , d E ± and d E ± ) and magnetic (γ ⊥ B ± and γ B ± ) terms allowed by symmetry (See Fig. 1).We derive a Lindblad dynamical equation [46] for our qutrit containing eight different Lindblad operators, which captures the resulting dephasing and relaxation processes.We then calculate the population and dephasing dynamics of the spin, and describe several different regimes and scenarios for these processes, suggesting improved conditions for the experimental utilization of qutrits.For example, the threestate character of our spin-1 qutrit can improve the abil-ity to probe "flat" (frequency-independent) regions of the spectral noise density.We also show that the relaxation between the spin-defect states |T ± and |T 0 , usually attributed to magnetic noise (via the γ ⊥ B ± term), also has a significant contribution from the charge noise via the commonly-ignored dipole term d E ± .
We then apply this theory to study the surface noise arising from the fluctuations of charges and magnetic moments on the diamond surface.For both hydrogen (H) and oxygen (O) terminated diamond, the bonds between carbon atoms of the diamond and these other atoms effectively create either acceptor (hydrogen) or donor (oxygen) levels at the diamond surface [47,48].These acceptor (donor) levels are occupied by the electrons (holes) provided by nitrogen dopants (so-called P 1 centers) in the diamonds which also contain negativelycharged nitrogen-vacancy (NV − ) centers.The electrostatic effects of this charge transfer bends the host bands and creates an effective, very low mobility, surface twodimensional (2D) hole (electron) gas [13,[48][49][50][51][52][53][54].Similar effects also emerge from imperfections in the crystal termination [53].The hopping motion of trapped electrons (holes) and the charge motion within the confined 2D surface hole (electron) gas produce fluctuating electric and magnetic noise that influences our shallow defects, causing relaxation and decoherence of the spin center's quantum state.
The distinct character of the sources of charge and their fluctuations require different theoretical descriptions of their effect on the defect spin.The trapped charges can be modeled as electric dipole fluctuators [34,42], whereas the confined surface charged gas should be treated as the fluctuation of point-like charges [2,13,30,31,42] (with charge neutrality maintained by the fixed charges of donors).We derive analytical formulas for both types of fluctuating electric fields as a function of the areal density of fluctuating dipoles or fluctuating charge densities, of dipole length and defect depth.We analyze the competition between these two sources of charge noise, and compare them with their bulk noise counterparts [35,45].For completeness, we include as well the magnetic noise produced by both the fluctuations of the spins' magnetic moments and the movement of charged particles (Biot-Savart law) [35][36][37][38][39][40][41].We also identify the scenarios for which the magnetic noise dominates.Our quantitative theory for both magnetic and charge noise enables the study and analysis of the competition between electric (charge) and magnetic noise in different scenarios and environmental conditions.
Finally, combining our quantitative theory for the surface charge and magnetic noise with our complete formalism for relaxation and dephasing of a spin-defect with C 3v point-group symmetry, we report calculations of the decoherence and relaxation of the spin center as a function of the surface charge density, the defect depth and the frequency separation between the spin center's energy levels.The dependence of the decoherence and relaxation on the spin center's energy levels will thus allow us to identify the dominant source of noise if these quantities are studied as a function of magnetic field.Our results show good agreement with experimental reports of the dependence of the decoherence time on the defect depth [31].Thus we propose such studies will enable the dominant noise sources to be assigned for various surface treatments.
Section II presents the ground state Hamiltonian for a spin-1 defect with C 3v point group symmetry, along with its coupling to external electric and magnetic fields.Assuming that the magnetic and charge noise will manifest as classical electric and magnetic fields, we derive the Lindblad operators followed by the Lindblad dynamical equation for this Hamiltonian.The general expressions produced for the spin population dynamics yield the different relaxation times and decoherence times associated with the loss of information among different spin-1 subspaces.Section III focuses on the specific case of point-like and dipole charge noise, and explores the competition between these two sources.Different sources of magnetic noise are also calculated, and their contribution compared to that of charge noise.From this we clarify the charge and magnetic noise dependence of the relaxation and decoherence rates on the frequency separation between the spin center's levels.Section IV compares our theoretical findings with experimental results for the decoherence of shallow NV − centers.

II. SPIN-1 (QUTRIT) DECOHERENCE FOR C3v
POINT GROUP SYMMETRY Here we establish the general features required for a calculation of the decoherence and relaxation of the quantum state of a spin-1 qutrit due to electric and magnetic noise.We first present the complete Hamiltonian for qutrits with C 3v point group symmetry in the presence of electric and magnetic fields.We then find eight Lindblad operators that produce decoherence and relaxation of the three states of the qutrit.The dynamics of the spin in the presence of fluctuating electric and magnetic fields are obtained from a Lindblad [46] dynamical equation.Finally, we identify new features of the population dynamics of our spin-1 qutrit, including relaxation rates which have been previously interpreted as due to magnetic noise dominance, but which may be due to electric noise.

A. Spin center Hamiltonian
Semiconductor spin-1 centers with C 3v point group symmetry [2-4, 14, 55-70] do not possess inversion symmetry and therefore permit linear coupling of the spin's energy levels to an electric field (Stark effect [2, 14-16, 61, 64, 71-73]) and to strain.The spin is also coupled to a magnetic field via the Zeeman effect with an anisotropic gyromagnetic ratio.The ground state (GS) Hamiltonian, with all these terms, in the triplet basis |T − , |T 0 , |T + (where +, 0 and − are defined along the symmetry axis), for spin-1 centers with a C 3v point group symmetry is, from a group theory analysis [56,62,65,74,75] where h is Planck's constant, γ is the gyromagnetic ratio tensor, S are the triplet spin-1 matrices, is the magnetic field, {A, B} = AB + BA, D is the zero energy splitting between the triplet states |T 0 and |T ± , and d , d ⊥ and d are electric dipole constants.The z direction here corresponds to the defect symmetry axis.
The temporal fluctuations of magnetic and electric fields generate decoherence and relaxation of the quantum state of the spin.To understand the role of individual terms within the corresponding decoherence and relaxation processes, we rewrite the Hamiltonian in matrix form, where E ± = E x ± iE y and B ± = B x ± iB y .From the Hamiltonian, Eq. ( 2), the magnetic field produces a frequency split ∝ γ B z between the |T ± states, in addition to a coupling ∝ γ ⊥ B ± between the |T 0 and |T ± states.
Similarly, the electric field yields a frequency splitting ∝ d E z between |T ± and |T 0 , in addition to a coupling ∝ d E ± between the |T 0 and |T ± states.However, unlike the magnetic field, the electric field also couples the |T − and |T + subspaces, with a strength proportional to d ⊥ E ± .These different terms appear schematically in Fig. 1.As relaxation processes (1/T 1 ) occur when different levels are coupled to each other through random temporal fluctuations, the d , d ⊥ and γ ⊥ terms will contribute to relaxation processes.Conversely, the dephasing processes can also occur due to any terms responsible for relative fluctuations of the energy of the levels, namely d and γ in addition to d , d ⊥ and γ ⊥ .We stress that although prior work has neglected the presence of the d electric dipole terms within the spin center's Hamiltonian, these are important when charge noise dominates.They are also important to characterize correctly processes involving photoluminescence and spin dynamics near the level anticrossing of the electronic ground state (GSLAC) [76][77][78][79][80][81][82][83][84][85], as well as for acoustical driving experiments of the |T 0 ↔ |T ± spin transition [75,86].Moreover, although up to this point there is no precise experimental verification for the value of d , Ref. [86]  To obtain the Lindblad dynamical equation describing decoherence and relaxation of our qutrit we begin by considering that the only nonfluctuating external field is the dc magnetic field B z , which controls the frequency separation between |T + and |T − subspaces [See Fig. 1].The other fields can fluctuate, so we rewrite our Hamiltonian, Eq. ( 2), as the sum of a time independent part and a time-dependent one, i.e., H = H 0 + V(t) with z − 2/3), and V(t) produced by the remaining terms of Eq. ( 2).In the absence of these (weak, relative to the value of D) fluctuating fields the spin center's frequencies are ω 0 = 0 and ω ± /2π = D ± γ B z .
In order to solve the dynamics of the spin center, we move to the interaction picture with respect to H 0 , which leads to where |ψ (t) I ≡ e i H 0 t |ψ (t) and The evolution of the density matrix in the interaction picture associated with Eq. ( 3) is where T is the time-ordering operator.The solution to Eq. ( 5) can be obtained through a perturbative (Dyson) expansion of the propagator Π(t, t 0 ), namely, with Π = Π 0 + ΠΣΠ 0 where Σ is the self-energy.This yields the general dynamical equation for ρI (t) [87,88] dρ Mapping this equation onto a Lindblad equation [46] requires some assumptions and approximations.The equation is expanded in a diagrammatic perturbation series to second order, and the average is taken over different realizations, namely, • • • .We further assume the correlation time of the noise, τ c , is much smaller than the time interval t−t 0 , i.e., t−t 0 τ c .A further Markovian approximation yields the result The identity is of great use; the first right hand side term produces new contributions to the coherent evolution e.g., Stark shifts and Lamb shifts, however the remaining ones produce the terms associated with Lindblad operators.Defining The rotating wave approximation with respect to the frequency separations between the spin center's energy levels, ω µν = ω µ − ω ν with µ, ν = {0, ±} simplifies Eq. (10).This approximation cannot be employed if two of the states become degenerate, and hence the results are only valid for lifted degeneracies.We use Novikov's theorem [89][90][91][92][93] together with the weak coupling between our spin center states and the fluctuating fields [90][91][92][93].We further assume temporal translational symmetry for the fluctuating fields (stationary regime), y, z, which follows for fluctuating fields lacking a preferential direction.The corresponding noise spectral densities are where due to the classical character of our fluctuating fields we have S E(B) (ω) = S E(B) (−ω).This holds for For the NV − the largest frequency split is ≈ 2.5 GHz, so the approximation holds for T 1 K. Another consequence of the frequency-symmetric noise spectral density is the absence of an effective coherent Hamiltonian arising from the noise, i.e., no Stark nor Lamb shift, so H eff,I = 0. Three rates are usefully associated with the charge noise spectral density, namely . Two rates are correspondingly associated with the magnetic noise spectral densities, with γ ⊥ = γ⊥ /2π and γ = γ /2π.Finally, all the considerations above yield the Lindblad dynamical equation with Lindblad operators in the interaction picture, L k,I , given by Here the operator S z S + (S − S z ) represents the raising (lowering) operator within the subspace spanned by {|T + , |T 0 }, while S + S z (S z S − ) represents the raising (lowering) operator within the subspace spanned by {|T 0 , |T − }.Additionally, the operator S 2 + (S 2 − ) is the raising (lowering) operator within the subspaced spanned by {|T + , |T − }.
Using the Lindblad operators [Eqs.( 19)-( 26)] within the Lindblad equation [Eq.( 18)], we can also obtain the following differential equation that governs the dynamics of the density matrix ρI (t) µν = ρ µν (t), namely, with µ, ν = {0, ±} and the corresponding Lindbladian or Liouvillian matrix, L 9×9 .For our case, L 9×9 is composed of a 3 × 3 block diagonal matrix that governs the relaxation process of our quantum states, and a diagonal 6 × 6 matrix governing the dephasings between different subspaces.Both processes will be investigated in the next two subsections.

Spin center relaxation
The part of the Lindbladian governing the relaxation process is described by the evolution of the diagonal elements of ρI (t) , namely where we define Γ γd (ω) = Γ d (ω)+Γ γ ⊥ (ω).The solution of this equation is obtained through the ansatz Here, λ i are the three eigenvalues of the 3 × 3 matrix within Eq. ( 28), and a + i a 0 i a − i T are the corresponding eigenvectors.In principle, the three eigenvalues (λ i ) define three different relaxation rates, namely, with ), defined similarly to Refs.[38,42].Accordingly, the general solution for the evolution of the diagonal density matrix elements is Using Tr[ ρI (t) ] = 1 we obtain c 1 = 1/3, and the remaining coefficients c 2,3 are determined by the initial condition of the density matrix.For the initial condition ρ −− (t = 0) = 1, and Interestingly, and differently from the spin-1/2 case, here we obtain a bi-exponential relaxation.This biexponential behavior was already discussed and presented in the literature for spin center dephasing [45] and relaxation processes [38,42].
Now we analyze some cases of particular interest.The first one corresponds to ω +0 ≈ ω −0 ≡ ω ±0 , which produces Γ γd (ω +0 ) ≈ Γ γd (ω −0 ).Accordingly, this leads to Ω + = Ω − ≡ Ω, which was already discussed in the literature [42] and produces the relaxation times Importantly and contrary to Tetienne et al. [38] and Myers et al. [42], we can see that when we consider all symmetry-allowed terms within our Hamiltonian Eq. ( 2), the charge noise also contributes to the Ω rate.Therefore, this rate does have an exclusive magnetic noise origin.Associating Ω to the magnetic noise is only accurate for Γ γ ⊥ (ω ±0 ) Γ d (ω ±0 ).The identification of a significant role for electric field noise in this rates can only be obtained due to the inclusion of the d dipole terms allowed within the Hamiltonian for spin centers with C 3v symmetries.
Experimental comparison of the measured rates, Ω + and Ω − , indicates whether the noise spectral density is constant within the ω −0 < ω < ω +0 region.As a consequence, the difference between the nominal values of Ω + and Ω − can be used to obtain information about the flatness of any noise spectral density.Hence the spin-1 relaxation mechanisms of C 3v spin-defects can also be used to probe the presence of flat regions of the spectral noise density.
In Fig. 2, we plot the three level population dynamics of our spin center as a function of time.We use three different regimes of parameters with the same initial condition, ρ −− (t = 0) = 1.In Fig. 2(a), we have the corresponding dynamics for γ Ω ± .Although this regime is usually attributed to the dominance of the charge noise, the rate Ω ± also contains charge noise contributions via Γ d (ω ±0 ).As a consequence, this attribution is only accurate for Γ d (ω ±0 ) Γ d ⊥ (ω +− ).Since different works suggest d ≈ d ⊥ [65,86], this translates to small values of magnetic fields satisfying ω +− ω ±0 for S Ei (ω) ∝ 1/ω α .Nevertheless, in this case we see that a system initially prepared in the |T − state will start to increase the population of its |T + state due to the faster γ relaxation rate.The corresponding relaxation is characterized by the T + 1 timescale, Eq. ( 30), plotted as a vertical gray line.We see that in this process, the system first reaches an approximately equal population of |T − and |T + .After this, the Ω ± rate starts playing a role, and we see an increase in the population of |T 0 , with a characteristic timescale given by T − 1 , Eq. ( 31) (vertical gray line).Finally, for long times t T ± 1 , we obtain an equal population for all the levels.Similarly, the dotdashed lines represent the same process but for Ω − = Ω + .
In Fig. 2(b), the solid lines represent the plot for the population dynamics corresponding to γ = Ω ± = Ω.Here both the initially unpopulated levels, |T 0 and |T + , experience a population increase with equal rate, and with a characteristic timescale defined by T ± 1 = (3Ω) −1 [Eqs.(30) and ( 31)] (see vertical gray line).As a consequence, the populations of the states |T 0 and |T + in-crease equally.Similarly, the dot-dashed lines represent the same process but for γ = Ω + < Ω − .As a consequence of Ω − > γ, we see that the population of |T 0 increases faster than the population of |T + .
Lastly, in Fig. 2(c) we plot the population dynamics for γ Ω ± , usually associated with the dominance of magnetic noise.However, as Ω ± = Γ γd (ω ±0 ) also contains a charge noise contribution (since Γ γd (ω ±0 ) = Γ d (ω ±0 ) + Γ γ ⊥ (ω ±0 )) we stress this identification of magnetic noise dominance will only be accurate if Γ d (ω ±0 ) Γ γ ⊥ (ω ±0 ).We first see Ω − as responsible for the increase of population of |T 0 , with characteristic timescale T + 1 .At this point, as we have a finite population of |T 0 , the rate Ω + will begin to contribute, increasing the population of |T + .Depending on the relative value of γ and Ω + , γ can also contribute to the increase of the |T + population.19)-( 26)] lead to three different dephasing rates within the three corresponding subspaces {|T 0 , |T + }, {|T 0 , |T − }, and with corresponding dephasing times We stress that this is the first work to report the accurate expression for the decoherence times including all the symmetry-allowed fluctuating terms.We see from these expressions that dephasing within {|T µ , |T ν } is not solely originating from fluctuating terms within the same subspace.For example, even though d E ± (t) does not appear within the {|T − , |T + } subspace, the decoherence times T 0± 2 depend on d E ± (t).This shows that an indirect loss of coherence between coupled subspaces also happens.The same feature also occurs for the {|T 0 , |T ± } subspace, and was already discussed in Ref. [45].In short, this shows the importance of taking into account fluctuators over the whole spin-defect manifold when calculating dephasing times.

III. THEORY OF THE FLUCTUATING ELECTRIC AND MAGNETIC FIELDS
Unintentional impurities within crystals can either donate electrons or accept electrons, leading to free electrons or holes in the crystal.As already discussed in Refs [13,42,45], these particles do not distribute uniformly and are also non-static, due to the thermal fluctuations of the electron and hole position, collisions between them, continual absorption and release by donors or acceptors, among other processes.Moreover, as the read-out and initialization of the spin center's state are performed with laser illumination, the measurement process additionally agitates the particles, increasing these fluctuations.As different types of bulk noise were already discussed in Ref. [35,45], here we calculate the fluctuating surface fields arising from all the possible types of charge [2,[30][31][32][33][34][35]42] and magnetic noise [36][37][38][39][40][41][42].We also provide an analysis for the competition between surface and bulk contributions.We emphasize here that surface charge noise, δE(t), can occur due to: (1) electrons or holes trapped at the crystal surface, giving rise to a fluctuating dipole electric field; (2) confined hole or electron gases produced by band bending near the surface [48,[51][52][53], which produces a point-like fluctuating electric field; and (3) electrons that are excited to the conduction band and therefore also contribute to a pointlike fluctuating electric field.We calculate the magnetic noise arising from fluctuating magnetic moments at the surface, in addition to the magnetic noise produced by a random movement of charged particles (Biot-Savart law).For both charge and magnetic noise, we analyze the competition between their bulk and surface counterparts.

A. Coordinate system
The z axis used in both Eqs. ( 1) and ( 2) is defined with respect to the main symmetry axis of our defect spin center.This axis, however, can assume different directions with respect to the surface normal vector n.Examples include diamond with surfaces perpendicular to either [111] or [001] crystalographic axis [94], and for the different orientation of divacancies in SiC [45].For general results we keep an arbitrary direction of the spin center main axis with respect to the surface normal, defined by n = ẑ .Accordingly, for the spin center and the surface, we have axis {x, ŷ, ẑ}, and {x , ŷ , ẑ }, respectively.They are related to each other via a rotation of θ around the x axis, R x (θ), This distinction is important as the fluctuations of both charges and magnetic moments at the surface produce a large fluctuating field along ẑ.Accordingly, we define our surface as S = {(x , y , z ), −L/2 ≤ x , y ≤ L/2, z = 0}, and our defect position as r def = (0, 0, −z def ) = −z def ẑ .

B. Analytical calculation for the surface fluctuating point-like electric field
We now calculate the electric field correlation function E p,µ (t) E p,µ (0) arising from fluctuating surface pointlike charges [2,30,31,42].First, for the point-like charges we assume the electric field at r = r def produced by the i'th point-like charge is where Q i are both positive and negative trapped surface charges localized at r p i ≈ x p i , y p i , 0 with R p i = r p i − r def and R p i = |R p i |.Assuming now a total number of point-like charges given by N p , the total electric field experienced by the defect is , and the correlation is where we assume there is no correlation between electric fields produced by different point-like charges.The inclusion of the non-null correlation between different particles was extensively analyzed in Ref. [13].We can evaluate this expression by assuming a continuous probability distribution for the positions r i , p S = p S (r ), yielding where S is the surface containing the fluctuating charges, and n S (r ) = N p p S (r ) is the surface density of the particles i.Using an uniform probability distribution for the positions of the point-like charges, we obtain yielding These equations were obtained assuming the maximum value for the correlation of the fluctuating charge positions.For this type of noise, the frequency dependence of spectral noise density can be obtained from many different processess, e.g., absorption-release of electrons and holes by different traps and diffusion of electron and holes within our crystal.For those, the corresponding spectral noise density can be obtained similarly to Ref. [45], yielding We now assess an important characteristic of our system, namely, how does the surface charge noise compete with the bulk noise arising from the fluctuating charges within the bulk of the sample [45].This becomes important when spin centers are used for sensing, which requires maximizing the signal-to-noise ratio.The depth at which both contributions are nearly equal sets the optimal defect depth, z opt .To perform this analysis, we assume that the bulk noise comes predominantly from the near noise contribution of Ref. [45], namely where is the volume density of fluctuating charges and Assuming that the surface noise arises from point-like fluctuating charges, Eqs. ( 50)-( 52), we obtain the optimal defect depth, In Figs.Here, the solid lines represent the surface charge noise and the dot-dashed ones represent the bulk near noise.In Fig. 3(a) we plot both contributions for typical surface densities n S = 10 11 , 10 12 and 10 13 cm −2 together with bulk densities n V = 10 14 , 10 15 and 10 16 cm −3 .While the surface noise contribution shows a 1/z 2 def depth dependence, the bulk one shows no dependence.Assuming an approximately bulk density of fluctuators of 10 15 cm −3 , the surface noise always dominates, thus showing its critical importance for shallow defect implantation.Additionally, to have a better understanding of this competition, in Fig. 3(b) we also plot both contributions as a function of the surface density of point-like fluctuators for different defect depths, 5, 20 and 50 nm.We observe that for very shallow defects with z def ≈ 5 nm, the surface noise will dominate even for low surface densities n S ≈ 10 10 cm −2 .On the other hand, for z def = 50 nm and assuming n V = 10 15 cm −3 , we obtain a dominance of the surface noise only for n S > 10 12 cm −2 .

C. Analytical calculation for the surface fluctuating dipole electric field
Here we calculate the electric field correlation function E µ (t)E µ (0) arising from fluctuating surface dipole Point-like surface noise vs bulk near noise Dipole surface noise vs bulk near noise  50)-( 52)] and bulk near noise [Eq.( 54)] as a function of depth for different surface densities (a), and as a function of surface density for different defect depths (b).Similar to (a) and (b), but representing the competition between fluctuating electric field for dipole surface noise [Eqs.( 60)-( 62)] and bulk near noise [Eq.( 54)].In all the graphs, the fluctuating electric field due to the bulk near noise Eq. ( 54) is plotted for nV = 10 14 , 10 15 and 10 16 cm −3 charges [34,42].Hence, the electric field arising from the i'th dipole, corresponding to the displacement of charges Q i and −Q i located at r d i ≈ x d i , y d i , 0 and separated by the dipole distance, d i , is written as with Assuming a total number of dipoles given by N d , the total electric field experienced by our defect is E =

and we find
We first assume that the electric field produced by two different dipoles are not correlated, i.e., 0) .Secondly, we assume that the charge dipole displacements d i are randomly distributed along x , ŷ and ẑ with equal probability [45] and amplitude, yielding d i (t)d i (0) = d(t)d(0) for i = x, y, z.For large number of dipoles, the above expression can be calculated by assuming a continuous density for the discrete dipole charge positions r i , n S = n S (r), Finally, we can write For a position independent areal density n S (r) = n d S and considering both surface lengths much larger than |r def | = z def , we obtain From Eqs. ( 60)-( 62) and Eqs. ( 50)-( 52), we see that the point-like and dipole field contributions yield different dependences on the defect depth z def .Whereas for pointlike charge fluctuations we obtain a dependence z −2 def , the dipole one gives a weaker dependence z −4  def .This different defect depth dependence can be verified experimentally, and has already been discussed in the literature [40,42].In the last section of this paper we will discuss how different surface treatments result in different types of charge noise (point-like or dipole).
Here we also obtain the optimal defect depth that defines the depth in which the bulk near noise becomes comparable to the surface noise, namely Assuming a typical value dS = 0.5 nm, we plot in Fig. 3(c) and (d) the competition between the surface charge noise arising from the fluctuation of dipole charges, Eqs. ( 60)-( 62), and the bulk near noise Eq. ( 54).Due to the weaker character of the surface dipole charge noise, we can see that for a bulk density of fluctuators between 10 14 and 10 16 cm −3 and n d S = 10 11 , 10 12 and 10 13 cm −2 , the optimal depth is always smaller than 20 nm, showing the large contribution of the bulk noise to the dephasing and relaxation of shallow spin centers.
In Figs.3(c) and (d) the calculation of surface dipole charge noise contribution was done using dS = 0.5 nm.This, however, ignored the temporal dependence of the fluctuating dipoles.A more precise and accurate analysis includes the power spectral density of the dipole, yielding with Although it is still necessary to know and characterize S d (ω), different studies have already investigated that in a similar context [95].Here, we first assume that the electric dipole will be represented by a trapped charge at the surface, which behaves as a harmonic oscillator with frequency ω d and mass m * .Using the equipartition theorem, i.e., , we obtain The frequency ω d is commonly associated with the energy difference between two traps (E 0 and E 1 ), in which the charge can be imprisoned or trapped, yielding [95,96].Single dipole fluctuators are usually described by a Lorenztian spectral noise density, namely where Γ d is the rate associated with the transition between energy levels E 0 and E 1 .If we assume this rate is dominated by thermal activation, the rate will be given by Γ d ∝ e −(E1−E0)/k B T .When the temperature is not high enough to thermally activate the high-energy levels, Γ d is attributed to trapped electrons that diffuse at the surface through quantum tunneling over an energy barrier, E b , with a corresponding width b.For the case of a double-well potential described by a parabolic onedimensional potential, [97,98] with ω d = 2E b /m * b 2 assuming hopping between two 1D parabolic double wells with an effective electronic mass m * .

D. Competition between point-like and dipole surface charge noise
Here we analyze the competition between the charge noise due to point-like and dipole fluctuations, Eqs. ( 50)- (52) with Eq. ( 53), and Eqs. ( 60)- (62) with Eq. ( 68), respectively.The ratio defines the condition driving this competition.The characteristic lengths (dipole length and defect depth) and densities of these different mechanisms play a role in this competition, however the ratio of the correlation time of the fluctuations τ p and Γ −1 d also has a strong effect.If we assume that both fluctuations originate from the same physical mechanisms, we can assume τ p ≈ Γ −1 d .For this condition, the right-hand side of the equation above becomes n S z 2 def /n d S d 2 S .Hence, if the number of fluctuating point-like charges is approximately the same as the number of dipole fluctuators, i.e., n S ≈ n d S , the decoherence due to point-like fluctuations will dominate for z def dS , which is usually the case as z def 5 nm and dS 1 nm.Although this would imply that the fluctuation of the point-like charges are always the dominant ones, this result relies exclusively on τ p ≈ 1/Γ d and n S ≈ n d S .If we think of the fluctuating dipole charges as being described by trapped electrons (holes) due to surface acceptors (donors), and the point-like fluctuating charges as being produced by the surface hole (electron) gas produced by the electrons (holes) imprisonment, this would yield n S ≈ n d S .However, as different experiments have already shown the dominant dipole character of the surface noise, the trapped charges at the surface must be obtained either by the treatment of the surface crystal or its contact with the atmosphere [31].

E. Frequency dependence of the electric spectral noise density
In both previous two sections, we have used the fact that the frequency dependence of the spectral noise density is given by the Lorenztian Generally speaking, this follows when the dynamics of the fluctuation is characterized by only one characteristic time given by τ .As a consequence, we have a correlation function for the electric field produced by the fluctuators scaling as e −t/τ , whose Fourier transform yields Eq. ( 71).Although this is not incorrect, this form relies on a key assumption: that either dipole or point-like charges have the same timescale associated to their fluctuations (τ ).This is not typically the case, and as a consequence, the spectral noise density deviates from S 0 Ei (ω) in realistic physical situations.For instance, recent experiments on shallow NV centers suggest a S Ei (ω) ∝ ω −1 dependence [31,33,42].There are many studies [99][100][101] explaining the origin of the ω −1 spectral noise dependence.Here, we will use a similar approach, and apply it directly to our case.We first assume that τ has an origin in activation processes, either due to the continual trapping and release of electrons and holes, or due to the diffusion of them.We then assume τ = τ 0 e E/k B T where E is the energy associated with either the tunneling between different trap centers or the activation energy.Due to roughness of the surface crystal and the range of trap centers, we cannot assume only one particular value for E but rather a distribution for it, and here we take E 1 < E < E 2 .Taking this into account, our spectral density noise becomes with g(E) the density of states of τ with respect to E, and P (E) the weight associated to S 0 Ei (ω, E).By requiring ∞ 0 dEg (E) P (E) = 1 we obtain which yields the following three different frequency dependences which not only captures the ω −1 dependence but also shows a ω −2 dependence at higher frequencies.

F. Charges and dipoles at the interface
Within the last sections, we have assumed that all of our particles were always within the crystal region, thus experiencing a dielectric constant .Conversely, if these charges and dipoles are now placed right at the interface between our material and the external medium, we need to perform the following change in the equations above, where ext is the dielectric constant of the external environment.Hence, for defect spins embedded within a high dielectric constant material a reduction of the charge noise is expected, which was already verified through an enhancement of the spin coherence time [33].

G. Results for charge noise
In Fig. 4 we study the response of the rates (and associated relaxation times) generated by the charge noise as a function of the magnetic field (B z ) controlling the separation between the spin-defect energy levels [See upper plot of Figs.13)-( 15)], assuming the charge noise is dominated by fluctuations of point-like charges [Eq.(53)].First, we observe that the rate Γ d ⊥ (ω +− ) (yellow curve) decreases as a function of B z .This happens because the larger the B z , the larger the frequency separation between |T − and |T + (ω +− ), thus suppressing the relaxation rate between these states.Similarly, we see the same behavior for Γ d (ω +0 ), represented by the blue dotdashed line.Conversely, the rate Γ d (ω −0 ), indicated by the red plot, increases as a function of B z .This happens because the frequency separation between |T + and |T 0 decreases when we increase the magnetic field, until it reaches the degeneracy point defined by B c ≈ 0.105 T, where it starts to increase again.When the frequencies ω − and ω 0 becomes degenerate, Γ d (ω −0 ) reaches its maximum.In addition, we also plot the relaxation times 1/T 1± [Eqs.(30) and (31), here gray and black solid lines, respectively], with 1/T 1+ resembling a sum of all the rates previously discussed.Interestingly, the response of our levels to the magnetic field causes a non-monotonic behavior of 1/T 1+ , which was already discussed in the literature for different systems [102].Finally, we see that the only rate that does not change as a function of the magnetic field is Γ d (ω = 0).This rate does not give rise to an extra relaxation process, as it does not couple two different levels.Therefore, it will not depend on the frequency difference between the energy levels, and will only cause decoherence.
Similarly, in Fig. 4(c) we plot the same rates as a function of B z , but assuming the charge noise is dominated by a fluctuation of dipole charges.These are now calculated with Eqs. ( 60)-( 62) with a corresponding frequency dependent spectral noise density given by Eq. (68).For the chosen parameters, the charge noise arising from the fluctuation of point-like charges is stronger, thus generating faster relaxation rates.
Most importantly, this whole analysis shows that to achieve the maximum suppression of the charge noise in NV center spin defects, we must choose magnetic fields values B z ≈ 0.06 T as this represents the maximum value for T 1± .This value depends on the ratio of the dipole moments d ⊥ /d of our spin-defect, in addition to the assumption of a Lorentzian spectral noise density [Eq.(68)].

H. Analytical calculation for surface fluctuating magnetic fields
As we saw before, the fluctuation of the charged particles (electrons) creates the surface charge noise.However, in addition to the charge, electrons also contain spin, and hence a fluctuation of the orientation of the particle's intrinsic magnetic moment is also expected [35][36][37][38][39][40][41].Similarly, the magnetic field experienced by the defect at r def due to a magnetic dipole moment µ i located at r m i is given by with R m i = r m i − r def and |µ i | = hγ/2 = gµ B /2 for spin-1/2.Assuming that the density of charged particles is equal to the density of particles with magnetic moment, n S , we can also determine the temporal correlation of the magnetic field with Eq. ( 12), namely where we have assumed the µ i 's are randomly distributed along x , ŷ and ẑ with equal probability and amplitude, yielding µ i (t)µ i (0) = µ(t)µ(0) for i = x, y, z.Using the definition of power spectral density, we obtain Similarly to the result obtained in Eqs. ( 64)-( 66), here we also obtain a z −4 def dependence on the defect depth due to the dipole character of Eq. (76).Moreover, if we treat the magnetic moments as being predominantly due to particles with spin-1/2, the characteristic frequency of the system will be set by ∆ω µ = ω + µ − ω − µ with ω ± µ = ±(1/2)gµ B B = ±(1/2)γhB, and hence we can assume a power spectral density corresponding to a Lorenztian peaked at ω = ∆ω µ , Within the derivation above, the random character of the magnetic moment orientation was assumed.We emphasize this is only consistent at "high" temperatures defined by k B T gµ B B, which gives ≈ 1, where N ↑ , N ↓ are the populations of spin polarization.Hence our results are valid for temperatures T gµ B B/k B , which for g ≈ 2 and 0 ≤ B 0.1 T yields T 130 mK.For T 130 mK the spins starts to align and the magnetic noise starts to be suppressed.

I. Fluctuating magnetic field due to the movement of charged particles
In addition to the magnetic noise produced by the fluctuating magnetic moments, we also have a magnetic field noise produced by the movement of charged particles, e.g., the Johnson-Nyquist noise in metals [44,103,104].Accordingly, this can be obtained through Biot-Savart law where Q i is the charge of the i-th particle, v i is its velocity and N q is the total number of mobile charged particles.Since we are assuming an areal density of charges, our velocity components can be approximated to v i = v i,x x + v i,y ŷ , i.e., no velocity perpendicular to the surface.Similarly to the previous section, we assume no correlation between different particles' position, and a continuous probability distribution for the particles' positions n S (r), yielding where we have assumed a density n S of mobile charged particles behaving as a 2D Brownian-Drude model, and hence v µ (t) v ν (0) = δ µν v 2 ν e −t/τ for µ, ν = x , y , z , where τ is the relaxation time [100].Additionally, for Using now the definition of power spectral density, we write for µ = x, y, z, where σ 2D (ω) = n S e 2 m * τ 1+ω 2 τ 2 is the 2D conductivity.These results show a stronger dependence on the defect depth compared to the ones obtained from the fluctuation of magnetic moments, Eqs. ( 80)-( 82), and cannot be neglected.This type of noise was already studied in the literature in the context of NV-sensing the electrons within a metallic medium [44,105].

J. Competition between dipole magnetic noise and magnetic noise produced by charged particles
Here we compare the spectral noise density for the point-like and dipole surface charge noise, Eqs. ( 80)- (82) and Eq. ( 89).The ratio defines the condition driving this competition.This expression is very similar to Eq. ( 70).

K. Results for magnetic noise
Here we study and investigate the dependence of the rates (and associated relaxation times) generated by the magnetic noise as a function of the magnetic field (B z ) [See Fig. 5(a)].In Fig. 5(b), we plot (dot-dashed lines) the rates Γ γ ⊥ (ω ±0 ), and Γ γ (ω = 0) [Eqs.( 16) and ( 17)], assuming the magnetic noise is dominated by fluctuations of magnetic moments [Eqs.( 80)-( 82) and ( 83)] with ∆ω µ = 0. First, we observe that the rate Γ γ ⊥ (ω +0 ) (blue curve) decreases as a function of B z .This happens because the larger the B z , the larger the frequency separation between |T + and |T 0 (ω +0 ), thus suppressing the relaxation rate between these states.On the other hand, the rate Γ γ ⊥ (ω −0 ) (red curve) increases as a function of B z .This happens because the frequency separation between |T − and |T 0 decreases when we increase the magnetic field.Similarly to the case of charge noise, we see that the only rate that does not change as a function of the magnetic field is Γ γ (ω = 0).This rate does not give rise to an extra relaxation process, as it does not couple two different levels; thus it only causes decoherence.
In Fig. 4(c) we plot the same rates as a function of B z assuming ∆ω µ = γhB z .For this case, the maximum of the spectral noise density is achieved when the spin-defect frequency transition ω −0 matches the energy splitting associated to the magnetic moments causing the noise, i.e., ω −0 = ∆ω µ .This can be understood from the point of view of energy conservation, implying that the maximum relaxation rate happens when we have a relaxation of the spin-defect states (|T − → |T 0 ), followed by the excitation of the two-level systems defined by the spin-1/2 magnetic moments.This behavior was already seen experimentally, with magnetic spin noise given in terms of spin-1/2 nitrogen impurities [76].
times 1/T 1± [Eqs.( 30) and ( 31), gray and black solid lines, respectively].Similarly to the charge noise case, we observe a non-monotonic behavior of 1/T 1+ .Here, however, the maximum of the rates occurs for B z such that ω +0 = ∆ω µ .Finally, we also stress that the experimental B z -dependence of 1/T 1 shows the same trend as the results developed in Ref. [76].

IV. COMPARISON WITH EXPERIMENTAL DATA
In this section, we analyze and interpret the depth dependence of the decoherence times of NV-centers measured within the experimental work of Ref. [31].Using our theory developed in the previous sections, we analyze these experimental results and provide possible explanations for the different depth dependence originating from the treatment of different samples.
The experimental data of the NV-center decoherence time as a function of depth (z def ) is plotted in Fig. 6.Figs.6(a) and (b), contain the decoherence time for samples A and B, respectively, while Fig. 6(c) shows the decoherence for other samples.Many samples of Ref. [31], including sample A, present a depth dependence that is consistent with dipole fluctuations (e.g., fluctuating dipole charges or magnetic moments) [See Fig. 6(a)].This yields a T 2 ∝ z 4 def dependence as can be seen from Eqs. ( 60)-( 62) and ( 80)- (82).On the other hand, we see that for sample B, the decoherence time follows better the depth dependence related to fluctuation of point-like charges, described by Eqs. ( 50)-( 52) and ( 89), and T 2 ∝ z 2 def .Furthermore, we also see that sample B is the one with the longest coherence time for shallow NV-centers.Unlike sample A [Fig. 6(a)] and the other samples of Fig. 6(c), sample B is the only one that was subjected to high temperature annealing followed by oxygen annealing (even though it was was part of the same initial crystal as sample A).Accordingly, through the comparison between data of sample A and sample B, we understand that the roughness of the crystal surface tends to produce fluctuating dipole-like fields, which are the dominant source of noise.As sample B had its surface treated, this type of noise was suppressed, leaving only the contribution of point-like fluctuating noise due to either surface-confined electron or hole gases, or electrons at the conduction band excited during the laser readout or initialization.

V. CONCLUSION
In this work, we first present a complete theory for the decoherence and relaxation of qutrit (spin-1) spin centers with C 3v point-group symmetry.Accordingly, we obtain all the Lindblad operators arising from both charge and magnetic noise, and from that, we calculate the associated decoherence and relaxation times.We further present the relaxation dynamics for both charge and magnetic noise dominances.In the second part, we develop a microscopic theory for the charge noise arising from both point-like and dipole fluctuating charges, as well as for the magnetic noise arising from fluctuation of the magnetic moment and from randomness of the movement of charged particles.Using this quantitative theory, we study the evolution of the rates associated with both charge and magnetic noises as a function of the energy separation between the defect energy levels, which is produced by a finite magnetic field along the defect main symmetry axis.Finally, we use the theory developed in this work to analyze the depth dependence of decoherence times for samples with different treatments.60)-( 62) and ( 80)-( 82)].(b) Sample B, which has been treated with hightemperature annealing and oxygen annealing, shows T2 follows a depth-dependence that better explained with point-like fluctuations, i.e., 1/T2 ∝ 1/z 2 def [Eqs.( 50)-( 52) and ( 89)] (c) Other samples that were not annealed present T2 depth-dependence better explained by the dipole fluctuations.

FIG. 1 .
FIG. 1. Schematic plot of a spin-1 negatively-charged nitrogen-vacancy (NV − ) center within diamond in the presence of fluctuating surface charges (gray spheres) and magnetic moments (ruby arrows).(Lower right) NV − center's spin levels and response to the magnetic field.Dephasing and relaxation processes are indicated schematically by double arrows associated with the electric dipole terms d ⊥ E±, d Ez and d E± and the magnetic dipole terms γ ⊥ B± and γ Bz.
3(a) and (b) we plot the quantity |δE i | = E p i (t) E p i (0) for the surface charge noise arising from fluctuations of point-like charges at the surface, Eqs.(50)-(52), and the bulk near noise |δE b | arising from the fluctuation of point-like charges in bulk, Eq. (54).