Properties of nitrogen-vacancy centers in diamond: group theoretic approach

We present a procedure that makes use of group theory to analyze and predict the main properties of the negatively charged nitrogen-vacancy (NV) center in diamond. We focus on the relatively low temperatures limit where both the spin-spin and spin-orbit effects are important to consider. We demonstrate that group theory may be used to clarify several aspects of the NV structure, such as ordering of the singlets in the ($e^2$) electronic configuration, the spin-spin and the spin-orbit interactions in the ($ae$) electronic configuration. We also discuss how the optical selection rules and the response of the center to electric field can be used for spin-photon entanglement schemes. Our general formalism is applicable to a broad class of local defects in solids. The present results have important implications for applications in quantum information science and nanomagnetometry.


Introduction
Nitrogen-vacancy (NV) centers have emerged as promising candidates for a number of applications [1]- [4], ranging from high spatial resolution imaging [5] to quantum computation [6]. At low temperatures, the optical transitions of the NV center become very narrow and can be coherently manipulated, allowing for spin-photon entanglement generation [7] for quantum communication and all optical control [8]. A detailed understanding of the properties of this defect is critical for many of these applications. Although several studies have addressed this issue both experimentally and theoretically [9]- [20], not all aspects of the NV center are yet understood in detail. Furthermore, other atom-like defects can potentially be engineered in diamond [21] and other materials with similar, or perhaps better, properties for the desired application. Therefore, it is of immediate importance to develop a formalism to analyze and predict the main properties of defects in solids.
Here, we present a formalism based on standard group theory combined with high-level ab initio calculations. We apply group theory to identify the non-zero interactions between the given states, while the strength of the interaction will be calculated by using the wave functions and self-consistent potentials from atomistic ab initio calculations. We show that this combination is very powerful for explaining, or even predicting, several properties of the NV center. This treatment can be particularly useful in lowering the computational cost of point defects of such complexity that may be treated at the ab initio level, making a very detailed study of point defects possible.
While we focus on describing the NV center in diamond, the formalism can be applied to any point defect in solid state physics. The method takes advantage of the symmetry of the states to properly treat the relevant interactions and their symmetries. We apply group theory to find out not only the symmetry of the states but also their explicit form in terms of orbital 3 and spin degrees of freedom. We show that this is essential to build an accurate model of the NV center. In particular, we analyze the effect of the Coulomb interaction and predict that the ordering of the triplet and singlet states in the ground state configuration is { 3 A 2 , 1 E, 1 A 1 } and that the distance between them is of the order of the exchange term of the electron-electron Coulomb energy. This ordering has been discussed over the last few years [14,16,17], and our results agree with earlier [12] and recent ab initio calculations performed in bulk diamond [22].
The method is also used to analyze important properties of the center, such as polarization selection rules. The explicit form of the states allows us to identify a particularly useful lambdatype transition that was recently used for spin-photon entanglement generation [7]. We consider perturbations that lower the symmetry of a point defect, such as strain and electric field, and how they affect the polarization properties. We also show that the non-axial spin-orbit interaction discussed in [15] does not mix the basis states of the center in a given multiplet. Instead, we find that the electron spin-spin interaction is responsible for the spin state mixing of the excited state as a result of the lack of inversion symmetry of the center, contrary to what has been suggested in the literature [14,15,18]. Finally, we analyze the effect of electric fields via the inverse piezoelectric effect and compare it to experimental observations. We show that this effect can be used to tune the polarization properties of optical transitions and the wavelength of emitted photons, which is of direct importance for photon-based quantum communication between NV centers. The present study clarifies important properties of NV centers and provides the foundation for coherent interaction between electronic spins and photons in solid state.
This paper is organized as follows. In section 2, we present a general group theoretical formalism to calculate the electronic or hole representation of a point defect for a given crystal field symmetry and number of electrons contained in the defect. Using group theory and the explicit form of the states, we analyze the effect of the Coulomb interaction between electrons (section 3) and spin-spin and spin-orbit interactions for the NV center (sections 5 and 4, respectively). Next, the selection rules of the unperturbed defect are analyzed in section 6. Finally, in section 7, we analyze the effect of strain and electric field perturbations.

