Ab initio calculation of the spin lattice relaxation time $T_1$ for nitrogen-vacancy centers in diamond

We investigate the fundamental mechanism of spin phonon coupling in the negatively charged nitrogen vacancy center $(\mathrm{NV}^-)$ in diamond in order to calculate the spin lattice relaxation time $T_1$ and its temperature dependence from first principles. Starting from the dipolar spin-spin interaction between two electrons, we couple the spins of the electrons to the movements of the ions and end up with an effective spin-phonon interaction potential $V_{s-ph}$. Taking this time dependent potential as a perturbation of the system, a Fermi's golden rule expression for transition rates is obtained which allows to calculate the spin lattice relaxation time $T_1$. We simulate the color center with the Vienna Ab initio Simulation Package (VASP) to extract the figures necessary to quantify $T_1$. We investigate on the local phonon modes of the color center within the harmonic approximation using the small displacement method and extract the phononic density of states and bandstructure by diagonalizing the dynamical matrix using the phononpy package. We show that our model allows to calculate $T_1$ in good agreement with experimental observations.


I. INTRODUCTION
The negatively charged nitrogen vacancy center (NV − ) is an important colour center in diamond 1 that consists of a substitutional nitrogen atom adjacent to a vacant lattice site. Six electrons are located at the center, which exhibits C 3v -symmetry, and they form an electronic ground state spin triplet transforming according to the A 2 representation. This state is further split by the dipolar spin-spin interaction into a m s = 0 ground state and two degenerate excited m s = ±1 states with a zero field splitting constant D/h = 2.88 GHz. 2,3 The spin of the system can be prepared and read out optically 4 which leads to many applications in magnetometry, [5][6][7][8] biolabelling, 9 nano-sensing 10,11 and makes it a promising candidate for a solid state quantum bit. [12][13][14] Since the spin is the property to be played with in applications, a proper understanding of spin relaxation is of utmost importance. In this paper we deal with the longitudinal spin relaxation in the 3 A 2 ground state triplet caused by the interaction of the electron spins with the phonons of the crystal. Experimental studies 15,16 have suggested that the temperature dependence of the spin-lattice relaxation rate in a range between 10 K to 500 K is well described by a twophonon Raman process and an Orbach process, 17 however there are measurements where a different behaviour was observed. 18 Also, the measured relaxation rates differ by one order of magnitude for different samples. To un-derstand and predict spin-lattice relaxation times quantitatively in this system the fundamental mechanism of spin-phonon coupling has to be investigated. Insight into this coupling mechanism is most easily achieved by considering first-order processes, which are dominant at low temperatures around the spin transition energy D of the spins (T = 138 mK). At these temperatures the phonon spectrum is frozen out and thus higher-order processes are suppressed. In a recent paper 19 a direct single phonon relaxation process and spin lattice relaxation times T 1 of up to 8 h in this temperature regime were observed using a cavity QED protocol. Since no higher order processes were observed, the measured data are suitable to obtain a fundamental understanding of the spin-phonon coupling mechanism in this system. This paper is organized as follows: In Sec. II we derive an effective spin-phonon interaction V s−ph starting from the dipolar spin-spin interaction between two electrons and we give an expression for the spin-lattice relaxation rate Γ 1 . In Sec. III we explain the computational methods used to calculate Γ 1 ab initio by modelling both the electronic and phononic properties of the center by means of density functional theory. In Sec. IV the influence of lattice defects on Γ 1 is investigated and a comparison of our results with experimental data is presented followed by the conclusion in Sec. V.