State representation
We are particularly interested in quasi-static properties of defects in crystals where the complex electronic structure can be observed spectroscopically. In this limit, one can apply the Born-Oppenheimer approximation to separate the many-body system of electrons and nuclei. This approximation relies on the fact that nuclei are much slower than electrons. In this approximation, the nuclei are represented by their coordinates and the physical quantities of the electrons depend on these coordinates as (external) fixed parameters. A defect in a crystal breaks down the translational symmetry, reducing the symmetry of the crystal to rotations and reflections. These symmetries form a point group that in general is a subgroup of the point group of the lattice. The loss of translational symmetry indicates that the Bloch states are no longer a good approximation to describe the point defect. In fact, some states can be very well localized near the point defect. These defect states are particularly important in semiconductors and insulators when they appear within the fundamental band gap of the crystal.
In the tight-binding picture, the electron system of the diamond crystal may be described as the sum of covalent-type interactions between the valence electrons of two nearest neighbor atoms. When defects involve vacancies, the absence of an ion will break bonds in the crystal, producing unpaired electrons or dangling bonds, σ i , which to leading order can be used to 4 represent the single electron orbitals around the defect [10]. The particular combination of dangling bonds that form the single electron orbitals {ϕ r } is set by the crystal field of the defect and can be readily calculated by projecting the dangling bonds on each irreducible representation (IR) of the point group of the defect [23], where P (r ) is the projective operator to the IR r , χ (r ) e is the character of operation R e (element) for the IR r , l r is the dimension of the IR r , and h is the order of the group (number of elements). Applying equation (1) to the NV center leads to the following single electron orbital basis (see appendix A): [10,12]. The non-degenerate orbitals a 1 (1) and a 2 (1) are totally symmetric and transform according to the one-dimensional (1D) IR A 1 ; meanwhile, {e x , e y } are two degenerate states that transform according to the 2D IR E. At this stage, group theory does not predict the energy order of these states. However, a simple model of the electron-ion Coulomb interaction can be used to qualitatively obtain the ordering of the levels [24]. In appendix A, we model the effect of this interaction on the single electron orbitals, ϕ r , for the case of the NV center and find that the ordering of the states (increasing in energy) is a 1 (1), a 1 (2) and {e x , e y }. Indeed, ab initio density functional theory (DFT) calculations revealed [12,25] that the a 1 (1) and a 1 (2) levels are lower than the e x and e y levels, which demonstrates the strength of the combination of group theory and simple models of interactions for qualitative predictions.
The dynamics of the center are set by the number of electrons available to occupy the orbitals. In the case of the negatively charged NV center, each carbon atom contributes one electron, the nitrogen (as a donor in diamond) contributes two electrons, and an extra electron comes from the environment [25,26], possibly given by substitutional nitrogens [27]. The ground state configuration consists of four electrons occupying the totally symmetric states and the remaining two electrons pairing up in the {e x , e y } orbitals. In this single particle picture, the excited state configuration can be approximated as one electron being promoted from the a 1 (2) orbital to the e x,y orbitals [12].
If two more electrons were added to any of these configurations, the wavefunction of the defect would be a singlet with a totally symmetric spatial wavefunction, equivalent to the state of an atom with a filled shell [28,29]. Therefore, the electronic configuration of this defect can be modeled by two holes occupying the orbitals e x,y in the ground state (e 2 electronic configuration) and one hole each in the orbitals a 1 (2) and e x,y for the excited state (ae electronic configuration). A third electronic configuration, a 2 , can be envisioned by promoting the remaining electron from the orbital a 1 (2) to the orbitals e x,y . Hole and electron representations are totally equivalent and it is convenient to choose the representation containing the smallest number of particles. If a hole representation is chosen, some care must be taken, as some interactions reverse their sign, such as the spin-orbit interaction [28]. In what follows, we choose a hole representation containing two particles (instead of an electron representation containing four particles), since it is more convenient to describe the physics of the NV center. However, the analysis can be applied to electrons as well.
The representation of the total n-electron wavefunction, including space and spin degrees of freedom, is given by the direct product of the representation of each hole hn and its spin = n ( hn ⊗ D 1/2 ), where D 1/2 is the representation for a spin-1 2 particle in the corresponding point group. The reduction or block diagonalization of the representation Table 1. Partner functions of each IR for the direct product of two holes. The first column shows the electronic configuration and in parentheses their triplet (T) or singlet (S) character. The third column gives the symmetry of the row of the IR. α (β) stands for ↑ (↓) and E ± = |a 1 e ± − e ± a 1 , where e ± = ∓(e x ± ie y ), |X = (|E − − |E + )/2 and |Y = (|E − + |E + )i/2.

Configuration
State Symmetry gives the basis states of the Hamiltonian associated with the crystal field potential and any interaction that remains invariant under the elements of the point group in question. These interactions include spin-orbit, spin-spin and Coulomb interactions, as well as expansions, contractions and stress where their axes coincide with the symmetry axis of the defect. The states can be found by projecting any combination of the two-electron wavefunction onto the IRs of the group [23,30], where ϕ i can be any of the orbitals in equation (1), χ i represents the spin wavefunction, and the subindex i refers to the hole i. U e is the corresponding operator of element e in the SU (2) representation. In the case of the NV center, it is illustrative to note that the spin representation for the two particles can be reduced to D 1/2 ⊗ D 1/2 = A 1 + A 2 + E, where A 1 corresponds to the singlet state, and A 2 and E to the triplet state with zero and nonzero spin projections, respectively. A list of the states and their symmetries for the two-hole representation can be found in table 1 for the ground state (e 2 ) and the excited state (ae). For completeness, we include the doubly excited state (a 2 ) electronic configuration, although this state is not optically accessible in the excitation process of the NV center in experiments. Note that each electronic configuration might have singlet and triplet states. We have written the degenerate ground state triplet as |αα and |ββ instead of the natural basis |αα ± |ββ as they immediately split into |αα and |ββ , even for small magnetic fields. We have used the complex notation E ± to identify 6 the orbital parts that factorize the spin parts on the excited triplet states. The orbital part will be important in the discussion of the optical transitions in the limit of low strain in section 6. For the degenerate states E x,y (ae), we write the orbital part as X, Y as they split (or get mixed), even for low strain. Note that each electronic configuration might have singlet and triplet states. Group theory can explain why the hyperfine interaction with the nuclear spin of the nitrogen in the excited state is more than an order of magnitude larger than in the ground state for both nitrogen species: the non-zero spin density in the ground state wavefunction of the NV center is mostly concentrated in the orbitals e x,y , which have no overlap with the nitrogen atom. On the other hand, in the excited state, when one electron is promoted from the a 1 (2) orbital to one of the e x,y orbitals, the non-zero spin density comes now from unpaired electrons occupying the orbitals a 1 (2) and e x,y . As the orbital a 1 (2) is partially localized on the nitrogen atom, a sizable contact term interaction between the electronic spin and the nuclear spin of the nitrogen is expected [25,31,32].
Until now, states inside a given electronic configuration have the same energy, but the inclusion of the electron-electron Coulomb interaction will lift the degeneracy between triplets and singlets. The resulting energy splitting can be of the order of a fraction of an eV, and it is analyzed for the ground state configuration of the NV center in section 3. Furthermore, the degeneracy of triplet states is lifted by spin-orbit and spin-spin interactions of the order of GHz, where the crystal field plays an important role. These interactions will be treated in sections 4 and 5.

Ordering of singlet states
For a given electronic configuration, the most relevant interaction is the electron-electron Coulomb interaction, which is minimized when electrons are configured in an antisymmetric spatial configuration. As the total wavefunction must be antisymmetric for fermionic particles, the spin configuration must be symmetric. As a result, the state with the largest multiplicity lies lower in energy. This analysis, known as the first Hund's rule, predicts that the ground state of the NV center should be the triplet 3 A 2 state [11,12]. We now address the question related to the order of singlets in the ground state electronic configuration e 2 . The order of singlet states has great significance in understanding the spin-flipping fluorescence of the NV center and ab initio DFT calculations were unable to address this issue properly due to the many-body singlet states. Since we have the explicit form of the wavefunctions, we can work out the ordering of the singlets in a given electronic configuration by analyzing the expectation value of the Coulomb interaction, which can be written in the general form Using this expression, we find that in the ground state electronic configuration (e 2 ), the Coulomb interactions for these states are 7 where x, y correspond to e x , e y states. From this set of equations, we find that the spacing between the singlets 1 A 1 and 1 E 2 is equal to the spacing between the singlet 1 E 1 and the ground state 3 where the difference is the exchange Coulomb interaction, which is always positive [33]. In addition, as 1 E 1 and 1 E 2 belong to the same IR E, it can be shown that . Under this consideration, the ordering of the states is It should be noted that, in this case, the most symmetric state has higher energy since the Coulomb interaction between two electrons is repulsive. This picture might be modified by the following effect. Since the Coulomb interaction transforms as the totally symmetric IR, the matrix elements between states with the same symmetry are non-zero. The states 1 E(e 2 ) and 1 E(ae) can couple via the Coulomb interaction, increasing the gap between them. A similar effect happens with the states 1 A 1 (e 2 ) and 1 A 1 (a 2 ). In equation (3), we did not take into account the effect of the other electrons present in the system. Nevertheless, our basic results here serve as a qualitative estimate for the energy of levels and provide useful insights into the structure of the NV center. The results of very recent calculations based on many-body perturbation theory (MBPT) [22] support our conclusion.