II. THEORY
The idea to couple the spins to the phonons starting from the dipolar spin-spin interaction goes back to Waller 20 and was the first impact on spin-lattice relaxation in general. The relaxation mechanism is depicted in Fig. 1 1: (a) In the static case (no phonons in the system) the spins of the NV − center interact via the static spin-spin interaction, which is responsible for the fine structure splitting of the 3 A 2 ground state. (b) In the case the electron positions are coupled to phonons (right side of (b)), H ss becomes phonon dependent and can induce a spin flip. This allows to exchange energy between the spin system and the lattice and therefore allows the spin-system to equilibrate with the phonon bath.
phonon is excited, the dipolar spin-spin interaction between the i-th and j-th electron is altered, because the electronic distance vector r ij depends on the displacements of the ions {Q m }. Here µ 0 denotes the vacuum permeability, g e the g-factor 21 of the electron which is close to that of a free electron in the NV − center with a value of 2.0028, µ B is the Bohr magneton and S i and S j are the spin vectors of the i-th and j-th electron. In his original work Waller neglected the orbital character of the electrons and treated them as point sources located at the positions of ions. This assumption will be dropped in the following derivation. The change of the position of the electron with the ionic motion has to be modelled to couple the electronic spin vectors in H ss to the ionic movements. This is achieved by defining a region Ω around each ion, in which the electronic orbital follows the movement of the ions rigidly, resulting in the electronic distance vector For our calculations we use the Wigner-Seitz cell for Ω, dividing space geometrically. Since the ionic displacements in the low temperature regime are very small (mean square displacements are in the order of 10 −4Å −2 ), a Taylor expansion to first order in Q m is sufficient to calculate the transition rates between a m s = ±1 and a m s = 0 state. Thus, the spin-phonon interaction reads: To extract the relevant matrix elements responsible for a transition between the 3 A 2 levels, the spin operators S are expanded in raising and lowering operators. The Hamiltonian in (3) contains terms (a · S i )(b · S j ), where a and b are elements of {r ij , Q m }, and a term S i S j , which can be rewritten as single spin flip events and The only matrix elements, which can cause a transition in the ground state triplet, are those, that contain only a single raising or lowering operator and are underbraced in Eq. (4), the remaining terms account for double spin-flip or no spin-flip events. Taking only the spin-flip matrix elements of V s−ph in Eq. 3 into account we obtain the spin-flip potential Likewise the ionic displacements Q m are written in second quantized form 22 Here M m denotes the mass of the ion, N the number unit cells, a † /a the raising/lowering operator, m q,ρ the polarization vector of the m−th ion in the mode and R 0 m is the equilibrium position of the ion. Substituting Eq. (7) into Eq. (3) and taking V s−ph as a time dependent perturbation of the system leads to a Fermi's golden rule expression for a transition between the m s = ±1 and the m s = 0 states. The overall transition rate Γ f ←i from an initial to a final electronic state is obtained by a summation of the matrix elements of all final phonon states obeying energy conservation The raising and lowering operators acting on the initial phononic state |N i give N ph + 1 and N ph as eigenvalues of the particular state, where N ph is the occupation number of the phonons. We assume N ph to be the thermal occupations following the Bose-Einstein distribution. Putting everything together and considering the fact that only phonons with a single frequency at the spin-transition energy D can take part in this process, the transition rates for emission and absorption of a phonon Γ f ←i read To emphasize the temperature dependence of the relaxation rate this is rewritten as Γ f ←i = (N ph + 1)Γ 0 for emission of a phonon N ph Γ 0 for absorption of a phonon (10) with Γ 0 being the transition rate at zero temperature. To simulate an ensemble of spins relaxing from a nonequilibrium spin-distribution to equilibrium with the environment, both deexcitations and excitations of spins have to be considered 23 and the following rate equations have to be solved for our system with a degenerate excited state: The solution is straightforward by introducing the occupation difference ∆N = N ms=±1 − N ms=0 , since it obeys a simple exponential decay law to its thermal equilibrium value according to The calculation of the zero temperature transition rate Γ 0 between the 3 A 2 sublevels is sufficient to extract the transition rates in the low temperature regime, where single phonon processes are dominating over two phononprocesses.

III. METHODS
The calculation of Γ 0 requires the spin-polarized electronic orbitals as well as the phononic bandstructure, density of states and the polarization vectors for all the modes. We perform ab initio calculations using density functional theory on a supercell containing 64 lattice sites with one NV − center employing the Vienna Ab initio Simulation Package (VASP 24 ) using projector augmented wave pseudopotentials. 25 We use the local density approximation and a generalized gradient approximation using the PBE 26 exchange correlation potential for structural relaxations and the HSE functional 27,28 on top to obtain a proper bandstructure. Plane waves up to a cutoff of 800 eV are included and the first Brillouin zone is sampled with a 5 × 5 × 5 Monkhorst Pack grid. 29 A subtle relaxation of the ions resulting in forces on the atoms of less than 1 meV A shows that the neighbouring carbon atoms and the nitrogen atom move away from the vacancy, where the nitrogen atom is further displaced in accordance with an earlier study. 30 Since we are interested in the spin-polarized orbitals, we use the relaxed positions to calculate the electronic band structure. It is found that the a 1 ,e x and e y orbitals are located inside the bandgap and that e x and e y are the spin-polarized orbitals (see Fig. 2a). This familiar result 21,30-32 allows to extract these orbitals by applying the wannier90 package 33,34 to obtain the maximally localized orbitals on the nearest neighbour atoms of the vacancy. By considering the symmetry of the defect we add up the maximally localized orbitals to fulfill the C 3v symmetry constraints and end up with the spin polarized e x and e y orbitals of the 3 A 2 groundstate shown in Fig. 2b. Building Slater determinants with these orbitals we cal- |m i s , which occur for every phononic polarization in Eq. (10). The phonons are modelled by using the small displacement method within the harmonic approximation similar to a previous study. 35 We use the PHONOPY package 36 to extract the necessary displacements to build the dynamical matrix and apply it to the diagonalization thereof. We sample the Brillouin zone with a very dense mesh to extract 10 000 phonon polarization vectors at the transition frequency per band and the re-

IV. RESULTS
After carrying out the calculations, we end up with a theoretically predicted temperature dependent relaxation rate Γ, which can be compared to the experiment. As shown in Fig. 3a we find a direct single phonon process 23,37 at temperatures above the spin-transition T D/k B where thermal phonons excite and deexcite the spins by induced emission or absorption resulting in a linear dependence of Γ on T . This temperature dependence stems from the high temperature limit of the Bose-Einstein distribution, where N ph ∝ T in Eq. (12). At temperatures below the spin-transition the 2.88 GHz phonons start to freeze out and the only decay channel left for a spin-transition is the temperature independent spontaneous emission of a phonon occuring with a rate Γ 0 which results in the observed plateau in the low temperature regime.To compare the calculated rates with experiment, the treatment of the samples has to be explained: To create NV − centers in diamond, samples with a high initial nitrogen concentration (type Ib diamond) are irradiated by electrons, neutrons or ions in order to obtain vacancies followed by an annealing procedure. [38][39][40] The influence of the radiation damage on the phononic density of states is essentially unknown, but irradiation will create point defects, which can shift the phononic density of states towards lower energy excitations. 41 We simulate this effect and calculate the density of states for diamond crystals with defects and compare them with a perfect crystal. Introducing point defects (substitutional nitrogens and vacancies) in the diamond structure the phononic density of states indeed shifts towards lower frequencies as illustrated in Fig. 3b . However, we can only estimate the phononic density of states in the irradiated crystals. We model the phonons with the phononic densities of states for the simulated cells and the calculated relaxation rates Γ 0,ab initio = 2.1 × 10 −5 s −1 to 3.4 × 10 −5 s −1 are close to the lowest experimental values (Γ 0,exp = 3.47(16) × 10 −5 s −1 ). 19

V. CONCLUSION
In this paper we have shown that the very low spinlattice relaxation rates of the NV − center in diamond can be explained by the change of the dipolar spin-spin interaction induced by the movement of the ions as proposed originally by Waller in 1932. 20 We coupled the electronic distance vector r ij to the ionic movement by a first order Taylor expansion in the ionic displacement vectors {Q m } and ended up with an effective spin-phonon interaction V s−ph . We applied this interaction as a perturbation of the system to calculate the transition rates between the 3 A 2 ground state spin triplett ab initio using density functional theory by modelling the electronic wavefunctions and the phonons in a supercell containing 64 atoms. The calculated relaxation rates are comparable to the ones measured for samples that show little crystal damage. We propose that the deviation to samples with a strong irradiation damage is caused by the difference in the phononic density of states due to the irradiation treatment. Knowing the fundamental mechanism of spin-phonon interaction in this system will allow us to further investigate on higher order two phonon Raman processes and Orbach processes at higher temperatures for the spins of the NV − center in diamond.