Spin-orbit interaction
The spin-orbit interaction lifts the degeneracy of multiplets that have non-zero angular momentum, and is also responsible for transitions between terms with different spin states [28]. It is a relativistic effect due to the relative motion between electrons and nuclei. In the reference frame of the electron, the nuclear potential, φ, produces a magnetic field equal to ∇φ × v/c 2 . In SI units, this interaction is given by [23] where V = eφ is the nuclear potential energy, m e is the electron mass and p k (s k ) is the momentum (spin) of electron k. The presence of the crystal field breaks the rotational symmetry of this interaction. Since φ is produced by the nuclear potential, it transforms as the totally symmetric representation A 1 , and therefore ∇V = (V x , V y , V z ) transforms as a vector, where Since p also transforms as a vector, it is possible to identify the IRs to which the orbital operator components belong. In C 3v , the components of ∇V and p transform as (E 1 , E 2 , A 1 ) and therefore O transforms as the IRs  [11], other matrix elements, such as e x |O x |e y , might be different from zero as they contain the IR A 1 . However, it is clear that they are zero by considering the Hermitian character of the interaction. In this case, the spin-orbit interaction can be written in terms of the angular momentum operators l i and takes the following form, where λ x,y (λ z ) denotes the non-axial (axial) strength of the interaction. In a system with T d or spherical symmetry, A = B (therefore λ x,y = λ z ) and the usual form (S · L) of the spin-orbit interaction is recovered. It is also useful to think about e ± as p ± orbitals and a 1 (2) as a p z orbital, where the angular momentum operators satisfy l ± a 1 (2) ∝ e ± [34].
Once it is known how the spin-orbit interaction acts on the orbitals, e x , e y and a, it is possible to calculate the effect of this interaction on the 15 states given in table 1. An important effect is the splitting in the excited state triplet between the states A 1 , A 2 and E x , E y and between states E x , E y and E 1 , E 2 [11]. The spin-orbit interaction can be written as Another effect, relevant when treating non-radiative transitions, is that the axial part of the spin-orbit interaction (λ z ) links states with m s = 0 spin projections among states of the same electronic configuration, while the non-axial part (λ x,y ) links states with non-zero spin projections with singlets among different electronic configurations. In figure 1, we show the states linked by the axial and the non-axial parts of the spin-orbit interaction, for which non-radiative transitions might occur. In addition to the well-known transition between A 1 (ae) → 1 A 1 (e 2 ), we find that this interaction might also link E 1,2 (ae) → 1 E 1,2 (e 2 ) and in particular E x,y (ae) → 1 E x,y (ae). The latter transition may play an important role, as recent ab initio calculations have shown that the singlets 1 E x,y (ae) might lie very close in energy to the excited state triplet [22]. In our model, the non-axial part of the spin-orbit interaction, λ x,y (l + s − + l − s + ), does not mix the states of the excited state triplet with different spin projections because the raising and lower operators, l − and l + , link states of different electronic configurations. In principle, this interaction can mix the states of the excited state triplet to higher order; however, this mixing will be largely suppressed by the large gap that separates different electronic configurations. Likewise, the transverse part of spin-orbit cannot mix the singlet 1 E(e 2 ) with the ±1 spin projections of the ground state 3 A 2 . In this sense, there is no conflict with the singlet 1 E(e 2 ) state to lie lower in energy than the singlet 1 A 1 , as suggested in [14], because the inter-system crossing, presumably allowed by spin-orbit [28], will be weak and will not preferably populate the ±1 spin states of the ground state. We have numerically evaluated the ratio between the axial part and the transverse part of spin-orbit, λ z /λ x y = B/A = 0.75 using the functions e x and e y and a 1 (2) from ab initio calculations (see appendix D). This suggests that if the axial part of spin-orbit is 5.5 GHz [17], the non-axial part should be of the order of λ x y = 7.3 GHz and only couples singlet with triplet states, as shown in figure 1. Although the transverse spin-orbit is much larger than the value it was incorrectly assigned in [14] (0.2 GHz), we have shown that it does not cause mixing of different spin projections in the excited state triplet. The transverse spin-orbit might be important to consider in inter-system crossing transitions [28]. Figure 1. Energy diagram of the unperturbed NV center in diamond. Note that each electronic configuration can contain triplets (left column) as well as singlets (right column), which have been drawn in separate columns for clarity. Red arrows indicate allowed optical transitions via electric dipole moment interactions. The circular arrows between the states E 1,2 and E x,y represent the mixing due to spin-spin interaction (see figure 2). Dashed lines indicate possible non-radiate processes assisted by spin-orbit interaction. In the ground state (e 2 configuration), the distance between singlets and triplets is equal to the exchange energy of Coulomb interaction (2J ). The horizontal dashed blue line represents the orbital energy of the ground state (without including spin-spin interaction).

Spin-spin interaction
The spin-spin interaction between electrons is usually not present in systems with spherical symmetry, due to the traceless character of the magnetic dipole-dipole interaction. However, if the electron wavefunction is not spherically distributed, this interaction does not average out. Here we describe its effect on the excited state triplet of the NV center and we provide a numerical estimation of its strength. The spin-spin interaction can be written (in SI units) as where s i = 1 2 σ x , σ y , σ z are the spin operators of particle i and σ j ( j = x, y, z) are the Pauli matrices, β is the Bohr magneton, g is the Landé factor for the electron and µ 0 is the magnetic permeability of free space 8 . To analyze the effect of this interaction on the defect, it is useful to write the spatial and spin parts separately in terms of the IRs of the point group [11,23]. Then, it is straightforward to express this interaction in terms of the basis states of the defect (see appendix B), where the gaps between the m s = ±1 and m s = 0 projections and between A 1 and A 2 states are given by These parameters have been analyzed in [11]. However, the mixing term given by [19] = µ 0 4π was not considered in [11]. Figure 2 shows the effect of spin-orbit and spin-spin interactions on the excited state manifold. In particular, we find that the state A 2 has higher energy than the state A 1 (2 > 0), contrary to previous estimations [11] 9 . In addition, we find that the spin-spin interaction mixes states with different spin projections. This effect is the result of a lack of inversion symmetry of the NV center and it is not present in systems with inversion symmetry, such as free atoms or substitutional atoms in cubic lattices. This does not contradict group theoretical estimates as the mixed states transform according to the same IR (e.g. the E 1 and E y states both transform according to the IR E 1 ; see table 1). The mixing term is responsible for the observed lambda transitions [8].
We estimated these parameters using a simplified model consisting of the dangling bonds given in figure A.1 (in the appendix) for the three carbons and the nitrogen atom around the vacancy. The dangling bonds are modeled by Gaussian orbitals that best fit the wavefunction obtained by an ab initio DFT supercell calculation (see appendix D). The distance between atoms is also taken from these simulations. To avoid numerical divergences when r = 0, we estimate equations (9)-(11) in reciprocal space following [35]. The values for the zero field splitting ( es = 3 ), the gap between states A 1 and A 2 (4 ) and the mixing term between states E 1,2 and E x,y ( ) are given in figure 2(b). As ab initio calculations cannot accurately estimate the nitrogen population p N = |β| 2 in the single orbital state a 1 (2) (see appendix A for a definition of parameter β), we have plotted in figure 2 Spin-orbit Spin-spin Spin-spin (mixing part) Figure 2. Splitting due to spin-orbit and spin-spin interactions in triplet ae. (a) The axial part of the spin-orbit interaction splits the states {A 1 , A 2 }, {E x , E y } and {E 1 , E 2 } by λ z . The spin-spin interaction splits states with different spin projections and also splits the A 1 and A 2 states. Our theory predicts the A 2 state at higher energy than the A 1 state and that the states (E 1,2 ) and (E x,y ) are mixed. As the state A 1 has an additional non-radiative decay channel, it is possible to confirm this finding by measuring the lifetime of the state. Note that the splitting between A 1 and A 2 is a direct consequence of spin-orbit mixing of the spatial and spin parts of the wavefunction. (b) Values for the zero field splitting (3 ), the gap between the states A 1 and A 2 (4 ) and the mixing term ( ) due to spin-spin interaction in the excited state as a function of the nitrogen population, p N , in the state a 1 (2). The shadowed areas indicate the possible values for these parameters when the distance between the vacancy and the three carbons is increased between 0 and 3%, and the distance between the vacancy and the nitrogen is decreased between 0 and 4% of their excited state configuration. The solid lines correspond to the maximum (minimum) distance between the carbons (nitrogen) and the vacancy.
variations in the relative distance between the three carbons, the nitrogen and the vacancy. The distance between the carbons and the vacancy is increased between 0 and 3%; meanwhile, the distance between the nitrogen and the vacancy is decreased between 0 and 4% relative to their excited state configuration (solid lines). This shows how the spin-spin interaction depends on the distance between the atoms. For a value of p N around 20%, our estimations are in good agreement with experimental observations [15,17,31].

Selection rules and spin-photon entanglement schemes
Transitions might be dipole allowed if the matrix element in the length representation contains the totally symmetric IR, ϕ f |ε · r|ϕ i ⊃ A 1 , where ε is the polarization of the electric field. The matrix elements a|x|e x and a|y|e y are non-zero, from which it is straightforward to calculate the selection rules among the 15 states given in table 1 for the unperturbed center. This is shown in table 3. In addition to the well-known triplet-triplet transition [36], transitions are also Table 3. Selection rules for optical transitions between the triplet excited state (ae) and the triplet ground state (e 2 ), the singlets (ae) and the singlets (e 2 ), and the singlet (a 2 ) and the singlets (ae). Linear polarizations are represented byx andŷ, while circular polarizations are represented byσ ± =x ± iŷ. As an example, a photon with σ + polarization is emitted when the electron decays from state A 2 (ae) to state 3  allowed between the singlets 1 E(ae) ↔ 1 A 1 (e 2 ), 1 E(ae) ↔ 1 E(e 2 ) and 1 A(a 2 ) ↔ 1 E(ae). Other non-zero matrix elements are e y |y|e x = e y |x|e y = e x |y|e y = − e x |x|e x , which might allow the transition between the singlets 1 A 1 (e 2 ) and 1 E(e 2 ). Recent experiments by Rogers et al identified an emission from singlet to singlet [18], attributing it to the 1 A 1 (e 2 ) → 1 E(e 2 ) transition [20] in the infrared energy range. The 1 A 1 (e 2 ) → 1 E(e 2 ) singlet-singlet transition might also be allowed by phonons or mixing between these states and singlets of different electronic configurations via Coulomb interaction, as discussed in section 3. However, we note that a recent sophisticated MBPT calculation [22] indicated that the infrared transition at about 1.1 eV belongs to the 1 E(ae) → 1 A 1 (e 2 ) transition. The polarization dependence is the same for both types of singlet-singlet transitions, and the association of the near-infrared emission of the NV center with the 1 E(ae) → 1 A 1 (e 2 ) transition does not contradict the nature of the emission in the original [18] and recent [37] measurements including the polarization dependence study.
Once the selection rules are known for the defect, it is possible to realize interesting applications, such as spin-photon entanglement generation [38]. In the case of the NV center, the system can be prepared in the A 2 (ae) state. Next, the electron can spontaneously decay to the ground state 3 A 2− by emitting a photon with σ + (right circular) polarization or to the state 3 A 2+ by emitting a σ − polarized photon (see figure 3). As a result, the spin of the electron is entangled with the polarization (spin) of the photon. The implementation of this scheme is sensitive to strain, which will be analyzed in section 7.

The effect of strain
Strain refers to the displacement u of the atomic positions when the crystal is stretched by x [34]. It is a dimensionless tensor expressing the fractional change under stretching, e i j = ∂δ R i ∂ x j , and it can be produced by stress (forces applied to the solid structure), electric field or temperature [39]. A systematic study of strain can be used to unravel the symmetry of defects and explore their properties [40]. Strain can shift the energy of the states as well as mix them. It can reduce the symmetry of the crystal field by displacing the atoms. However, not all nine components of strain change the defect in a noticeable way. The antisymmetric part of e i j transforms as a generator of the rotational group and therefore only rotates the whole structure. The symmetry and energies of the unperturbed states do not change upon rotation. Only the symmetric part of strain, = e + e T , affect the structure of a defect [34]. As with any other element of the theory, strain can be expressed in terms of matrices that transform according to the IRs of the point group under consideration. These matrices can be found by projecting a general strain matrix on each IR, In appendix C, we show in detail how to determine the effect of strain on the basis states of the defect. For simplicity, in the case of the NV center, we only write the effect of strain in the manifold {e x , e y , a}, where δ a A1 = (e x x + e yy )/2, δ b A1 = e zz , δ a E1 = (e x x − e yy )/2, δ a E2 = (e x y + e yx )/2, δ b E1 = (e x z + e zx )/2, δ b E2 = (e yz + e zy )/2 and in the manifold {e x , e y , a}. The effect of strain on the orbitals a, e x , e y is easy to see. A a 1 will shift equally the energies of the states e x and e y , while A b 1 will shift the energy of states a. Note that both describe axial stress: the former leaves the e 2 electronic configuration unaffected and the latter leaves the ea configuration unaffected. Either one produces relative shifts between both configurations, resulting in an inhomogeneous broadening of the optical transitions. However, they do not change the selection rules. Only the stress A a 1 + A b 1 , corresponding to either expansion or contraction, leaves all relative energies unaffected. E a 1 splits the energy between e x and e y , and E a 2 mixes the two states. Finally, E b 1 and E b 2 mixes the states e x and a and the states e y and a, respectively. In the case of the NV center, the effect of the matrices E b 1,2 can be neglected thanks to the large gap between the orbitals a and e x,y . Therefore, in what follows, we do not consider them further. Recently, studies have been performed to analyze how strain affects the excited state structure of the NV center [15,17,18]. Here we derive the explicit form of strain affecting the different electronic configurations and look at how strain affects the selection rules described in section 6. The relevant strain matrices that will lower the C 3v symmetry in each electronic configuration are E a 1 and E a 2 , for which the Hamiltonian is This mostly affects the singlet and excited state configurations in the following form, for the manifolds The ground state, due to its antisymmetric combination between e x and e y , is stable under the perturbation H strain . This can be checked by applying equation (15) to the ground state given in table 1. The effect on the excited state triplet can be seen in figure 4(a), where the unperturbed states are mixed in such a way that, in the limit of high strain, the excited triplet structure splits into two triplets with spatial wavefunctions E x and E y . When strain overcomes the spin-orbit interaction (δ a E1 > 5.5 GHz), the spin part decouples from the spatial part and the total angular momentum is no longer a good quantum number. Transitions from the excited state triplet to the ground state triplet are linearly polarized, where the polarization indicates the direction of strain in the x y-plane. Figure 4(c) shows how the polarization of the emitted photon from the state A 2 to the ground state 3 A 2− varies from circular to linear as a function of strain. In the case of δ a E2 strain, the effect is similar but now the mixing is different. As shown in figure 4, A 2 mixes with E 1 and the photons become polarized along x-y. Note that, in the limit of low strain, in both cases the polarization remains right circularly polarized for the transition between the excited state A 2 (ae) to the ground state 3 A 2− (e 2 ), while the polarization remains left circular for the transition between the excited state A 2 (ae) to the ground state 3 A 2+ (e 2 ). The fact that at lower strain the character of the polarization remains circular has been successfully used in entanglement schemes [7]. The polarization properties of the states E 1,2 are similar to those of the states A 1,2 but with the opposite polarization.

Strain and electric field
The application of an electric field to a defect leads to two main effects. The first effect, the electronic effect, consists of the polarization of the electron cloud of the defect, and the second one, the ionic effect, consists of the relative motion of the ions. It has been shown that the two effects are indistinguishable, as they have the same symmetry properties [41]. The ionic effect is related to the well-known piezoelectric effect. When a crystal is under stress, a net polarization P i = d i jk σ jk is induced inside the crystal, where d i jk is the third-rank piezoelectric tensor and σ jk represents the magnitude and direction of the applied force. Conversely, the application of an electric field might induce strain given by jk = d i jk E i , where E i are the components of the electric field [39]. The tensor d i jk transforms as the coordinates x i x j x k and, therefore, group theory can be used to establish relations between its components for a given point group. In particular, the non-zero components should transform as the IR A 1 . By projecting d i jk (or x i x j x k ) onto the IR A 1 , we can determine the non-zero free parameters of the tensor d and determine the effect of electric field on the basis states of the unperturbed defect (see appendix C). In the case of the NV center, the effect on the excited state triplet is given by the following matrix, In the presence of strain along theŷ direction (dashed lines), the response is quadratic due to the splitting between E x and E y states in the excited state. Our numerical results are in fair agreement with experimental results [43].
in the basis {A 1 , A 2 , E x , E y , E 1 , E 2 }, while the effect on the ground state triplet is in the basis The parameters a, b and d are the components of the piezoelectric tensor d i jk and g is the coupling between the strain tensor e and the NV center. Comparing equations (17) and (18), we note that the linear response of the excited state and ground state are, in principle, different. An electric field along theẑ (NV axis) can be used to tune the optical transition without distorting the C 3v symmetry of the defect, provided b = d. In figure 5, we show the linear response of NV centers under an electric field parallel to the NV axis. In this case, the linearity is not affected by the presence of strain. Our estimates for the ionic effect, based on the response of the lattice defect to electric field and the response of the orbital energies to strain (see appendix C), indicate that the relative shift between the ground and excited states is about 4 GHz MV −1 m −1 . This electric field control [15] could be very important in schemes to entangle two NV centers optically as the wavelength of the photons emitted from each NV center need to overlap [42]. In figure 5, we show the response of optical transitions under an electric field perpendicular to the NV axis. In this case, the response is linear if strain is absent and quadratic if strain is non-zero. Dashed lines show the response to an electric field when the defect experiences a 0.3 GHz strain along the [01-1] axis. Our estimates can be used to interpret the Stark shift observations of Tamarat et al [43].

Conclusions
We have used group theory combined with ab initio calculations to identify, analyze and predict the properties of NV centers in diamond. This analysis can be extended to other deep defects in solids. A careful analysis of the properties of a defect using group theory is essential for predicting spin-photon entanglement generation and for controlling the properties of NV centers in the presence of perturbations, such as undesired strain. We have shown that group theoretical approaches can be applied to determine the ordering of the singlets in the (e 2 ) electronic configuration and to understand the effect of spin-orbit, spin-spin and strain interactions.  }, which transform as the IR A 1 , the rotation operator R z as A 2 , and the pair of functions {(x, y), (R x , R y ), (x y, x 2 − y 2 ), (yz, x z)} as E. The spin projections {α(↑), β(↓)} transform as the IR E 1/2 (or D 1/2 ), while the functions ααα + iβββ and ααα − iβββ transform as the IRs 1 E 3/2 and 2 E 3/2 , respectively.
Next, we model the electron-ion interaction to find out the ordering of these states. This interaction can be written in the basis of the dangling bonds σ i as where v i < 0 is the Coulomb interaction of orbital σ i at site i, h c is the expectation value of the interaction between orbitals σ i and σ i+1 at site i = {1, 2, 3}, v n = σ N |V |σ N and h n = σ i |V |σ N .
This interaction, which transforms as the totally symmetric IR A 1 , not only sets the order of the orbitals but also mixes orbitals a N and a C . This is a consequence of the important concept that whenever a matrix element contains the totally symmetric representation, its expectation value might be different from zero [23]. Since both wave functions as well as the interaction between them transform as the totally symmetric representation A 1 , the representation for the matrix element also transforms as A 1 : This interaction leads to the new basis [10,12,25] (1) . We see that the most symmetric state is lowest in energy, which is usually the case for attractive interactions.

Appendix B. Spin-spin interaction
To analyze the effect of spin-spin interactions (equation (7)) from the perspective of group theory, we first rewrite this interaction to identify spatial and spin terms that transform as IR objects in the point group, wherex,ŷ andẑ are directional cosines and s ± = s x ± is y . In the case of C 3v , for the unperturbed center, the expectation values of the fourth and fifth terms are non-zero in the spatial manifold of the excited state {|X , |Y } because the center lacks inversion symmetry. However, these terms might be neglected when considering other defects with inversion symmetry. We note now that the spatial part of the first term transforms as the totally symmetric representation A 1 , while the second and third terms transform as the IR E. The reader can check which IR these combinations belong to by looking at the character table in appendix A. Therefore, their expectation values can be written as where |X and |Y are the two-electron states given in table 1. Note that, for symmetry reasons, the second and third relations are characterized by the same parameter , while the last two relations are characterized by the same parameter . Similarly, it is possible to write the spin operators in the spin basis of the two holes, {|αα , |αβ , |βα , |ββ }. For example, s 1+ s 2− = |αβ βα|. Using these relations and equation (B.1), the Hamiltonian in the fundamental bases of the excited state of the NV center is Finally, we can write H ss in terms of the basis states of the unperturbed defect (see table 1). This leads to equation (8).

Appendix C. Strain and electric field
The effect of strain on the electronic structure of the defect can be obtained from the effect of the electron-nuclei Coulomb interaction on the basis states of the defect. In our example, the Coulomb interaction is given by equation (A.1). However, when the positions of the atoms are such that the symmetry of the defect is reduced, we should allow for different expectation values of the matrix elements: h i j = σ i |V |σ j and h in = σ i |V |σ N . We have assumed that the self-interactions, v c and v n , do not change as the electrons follow the position of the ion according to the Born-Oppenheimer approximation. To relate the matrix elements to the ionic displacements, we can assume as a first approximation that the electron orbitals are spherical functions and therefore the matrix elements can be parameterized by the distance between ions, h i j (q i , q j ) = h i j (|q i − q j |), so that we can write The change in the matrix elements is linear in the atomic displacements. In turn, the atomic displacements are related to the strain tensor by δq i = eq i , and therefore the change in the matrix element is given by Under these considerations, it is straightforward to calculate the effect of strain on the basis states of the defect. For simplicity, we write here only the effect of strain on the degenerate orbitals, e x and e y , δV = −g e x x e x y e x y e yy , (C.3)

21
where g = 8q 3 ∂h i j ∂q i j and q is the nearest neighbor distance between atoms. Using the electron wavefunction obtained from ab initio calculations (see appendix D), we estimate that g ≈ 2 pHz (p = peta = 10 15 ).
The effect of electric field on the states of the defect can be analyzed by the inverse piezoelectric effect, as described in section 8. In this appendix, we show how group theory can determine the nature of the piezoelectric tensor. By projecting d i jk (or x i x j x k ) onto the IR A 1 , we can build the following relations, and the d tensor can be written in the following short notation (contracted matrix form) [39], For a given electric field, we have a strain tensor of the form To evaluate the magnitude of the piezoelectric response, we have used first-principles calculations, as described in appendix D. The values for the components of the piezoelectric tensor due to the ionic effect are a ≈ b ≈ c ≈ 0.3 µ(MV m −1 ) −1 and d ≈ 3 µ(MV m −1 ) −1 .

Appendix D. Information about the first-principles methods applied in our study
To determine the values of the constants a, b, c and d introduced in appendix C, we applied DFT [45] calculations within a generalized gradient approximation PBE (Perdew-Burke-Ernzerhof) [46]. In the study of spin-orbit and spin-spin interactions, we used a 512-atom supercell to model the negatively charged NV defect in diamond. In particular, we utilized the VASP code [47,48] to determine the geometry of the defect that uses the projector augmented wave method [49,50] to eliminate the core electrons, while a plane wave basis set is employed for expanding the valence wavefunctions. We applied the standard VASP projectors for the carbon and nitrogen atoms with a plane wave cutoff of 420 eV. The geometry optimization was stopped when the magnitude of the forces on the atoms was lower than 0.01 eV Å −1 . We calculated the geometry of both ground and excited states. We applied the constrained DFT method to calculate the charge density of the excited state, that is, by promoting one electron from the a 1 (2) orbital to the e x and e y orbitals, as explained in [32,51]. This procedure is a relatively good approximation, as confirmed by a recent MBPT study [22]. The obtained geometries from VASP calculations were used as starting points in the calculations of spin-orbit and spin-spin interactions. The spin-orbit energy was calculated by following equation (4) in this paper. Since the spin-orbit interaction is short range, we applied all-electron methods beyond the frozencore approximation. We utilized the CRYSTAL code [52] for this calculation using the PBE functional within DFT. We took the geometry as obtained from the VASP calculation. We applied the 6-31*G Gaussian basis set for both the carbon and nitrogen atoms. The calculated properties (such as the position of the defect levels in the gap) agreed well with those from planewave calculations. We obtained the all-electron single-particle states and the corresponding Kohn-Sham potentials on a grid and calculated the spin-orbit energy numerically.
Finally, we also studied the piezoelectric effect. In this case, an external electric field was applied along the NV axis and perpendicular to it. For this investigation, only a finite size model can be used; thus we modeled the negatively charged NV defect in a molecular cluster consisting of 70 carbon atoms and one nitrogen atom. The defect was placed in the middle of the cluster. The surface dangling bonds of the cluster were terminated by hydrogen atoms. In our previous studies, we showed [25] that the defect wavefunctions are strongly localized around the core of the defect; thus our cluster model can describe reasonably well the situation occurring in the bulk environment. For this investigation, we again applied DFT with the PBE functional as implemented in the SIESTA code [53]. We used the standard double-ζ polarized basis set and Troullier-Martins norm-conserving pseudo-potentials [54]. This method gives identical results to those obtained from plane wave calculations regarding the geometry and the wavefunctions in supercell models [25]. We fully optimized the defective nanodiamond with and without the applied electric field. In this case, we applied a very strict limit to the maximum magnitude of forces on the atoms, 0.005 eV Å −1 . We applied six different values of the external electric field along the NV axis and in perpendicular directions to it, where we could clearly detect the slope of the curvature of atomic displacements versus the applied electric field. The resulting values for the atom displacements in the presence of 1 MV m −1 electric field are of the order of a few 0.1 µÅ .