Nitrogen-vacancy (NV) complex in diamond is one of the most prominent solid state defects as the negatively charged NV defect (NV ) is a leading contender for quantum technologies. In quantum information processing applications, NV is photoexcited that often leads to photoionization to neutral NV defect, NV 0, and re-ionization back to NV should occur to control the S = 1 spin of NV . As a consequence, understanding the photophysics of NV 0 is crucial for controlling NV . Furthermore, recent studies have shown that the S = 1 / 2 electron spin of NV 0 can also be initialized and read out at certain conditions that turns single NV 0 a potential quantum bit. Quantum optics protocols rest on detailed knowledge on the electronic structure of the given system, which is obviously missing for NV 0 in diamond. In this study, we combine the group theory and density functional theory calculations toward exploring the nature of the ground and excited states of NV 0. We show that the effective three-electron system of NV 0 leads to high correlation effects that make this system very challenging for ab initio simulations.

Nitrogen-vacancy (NV) center1 in diamond has been receiving a continuous interest due to its broad applicability for quantum technologies.2–8 Most of these applications utilize its negative charge state (NV ) that belongs to the 637 nm1 ( 1.945 eV) zero-phonon-line (ZPL) photoluminescence center. The NV color center possesses S = 1 electronic spin that can be initialized, coherently driven, and read out, thus NV may act as an exemplary qubit operating at room temperature.9–12 An early theoretical work suggested13 that the neutral counterpart, NV 0, may also act as a qubit. This has been recently realized in experiments,14–17 where NV 0 possess S = 1 / 2 electronic spin in the ground state with the corresponding 575 nm ( 2.156 eV) ZPL emission (Fig. 1). However, in stark contrast to the extensively studied NV , various magneto-optical properties of NV 0 are vaguely studied and the underlying physics is left unknown.

FIG. 1.

Electronic structure of the neutral nitrogen-vacancy ( N V 0) center of diamond. (a) Defect levels inside the bandgap of diamond. (b) Single particle (Kohn–Sham) orbitals of defect levels.

FIG. 1.

Electronic structure of the neutral nitrogen-vacancy ( N V 0) center of diamond. (a) Defect levels inside the bandgap of diamond. (b) Single particle (Kohn–Sham) orbitals of defect levels.

Close modal

Both 1.9452,18–20 and 2.156 eV21,22 color centers are active in both absorption (ABS) and photoluminescence (PL) processes. NV defect’s 1.945 eV ZPL is accompanied19 by a broad phonon sideband both in PL and ABS exhibiting only minor asymmetry [cf., Figs. 2(c) and 2(d)]. The sideband mainly couples with a ω = 64 meV quasilocal vibration and the strength of coupling can be depicted by a S 3.5 Huang–Rhys (HR) factor. Consequently, the ZPL exhibits a Debye–Waller factor of DW = e S 3 %, thus only 3% of total intensity resides in the ZPL in agreement with previous ab initio simulations within23,24 and beyond the HR approach.25 

FIG. 2.

Comparison of theoretical and experimental phonon sidebands for the negative and neutral NV defects. (a) and (c) Absorption line shape of the 2.156 and 1.945-eV peaks of NV 0 and NV , respectively. We scaled the experimental and theoretical line shapes to exhibit the same area under their curves between the ZPL and 100 meV, 200 meV for NV 0, NV , respectively. (b) and (d) Luminescence line shapes for 2.156-eV NV 0 and 1.945-eV NV . We normalized the experimental (expt.) and theoretical line shapes in their full range. Expt. data are taken from different sources as follows: [a] Fig. 1(b) from Ref. 21, [b] Fig. 2(a) data from Ref. 29 at 1400 °C anneal, [c] Ref. 32, [d] Fig. 4 in Ref. 33, [e] Fig. 6(b) in Ref. 19, and [f] Fig. 3(b) in Ref. 20.

FIG. 2.

Comparison of theoretical and experimental phonon sidebands for the negative and neutral NV defects. (a) and (c) Absorption line shape of the 2.156 and 1.945-eV peaks of NV 0 and NV , respectively. We scaled the experimental and theoretical line shapes to exhibit the same area under their curves between the ZPL and 100 meV, 200 meV for NV 0, NV , respectively. (b) and (d) Luminescence line shapes for 2.156-eV NV 0 and 1.945-eV NV . We normalized the experimental (expt.) and theoretical line shapes in their full range. Expt. data are taken from different sources as follows: [a] Fig. 1(b) from Ref. 21, [b] Fig. 2(a) data from Ref. 29 at 1400 °C anneal, [c] Ref. 32, [d] Fig. 4 in Ref. 33, [e] Fig. 6(b) in Ref. 19, and [f] Fig. 3(b) in Ref. 20.

Close modal

Both 1.945 and 2.156-eV ZPLs are attributed to transitions between the same a 1 and e electronic orbitals [Fig. 1(b)] that would imply comparable optical features. Surprisingly, however, there is a significant asymmetry in NV 0 defect’s PL and ABS phonon sidebands as plotted in Figs. 2(a) and 2(b), respectively. S = 2.35 HR factor of the PL spectrum was reported for NV 0 defect26 and another study reported S = 1.9 ( 1 ) at 77 K (Ref. 15). ABS spectra exhibit much broader sideband quantified by an S = 3.9 HR factor.21 Interestingly, the PL signal is coupled with phonons with significantly smaller energy ω 42–45 meV15,27 than those of NV . The radiative lifetime of the excited state of NV 0 is reported as 13.2(2),28 19(2),29 17(1),30 and 20(1)31 ns in different studies.

The 2.156-eV ZPL of NV 0 occurs between the 2 E orbitally degenerate ground state to an orbitally non-degenerate 2 A 2 confirmed by uniaxial stress measurements.21,34 Furthermore, 2 E splits to two Kramers doublets ( E 1 / 2 and E 3 / 2) due to spin–orbit coupling. In this context, we note that NV defect’s spin–orbit splitting was already reported as λ = 5.3 GHz35 or λ = 5.33 ( 3 ) GHz36 by various experiments more than a decade ago, yet the λ for NV 0 has not been reported until recently. The exact spin–orbit splitting value37 is still ambiguous as λ = 4.48 ( 10 )14 and 9.8(8)16 GHz were reported by two different research groups. Additionally, we note that electron paramagnetic resonance (EPR) for NV 0 defect’s S = 1 / 2 spin has not yet been reported, to our best knowledge. The lack of EPR signal is attributed to linewidth broadening caused by the Jahn–Teller effect.38,39

On the other hand, under continuous green illumination at 532 nm, an S = 3 / 2 spin associated with the metastable 4 A 2 state of NV 0 was indeed observed in EPR38 where early ab initio calculations confirmed the model based on the favorable comparison of the observed and calculated hyperfine signatures between the S = 3 / 2 electron spin and I = 1 / 2 15N [ A = 35.7 ( 3 ) MHz and A = 23.8 ( 3 ) MHz in experiments] and 13C nuclear spins.13 The observed near free electron g-factors ( g = 2.0029 ( 2 ), g = 2.0035 ( 2 )) imply very small second-order spin–orbit interaction associated with this state. In this case, the C 3 v symmetric axial field causes a zero-field splitting between the E 1 / 2 and E 3 / 2 levels at D = 1685 ( 5 ) MHz as observed in the EPR spectrum.38 

The existence of the metastable 4 A 2 state has been confirmed by ab initio calculations with various levels of theory.13,40–42 However, the level position referenced to the 2 E ground state scatters between 0.4 and 0.86 eV (see Table I). According to a recent theoretical study with combining selection rules and ionization potentials of NV (Ref. 42), successive two-photon absorption of NV via the 3 E excited state could populate the m s = ± 1 / 2 state over the m s = ± 3 / 2 state of the 4 A 2 state, which would explain the occupation of the 4 A 2 metastable state and the EPR signal of this metastable state under green optical pumping that is enhanced by the population difference in their spin levels. This theory sets the position 4 A 2 level over the 2 E level at around 0.4 eV.

TABLE I.

Electronic excited states for NV0. All energies are in eV units. We note that we optimized the atomic positions of atoms in our present work (p.w.) with HSE06, PBE, LDA functionals. However, the atomic positions were not relaxed within configuration interaction (conf. int.)c and Hubbard modeling (Hubbard)d results.

Method |4A2⟩ |2A2⟩ |2E*⟩ |2A1⟩ 
expt. ∼0.4a 2.156 n.a. n.a. 
HSE06 0.48p.w. 2.36p.w. 2.33p.w. ∼2.8p.w. 
PBE 0.58p.w. 1.20p.w. 1.24p.w. n.a 
LDA 0.86b 1.22p.w. 1.26p.w. n.a 
conf. int. 0.68c 1.65c 2.04c 2.93c 
Hubbard 0.68d 2.64d 2.88d 3.29d 
Method |4A2⟩ |2A2⟩ |2E*⟩ |2A1⟩ 
expt. ∼0.4a 2.156 n.a. n.a. 
HSE06 0.48p.w. 2.36p.w. 2.33p.w. ∼2.8p.w. 
PBE 0.58p.w. 1.20p.w. 1.24p.w. n.a 
LDA 0.86b 1.22p.w. 1.26p.w. n.a 
conf. int. 0.68c 1.65c 2.04c 2.93c 
Hubbard 0.68d 2.64d 2.88d 3.29d 
a

Although directly not known, the difference between photoionization thresholds of |3A2⟩ → |2E⟩ and |3E⟩ → |4A2⟩ are 432 and 342 meV, respectively, which hints for ∼0.415 eV energy position for the |4A2⟩ multiplet. However, these photoionization thresholds were recorded at room temperature where phonon-assisted transition can also occur, thus ∼0.415 eV is still an estimate.

b

Ab initio taken data from Ref. 13.

c

Ab initio data by means of screened configuration interaction on top of LDA (local density approximation of DFT) from Ref. 41.

d

Ab initio data by means of generalized Hubbard Hamiltonian on top of B3LYP functional from Ref. 40 in a C71H85 diamond molecular cluster.

This example shows that it is crucial to know the lifetime of this metastable state so one can reliably set up quantum optics protocols for this potential qubit state. This is associated with the intersystem crossing (ISC) between 4 A 2 and 2 E as well as the photoionization cross section from 4 A 2 to the ionization bands. To our knowledge, this issue has not yet been resolved for NV 0 as well as the origin on the phonon sidebands of the absorption and luminescence spectra. In this study, we aim to understand the photophysics of NV 0 and the non-radiative decay routes that are all associated with the strong electron–phonon coupling of the system.

The paper is organized as follows. In Sec. II, we describe our ab initio methodology. Then, we depict the electronic structure of NV 0 and establish the nomenclature of our paper in Sec. III. Next, we discuss the interactions that govern the 2 E ground state and its optical excitation at 2.156-eV in Sec. IV. After understanding both spin–orbit and electron–phonon couplings of NV 0, we show and discuss our results on the ISC rates of 4 A 2 in Sec. V and the photoionization threshold energies and processes in Sec. VI. Finally, we will conclude our paper in Sec. VII.

We apply ab initio calculations with the use of simple cubic 512-atom supercells embedding NV defect within the framework of the spinpolarized density functional theory (DFT) as implemented in the vasp 5.4.1 code.43 The 512-atom model suffices to sample the Brillouin zone of the supercell at the Γ point. We converged the electronic structure by self-consistent cycles with 10 5 eV convergence threshold. The applied projector-augmentation-wave method (PAW)44,45 on the core electron orbitals are made possible to use a relatively low cutoff (370 eV) plane wave basis for electron wavefunctions. We converged the forces acting on ions below 10 3 eV/Å.

We applied the Heyd–Scuseria–Ernzerhof (HSE06) hybrid functional46,47 that reproduces the experimental bandgap and the charge transition levels of defects in group-IV semiconductors within 0.1 eV accuracy.48,49 We determined the excitation energies for excited states within 0.1 eV accuracy with the ΔSCF method50 that provides accurate ZPL energy and Stokes-shift for the optical excitation spectra of the triplets of the NV center23,24 and group-IV impurity-vacancy centers51,52 of diamond.

We used the HR theory23,53,54 to predict the phonon sidebands near the ZPL emission. We used the Perdew–Burke–Ernzerhof (PBE) semilocal functional55 to determine the vibration frequencies and eigenvectors as PBE is reported to accurately predict these features.23 However, we employed the computationally demanding HSE06 to determine the optical excitation energies by means of the ΔSCF method where the atomic positions upon excitation are relaxed.

We determined the spin–orbit coupling as implemented in the vasp code within the scalar-relativistic approximation.56 During the non-collinear calculations for spin–orbit coupling, we fixed the spin quantization axis along the trigonal symmetry axis of the defect. We determined the spin–orbit coupling parameters directly from matrix elements between single particle Kohn–Sham orbitals that reproduced the observed fine structure of defects within 20 % accuracy.24,51,52,57,58

Here, we define the basic nomenclature of the paper. We note that the orbitals and levels of NV from DFT calculations have been already published in several papers.40,59–66 Furthermore, the respective many-body states and spin–orbit coupling coefficients acting between them were also thoroughly analyzed.67,68 As we will show in this paper, the strong electron–phonon coupling should be explicitly considered for the many-body states of NV 0 that was lacking in previous studies. Therefore, instead of referring toward these papers, we explicitly reiterate the respective wavefunctions and interactions that we use in the entire paper. Therefore, we begin our discussion presenting electronic structure of NV 0 in Sec. III A. We describe the different forms of electronic correlations governing multiconfigurational character | 2 E , | 4 A 2 , | 2 A 2 , | 2 E multiplets of N V 0 in Secs. III BIII D. We summarize the multiplet energies within different levels of theory from the literature accompanied by our present results in Table I.

The NV defect introduces an a 1 and a double degenerate e level in the gap (see Fig. 1) that are occupied by three electrons in the relevant neutral charge state. The | 2 E ground state has ( a a e ) electron occupation. The | 2 E multiplet is split into two | 2 E 1 2 and | 2 E 3 2 Kramer’s doublets as
(1)
where we label them by their m j= ± 1 2 , ± 3 2 total angular momentum that includes both spin m s=( , ) and effective orbital momentum ( m l= ± 1) from e levels by introducing | e ± = 1 2 ( | e x ± i | e y ) complex combination of the e { x , y } real orbitals. Electronic excitation of NV 0 will promote an electron from the | a 1 orbital to an | e orbital, giving rise to the ( a e e ) electronic configuration. We construct these states as tensor products of NV with ( e e) configurations coupled with an ( a) hole. The latter can be labeled | 2 A 1 while the ( e e) gives rise the following | 3 A 2 | 1 E | 1 A 1 two particle wave-functions:
(2)
We then take the combinations with the additional | a 1 orbital, thus the optically excited states can be labeled
(3)
According to Hund’s first rule, | 4 A 2 with maximized spin is the lowermost excited state
(4)
The triplet configuration | 3 A 2 may couple with | 2 A 1 with opposite spins that results in
(5)
multiplets. Additionally, | 1 E singlet can couple with | a 1 that gives rise
(6)
When | 1 A 1 couples with | 2 a 1 , then it results in | 2 A 1 that can be written as
(7)
Finally, we show the high-energy | 2 E state for the sake of completeness when all three electrons occupy the e levels as ( e e e) electronic configuration,
(8)
However, these configurations can be excluded from the discussion as the corresponding levels lie higher than the ionization bands.
We recognize in our spinpolarized DFT calculations that the symmetry of | 2 E is spontaneously broken even if the geometry is constrained to C 3 v symmetry. That is, the orbitals follow the reduced C 1 h symmetry where only a single σ ^ mirror plane remains. Therefore, the | a 1 level is allowed to mix with | e x because they both belong to A representation while | e y cannot mix because it transform as the antisymmetric A representation. In order to identify electronic correlation in | 2 E , we decompose the broken symmetry single determinant obtained by DFT into symmetric orbitals respecting C 3 v symmetry as follows:
(9)
where u = + 0.998, v = 0.06, w = + 0.85, l = 0.53, thus ( u w ) 2 = 0.71, ( v w ) 2 = 0.002, ( u l ) 2 = 0.28, ( v l ) 2 = 0.001. Our results indicate that there is correlation in | 2 E . We note that the non-negligible | 2 A 2 contribution occurs in the broken-symmetry solution as an artifact. This term should be exactly zero with a multiconfigurational wavefunction method. Finally, the | 2 E contribution is only perturbative and can be again neglected. We additionally note that the | a 1 e x e y configuration is not a pure | 2 E (see Sec. III D). However, we neglect this issue for the sake of simplicity and assume that all | a 1 e x e y electronic character belongs to | 2 E . To sum up, we assume that the two orbital doublet multiplets are mixed by the expressions
(10)
where Γ runs over Γ = { ± 1 2 , ± 3 2 }, and we used the c 2 ( u w ) 2 = 71 % ( 85%), d 2 ( u l ) 2 = 28 % ( 13 %) parameters inherited from Eq. (9). We derived these parameters purely from DFT by means of HSE06 (PBE) functional. For comparison multiconfigurational methods report these parameters as c 2 = 90 % and d 2 = 8 % (see Table II in the supplementary materials of Ref. 41).

We quantify the effect of electronic correlation in energy units by removing the symmetry constraint on the wavefunctions and enforcing C 3 v symmetry only on the atomic positions. We report 12.5, 139, 891 meV energy gain by means of LDA, PBE, HSE06 functionals, respectively, where we enforce symmetrization by putting two electron halves on | e x , | e y orbitals resulting in the ( a 1 1 a 1 1 e x 0.5 e y 0.5 e x 0 e y 0 ) electron occupation. We note that we determined the vibration frequencies and modes within this smeared occupation by means of PBE functional to avoid broken symmetry wavefunctions as given in Eq. (9). However, this introduces an artificial self-interaction as the electrons on e x and e y repel each other due to exchange integrals present in hybrid functionals such as HSE06. In summary, we highlight that the electronic correlation between the two | 2 E , | 2 E multiplets is a significant effect that cannot be neglected. We set the reference point of our electronic total energy (for example, in Fig. 3) to this broken-symmetry solution on which we removed all possible symmetry constraints both on electronic and ionic degrees of freedom with the spinpolarized DFT method.

FIG. 3.

Jahn–Teller effect for the ground state of NV 0. (a) Adiabatic potential energy surface (APES) of NV 0 defect’s | 2 E ground state. We relaxed the atomic positions that of Jahn–Teller minima (green dots) and barrier saddle points (blue dots) C 1 h by means of the HSE06 functional. Additionally, we obtained the C 3 v case (red dot) by enforcing symmetry constraints only on atomic positions and thus electronic correlation effects that of Eq. (9) can appear. Finally, we interpolated the geometries and calculated the electronic structure along the purple line to ensure the continuous shape of the APES. Additionally, we depict purple route in (b) inset that shows the raw HSE06 data points while (a) depicts the APES defined by Eq. (9), thus X and Y are dimensionless therein. We note that only the green and blue data points that of (a) follow the axis on the top while the purple data points do not.

FIG. 3.

Jahn–Teller effect for the ground state of NV 0. (a) Adiabatic potential energy surface (APES) of NV 0 defect’s | 2 E ground state. We relaxed the atomic positions that of Jahn–Teller minima (green dots) and barrier saddle points (blue dots) C 1 h by means of the HSE06 functional. Additionally, we obtained the C 3 v case (red dot) by enforcing symmetry constraints only on atomic positions and thus electronic correlation effects that of Eq. (9) can appear. Finally, we interpolated the geometries and calculated the electronic structure along the purple line to ensure the continuous shape of the APES. Additionally, we depict purple route in (b) inset that shows the raw HSE06 data points while (a) depicts the APES defined by Eq. (9), thus X and Y are dimensionless therein. We note that only the green and blue data points that of (a) follow the axis on the top while the purple data points do not.

Close modal
Single configuration that resembles | 2 A 2 at largest extent that of Eq. (5) is the
(11)
determinant. In simple words, according to Eq. (11), the A | a 1 e e + configuration accessible within our DFT method corresponds to 66% of | 2 A 2 and 33% of | 4 A 2 character, respectively. In order to approximate the energy of multiconfigurational | 2 A 2 , we utilize a correction scheme similar to the method of von Barth originally proposed in 1979 (see Ref. 69). We multiply Eq. (11) in the left side by a 1 e e + | A H ^ e that results in E DFT ( A | a 1 e e + ) = a 1 e e + | A H ^ e A | a 1 e e + = ( 1 3 4 A 2 + 1 2 | 2 3 2 A 2 + | ) H ^ e ( 1 3 | 4 A 2 + 1 2 2 3 | 2 A 2 + ) = 1 3 E DFT ( | 4 A 2 ) + 2 3 E ( | 2 A 2 ), where H ^ e is the full electronic Hamiltonian and off diagonal matrix elements 2 A 2 | H ^ e | 4 A 2 are zeroes. Therefore, we can approximate the energy of multiconfigurational | 2 A 2 excited state by the following scheme:
(12)
where E DFT ( A | a 1 e + e ) = 1.22, 1.20, 1.74 eV with LDA, PBE, HSE06 functionals, respectively, and we allowed the atomic positions to relax and E DFT ( | 4 A 2 ) = 0.86 , 0.58 , 0.48 eV that of Table I. We note that this is only a crude approximation; nevertheless, we successfully applied this assumption in Eq. (6) of Ref. 70 and in Eq. (A5) in Ref. 52. In this case, we approximate the | 2 A 2 level at 2.36 eV by means of HSE06, which experiences a correction on its Stokes-shift49 by Δ E stokes = ω [ S ( | 2 A 2 ) S ( A | a 1 e e + ) ] = 0.07 eV, where S ( | 2 A 2 ) = 1.91, S ( A | a 1 e e + ) ] = 1.16, and ω = 93.9 meV that yields 2.22 eV ZPL energy close to the experimentally observed data at 2.156 eV.
Similar to considerations that of Sec. III C for | 2 A 2 , we estimate the energy of | 2 E excited state by left multiplying A | a 1 e x e y = 1 2 | E 3 2 + + 1 3 | 4 A 2 + + 1 6 | 2 A 2 + by a 1 e x e y | A H ^ e. Therefore, the approximated energy of | 2 E will be
(13)
where E DFT ( | a 1 e x e y ) = 1.26, 1.24, 1.72 eV for LDA, PBE, HSE06 functionals, respectively, where we allowed atomic positions to relax. We used Eq. (13) to determine excitation energies for Table I.

Upon depicting the electronic structure of the | 2 E ground state, we aim to describe the observed fine structure of the | 2 E | 2 A 2 optical transition. We analyze the adjoint Jahn–Teller (Sec. IV A) and spin–orbit coupling (SOC) (Sec. IV B) effects. Additionally, local and external strain around the defect may affect the optical lines. Therefore, we will discuss both effects to approximate the λ effective SOC parameter that splits the | 2 E multiplet into two electron–phonon coupled doublets: | E ~ 1 2 , | E ~ 3 2 . Finally, we will discuss the possible excitation schemes from the ground state in Sec. IV D.

The ground state of NV 0 is a Jahn–Teller (JT) system where the | 2 E degenerate orbital interacts with a degenerate E vibration mode as the double degenerate e orbital is occupied by a single electron. We map the adiabatic potential energy by means of HSE06 functional to derive the JT Hamiltonian as
(14)
acting between orbital σ ^ z = | 2 E x 2 E x | | 2 E y 2 E y | = ( 1 1 ), σ ^ x = | 2 E x 2 E y | + | 2 E y 2 E x | = ( 1 1 ) and vibrational X ^ = ( a X + a X ) / 2, Y ^ = ( a Y + a Y ) / 2 degrees of freedom. We determined the ω E frequency of the harmonic oscillator by a parabola fit in the adiabatic potential energy surface (APES). We formulated F and G electron–phonon interaction strengths71,72 from JT distortion energy as E J T = { F 2 / [ 2 ( ω E 2 G ) ] } and the barrier energy as δ JT = [ ( 4 E JT G ) / ( ω E + 2 G ) ]. The computed first vibronic level ( 14.0 meV) is in excellent agreement with the experimentally observed one21 at 13.6 meV (see Table II). Therefore, we expect that our reported the p Ham factor is reliable24,71,72 that is used throughout our study.
TABLE II.

Jahn–Teller (JT) coupling parameters for the |2E⟩ ground state by means of various DFT functionals: LDA, PBE, HSE06.

LDA73 LDAp.w.PBEp.w.HSE06p.w.expt.
ω E62.5 64.5 63.5 62.5  meV 
EJT73.2 61.2 49.9 30.1  meV 
δJT10 9.9 6.0 20.7  meV 
tun.a21.4 25.5 29.2 14.0 13.6b meV 
p factor: 0.118c 0.161 0.199 0.304   
LDA73 LDAp.w.PBEp.w.HSE06p.w.expt.
ω E62.5 64.5 63.5 62.5  meV 
EJT73.2 61.2 49.9 30.1  meV 
δJT10 9.9 6.0 20.7  meV 
tun.a21.4 25.5 29.2 14.0 13.6b meV 
p factor: 0.118c 0.161 0.199 0.304   
a

Energy of first vibronic excited state of JT system depicted by Eq. (14) known as tunneling splitting.71,72

b

Experimentally observed via uniaxial stress measurements21 on the ZPL where 110 cm−1 fine structure was reported. Baier et al. reported this splitting as 12(2) and 13(4) meV by fitting thermal activation energy for an Orbach process.16 

c

We applied Eq. (14) on data found in Ref. 73 to calculate the Ham reduction factor as expressed in Eq. (17).

Spin–orbit coupling (SOC) arises due to parallel or antiparallel spin and orbital momenta within | e , orbitals. Therefore, the single particle spin–orbit operator can be depicted as
(15)
where we show both parallel ( l ^ z | e ± = ± | e ± ) and perpendicular SOC components that flip the orbital character as l ^ ± | a 1 = 2 | e ± .74 We depict the electron spin by s ^ x = 1 2 ( 1 1 ), s ^ y = 1 2 ( ( i i ), s ^ z = 1 2 ( 1 1 ), s ^ ± = s ^ x ± i s ^ y standard spin operators.
One may recognize that the net effect of two e orbitals are completely canceled out for | E 1 2 / 3 2 configurations [Eq. (6)], thus it is unaffected by any spin–orbit coupling. Therefore, only the | E 1 2 / 3 2 configurations of Eq. (1) are subject splitting by λ 0 as
(16)
Therefore, the spin–orbit coupling for the correlated | E 1 2 / 3 2 SOC parameter will be partially quenched by the c 2 correlation parameter of Eq. (10) and p Ham reduction factor24,71,72 simultaneously as
(17)
where c 2 is the probability that the system is in the | 2 E configuration and p is the Ham reduction factor of the JT effect.

First, we determined the λ 0 SOC parameter by constraining a half–half electron into | e ± Kohn–Sham orbitals and leave | e ± orbitals empty. This way, we completely cancel out the symmetry breaking effect both JT and electronic correlation origins thus the wavefunctions will be complex valued, that is, | e ± , where only spin–orbit coupling can split the degeneracy. We approximated the spin–orbit coupling parameters as λ 0 = 290 (318) GHz with HSE06 (PBE) functional depicting negligible difference on the choice of the DFT functional.

We evaluate Eq. (16) by means of HSE06 functional, thus we used c 2 = 0.71, p = 0.304 and λ 0 = 290 GHz that results in λ t h e o r y = 31.4 GHz that is significantly larger than that of the experimentally derived ones at λ e x p t . = 4.48 GHz (Ref. 14) or 9.8 GHz (Ref. 16). The p and c 2 factors cannot correspond to this degree of discrepancy. Therefore we conclude that the ( a 1 1 a 1 1 e + 0.5 e 0.5 e + 0 e 0 ) constrained occupation that we forced upon our SOC calculations is not a reliable method. However, the broken symmetry solution of Eq. (9) cannot be used to this end as a and e orbitals are mixed in that case, which also brings λ l ^ ± s ^ matrix elements.

Nevertheless, we attempted to determine λ 0 from the correlated Eq. (9) wavefunction. Therefore, we determine λ 0 from the SOC matrix elements | u a 1 + v e x | H ^ S O C | e y | = v λ 0 λ 0 = 71.1 ( 74.7 ) GHz between unoccupied and occupied orbitals by means of HSE06 (PBE) functionals. However, it still overestimates experimental data as the reduced value is p c 2 λ 0 = 15.4 GHz with HSE06 functional. We note that the λ 0 value might be estimated by the non-correlated e orbital in the 3 E excited state of NV , which results in λ 0 = 31.6 GHz with HSE06 functional (Ref. 24). From this estimation the final spin–orbit gap yields p c 2 λ 0 = 6.8 GHz. In summary, our results cannot accurately determine the spin–orbit gap from the Kohn–Sham DFT level of theory because of the highly correlated nature of the | 2 E ground state.

Lastly, we note that we determine λ as a 1 | H ^ SOC | e = 2 λ = 2 × 107.3 GHz matrix element as calculated from the ( a 1 1 a 1 1 e + 0.5 e 0.5 e + 0 e 0 ) constrained occupation. We use this λ value for the ISC transition | 2 A 2 | 2 E (see Sec. V for details).

We determine the hyperfine and quadrupolar parameters to a central 14N nucleus as we described in Ref. 75. We used an increased (600 eV) plane wave basis for electron wavefunctions during the determination of spin parameters with HSE06 functional. The spin Hamiltonian that depicts the motion of 14N nuclear spin ( I = 1) goes as
(18)
where S and { σ ^ z = | 2 E x 2 E x | | 2 E y 2 E y | = ( 1 1 ), σ ^ x = | 2 E x 2 E y | + | 2 E x 2 E y | = ( 1 1 ) } depict the electron spin and orbital degrees of freedom, respectively. We show the hyperfine and quadrupolar parameters in Table III. We used the ( a 1 1 a 1 1 e x 0.5 e y 0.5 e x 0 e y 0 ) single electron occupation to determine the hyperfine parameters to avoid broken-symmetry spin densities. However, the electric field gradient related to the nuclear quadrupole interaction does not depend on the broken symmetry spin density but rather on the local potential that is a scalar quantity. Therefore, in contrast to the calculation of hyperfine tensors, we use the broken symmetry occupation of Eq. (9) to determine the quadrupolar tensor of 14N. By doing so the quadrupolar parameters become reduced upon removing the symmetry constraints from 5.2 MHz, 10.1 kHz, 11.0 kHz to 4.8 MHz, 9.6 kHz, 9.7 kHz for Q, Q 1, Q 2 matrix elements, respectively. However, the hyperfine parameters undergo an artificial change upon symmetry constraint removal: 2.8, 3.2, 0.1, 0.1 MHz parameters become 9.9, 6.3, 1.0, 0.8 MHz for A , A , A 1, A 2, respectively, due to the spin contamination error. Therefore, we use the symmetrical spin density results for hyperfine parameters in Table III for the sake of simplicity.
TABLE III.

Theoretical hyperfine and quadrupolar parameters for 14N within the |2 E⟩ ground state.

A = 2.8 MHz A ˭ −3.2 MHz 
A1 ˭ 0.1 MHz A2 ˭ −0.1 MHz 
Q ˭ −4.8 MHz  
Q1 ˭ 9.6 kHz Q2 ˭ −9.7 kHz 
A = 2.8 MHz A ˭ −3.2 MHz 
A1 ˭ 0.1 MHz A2 ˭ −0.1 MHz 
Q ˭ −4.8 MHz  
Q1 ˭ 9.6 kHz Q2 ˭ −9.7 kHz 
The dynamic JT also affects the hyperfine interaction of 13C nuclear spins that might be present around the defect with 1.1% natural abundance in diamond. Three carbon atoms with dangling bonds near the vacancy exhibit the largest hyperfine couplings (see Fig. 4). The spin Hamiltonian for such states can be constructed76 as follows:
(19)
where A is the symmetric part and A x, A y are the non-symmetric parts, and S = ( S ^ x S ^ y S ^ z ) T and I = ( I ^ x I ^ y I ^ z ) depict the electron and nuclear spin vectors, respectively. Similar to that of the central 14N nuclear spin, we determined the hyperfine matrices by means of the HSE06 hybrid functional on the ( a 1 1 a 1 1 e x 0.5 e y 0.5 e x 0 e y 0 ) electronic configuration for which we report the following three:
(20)
hyperfine matrices. We note that the carbon atoms with C N bonds and all other second, third, etc., neighbors exhibit negligible < 1 MHz hyperfine constants except for the three equivalent carbon atoms that are connected to the carbon atoms with dangling bonds that we show in Fig. 4(b) with gray lobes. These three 13C nuclear spins exhibit substantial isotropic hyperfine interaction ( W ^ = I A iso S effective Hamiltonian), e.g., only A is nonzero that of Eq. (20) and its diagonal is filled with the same A iso = 10 MHz parameter, and all other parameters of Eq. (19) are less than 1 MHz.
FIG. 4.

Atomic positions of the NV 0 defect from two different orientations. We show two different viewing directions: (a) depicts the system viewed along the [111] trigonal symmetry axis, while (b) depicts a viewing angle perpendicular to it. We depict the [111] oriented coordinate system in (i) used within Eqs. (18) and (20) for hyperfine and quadrupolar parameters larger than > 1 MHz. Additionally, we show the Cartesian coordinate system of diamond in (ii). We depict the | e x orbital by the blue/red colors for negative/positive parts of the wavefunction.

FIG. 4.

Atomic positions of the NV 0 defect from two different orientations. We show two different viewing directions: (a) depicts the system viewed along the [111] trigonal symmetry axis, while (b) depicts a viewing angle perpendicular to it. We depict the [111] oriented coordinate system in (i) used within Eqs. (18) and (20) for hyperfine and quadrupolar parameters larger than > 1 MHz. Additionally, we show the Cartesian coordinate system of diamond in (ii). We depict the | e x orbital by the blue/red colors for negative/positive parts of the wavefunction.

Close modal

As mentioned in the Methods section (Sec. II), we determine the optical sidebands of PL and ABS spectra by means of the Huang–Rhys method.23,53 First, we calculated the Huang–Rhys sideband by the geometry differences between the A | a 1 e e + configuration and the | 2 E correlated ground state as depicted in Eq. (9) with all symmetry constraints lifted. However, the predicted S = 1.16 Huang–Rhys factor significantly underestimates the intensity of the PL sideband as plotted in Fig. 2(b).

Therefore, we utilize the correction scheme by means of Eq. (11) as follows. The derivatives against all X i ( i = 3 N 3) configuration coordinates
(21)
gives an approximation for forces acting on nuclei for the | 2 A 2 state. In the case of linear electron–phonon coupling, the forces are proportional to the relaxed atomic displacements [see Eq. (7) in Ref. 23]. Therefore, the optimized geometry of | 2 A 2 can be estimated as
(22)
where X i ( Γ ) are relaxed atomic positions within electronic state Γ. We find that the correction scheme improves not only the ZPL energy (see Table I), it also predicts a PL sideband with S = 1.91 HR factor as depicted in Fig. 2(b).

The absorption spectrum is much broader than the PL spectrum for the NV 0 defect. We attribute this phenomenon to the presence of the | 2 E second excited state close to the | 2 A 2 first excited state. Furthermore, the ionization bands emerge at about 2.7 eV (e.g., Ref. 48), which overlaps with the phonon sideband of the | 2 A 2 related absorption. In Sec. IV E, we analyze the optical transition dipole moments toward the doublet excited states from the ground state.

We determine the optical transition dipole moments between occupied–unoccupied Kohn–Sham state pairs as implemented in the VASP package.77 For symmetry reasons, we placed half–half electrons on | e x , y orbitals of Fig. 1(a), thus populating the ( a 1 1 a 1 1 e x 0.5 e y 0.5 e x 0 e y 0 ) configuration to avoid broken-symmetry solutions shown in Eq. (9). We calculate the transition dipole moments within the spin majority channel | a 1 | d ^ | e x , y | = 4.22 Debye and | e x | d ^ | e y | = d / 2 = 1.81 Debye. Next, we determine the transition dipole moments acting between multiplets of Eqs. (5), (7), and (10) as follows:
(23)
where we show the dipoles for c = 0.84 and d = 0.54 taking the correlation effects into account whereas the raw data without correlation ( c = 1.0 and d = 0.0) are given in parentheses. Additionally, we give a rough estimate for d by determining the transition strength toward the [111] direction for the broken-symmetry solution as given in Eq. (9). Therefore, we evaluate the dipole integral from the occupied | u a 1 + v e x toward the unoccupied third spin- orbital as d = 1.07 Debye.
The radiative lifetime of dipole transitions ( d) can be evaluated78 as τ = 3 π ε 0 n ω 3 | d | 2 = α | d | 2, where n = 2.42 is the refractive index of diamond, ω 2.156 eV is the excitation energy and ε 0 is the permittivity of vacuum. Finally, α = 501.15 Deby e 2 ns is the conversion factor from Debye to nanosecond units. This allows us to evaluate the radiative lifetimes as
(24a)
(24b)
where the data in parentheses are calculated by neglecting electronic correlation effects. We note that there is an extra 1 2 factor for τ ( 2 A 2 ) because | 2 E is double degenerate unlike | 2 A 2 .

The calculated 6.6 ns radiative lifetime overestimates the optical transition’s strength as the experimentally observed excited state lifetime varies around 13–20 ns (Refs. 28–31). This result is surprising as our ab initio method correctly predicts radiative lifetime of NV excited state lifetime [13.7 ns in comparison to the experimentally observed 12.0 ns (Ref. 79)] by determining d N V = | a 1 | d ^ | e x , y | = 4.98 Debye. We tentatively explain the discrepancy with the complex nature of the electronic excited states. We conclude from our ab initio calculations that the bright excited state of NV 0 does not have a pure | 2 A 2 character because electron–phonon coupling can mix | 2 A 2 and | 2 E excited states with each other as the two energy levels are very close to each other. In such a case if the coupling between them is strong it could also explain why the second ZPL for | 2 E is missing in ABS as the two electronic excited states are merged into common vibronic bands. The calculation of the electron–phonon coupling, i.e., the pseudo-Jahn–Teller effect between | 2 A 2 and | 2 E is beyond the scope of this study as it is extremely difficult to achieve a reliable adiabatic potential energy surfaces and energy spacing for the two excited state doublets at high accuracy.

In summary, we suggest polarization dependent measurements both for the sideband both in ABS and PL to resolve this issue. The optical polarization already been resolved for the ZPL absorption80 that points to a | 2 E | 2 A transition, whereas later measurements confirmed34 that the excited state is a | 2 A 2 multiplet. In particular, an | 2 A 2 multiplet active through “ x” and “ y” directional d dipoles only, whereas | 2 E additionally active through the d dipole pointing toward “ z.” However, we propose that the optical sidebands contain not only the | 2 A 2 character, thus “z” dipole activity should be present toward the | 2 E multiplet at higher energies in the absorption spectrum. The variation of the observed optical lifetimes of NV 0 with various excitation techniques in bulk diamond might be also explained by the intertwined | 2 A 2 and | 2 E character of the excited state.

The | 4 A 2 is an orbitally non-degenerate state, thus the spin–orbit interaction is zero within | 4 A 2 , but it could occur between the doublets and | 4 A 2 as a second order interaction. The only doublet state close in energy is the ground state | 2 E , whereas the other states can be ignored because of large energy spacing. We introduce the SOC matrix elements between orbitals that are responsible for ISC. First, we determine the possible SOC matrix elements that connect | 4 A 2 to the ground state. One may notice that only the ( a a e ) configuration of | 2 E is connected by SOC, i.e., 4 A 2 | H ^ SOC | 2 E 0 because the 4 A 2 | H ^ SOC | 2 E = 0 for ( a e e ) configuration as two-particle excitation would be required to transform ( a e ± e ) configuration into ( a e ± e ± ), cf. Eqs. (1) and (6) for details. Therefore, the SOC matrix elements are
(25)
The first two equations in Eq. (25) are important to consider in the context. To summarize, in terms of electronic structure, the relaxation into | E 3 2 is thrice faster than that into | E 1 2 . Next, we use the HR theory as described in Refs. 24, 79, 81, 82, and 83 to determine the ISC rates that can be expressed as
(26)
where F is the Huang–Rhys spectral function between | 4 A 2 and | 2 E levels, and Δ is the energy difference between the two levels. We obtained F function by relaxing the atomic positions within | 4 A 2 and | 2 E configurations. We note that initially we applied C 3 v symmetry constraint on the geometry of the DFT calculation that results in HR factor S A 1 = 0.90 from totally symmetric ( A 1) modes. When we removed all symmetry constraints and relaxed to the JT minimum in Fig. 3, the HR factor increased to S = S A 1 + S E = 1.53 that includes the additional effect of E phonons that we depict the resulting F HR ( ω ) in Fig. 5 graphically.
FIG. 5.

Intersystem crossing (ISC) from the metastable 4 A 2 state of the NV 0 defect. (a) Mechanism of the ISC process. (b) and (c) ISC process within linear and log scale, respectively. We depict the theoretically viable region by green stripes in (c). We take an upper bound for the ISC rate ( τ EPR 1) from EPR measurements (Ref. 38, see Sec. V for details).

FIG. 5.

Intersystem crossing (ISC) from the metastable 4 A 2 state of the NV 0 defect. (a) Mechanism of the ISC process. (b) and (c) ISC process within linear and log scale, respectively. We depict the theoretically viable region by green stripes in (c). We take an upper bound for the ISC rate ( τ EPR 1) from EPR measurements (Ref. 38, see Sec. V for details).

Close modal

Additionally, we determine the spectral function for ISC by means of the JT theory: F ( ω ) = i c 0 , ( i ) 2 F A 1 ( ω ε i ), where F A 1 ( ω ) is the HR spectral function, exhibiting S A 1 = 0.90 HR factor and ε i/ Ψ i is the eigenspectrum/wavefuncion of JT Hamiltonian of Eq. (14). At zero temperature, there would be no E-phonons present within the metastable quartet, thus its electronic plus E x E y–phonon wavefunction would be a 4 A 2 | 0 , 0 | vacuum state for vibrations. Any | Ψ i can be expressed as | Ψ i = c 00 ( i ) | 2 E m j | 0 , 0 + n , m , m j n + m 1 c n , m , m j ( i ) | 2 E m j | n , m within the Born–Oppenheimer basis, where | n , m = ( a X ) n ( a Y ) m n ! m ! | 0 , 0 . Therefore, it can be concluded that the ISC process cannot target | Ψ i JT levels transforming as A 1 or A 2 representations because their c 00 ( i ) expansion coefficient is zero due to group theory constraints and | n , m s are thermally not activated. ISC toward to the lowermost four | 2 E m j relaxes with the c 00 ( 1 4 ) 2 = 0.516 coefficient and m j is conserved fully from the initial | 4 A 2 . However, if the system relaxes to the first vibronic E-level ( | Ψ 1 4 ) at ε 7 10 = 47 meV first by c 00 ( 7 10 ) 2 = 0.280 coefficient, a very fast relaxation will occur back to | 2 E m j or toward the | Ψ 5 6 JT tunneling state at ε 5 6 = 3 Γ tun . = 14 meV. We note that here m j = ± 1 2 3 2 transfers are allowed, and we expect that orbital relaxation84 would be the dominating process in the fs region.85 Nevertheless, the system would relax to | 2 E ± 1 2 from | 2 E 3 2 within μ s due to the same orbital relaxation.

Surprisingly, the correction from the JT theory is negligible. We calculated the ISC process fully by means of the HR theory only: F ( ω ) F HR ( ω ) = F A 1 ( ε ) F E ( ω ε ) d ε, we obtained almost the same transition rate as that of the JT theory [see Fig. 5(c)]. We used two estimates for λ . First, we used our λ = 107.3 GHz ab initio result as we discussed in Sec. IV B. However, our previous study about ISC processes for NV indicated that the off diagonal SOC matrix elements may overestimate24,82 the strength of the transition. Therefore, we used the previous result only as an upper bound, and we used λ = 1.2 × λ z NV , where λ z NV is the SOC parameter of NV in its | 3 E excited state as suggested by previous studies.24,79,81,82 Finally, we use the HSE06 c 2 = 0.71 correlation parameter as discussed in Sec. III B. We note that HSE06 DFT usually predicts the position non-correlated states within ± 0.1 eV precision. Therefore, we estimate the position of | 4 A 2 level as 0.48 ± 0.1 eV (see Table I). However, we note that configurational interaction wavefunction methods predict | 4 A 2 at 0.68 eV, so we use 0.7 eV as an upper bound for Δ. Therefore, all the necessary parameters are known to evaluate Eq. (26). By taking all the broad energy intervals into account, our theory predicts the lifetime of | 4 A 2 ± 1 2 within 40 μ s 13 ns within the Δ = 0.7 eV 0.4 eV interval (see Table I).

An another aspect of | 4 A 2 is that it is EPR active. As we mentioned in the Introduction, it can be occupied with ionizing NV . However, the ISC rate cannot be too fast, otherwise the optically induced EPR could not be observed. Indeed, we fitted Lorentzian peaks on the first top-left experimentally measured EPR spectra from Fig. 1(a) in Ref. 38 recorded at T = 10 K. We estimated the linewidth of the Lorentzian broadening as 0.068(3) mT that corresponds to 1.8(1) MHz in EPR frequency. Therefore, we can interpret the linewidth as natural line broadening to estimate the spin lifetime of | 4 A 2 as τ EPR 1 2 π × 1.8 MHz = 89 ( 3 ) ns. However, the broadening may arise from other sources such as hyperfine interaction broadening due to distant 13C isotopes, for example. Thus, we use this value only as a minimum estimation for Γ 1 2 1 [see the green area in Figs. 5(b) and 5(c)]. Our final conclusion is that the lifetime of | 4 A 2 is in the μs regime at cryogenic temperatures. We summarize the level structure of NV 0 and the possible transitions in Fig. 6. It is worthy to consider the route on which the system may end up in | 4 A 2 . However, it can be seen that the | 4 A 2 is separated from | 2 A 2 by 1.6 eV, thus ISC transition is negligible between the states. Therefore, in Sec. VI, we will present how photoionization of NV can populate the metastable spin quartet state, | 4 A 2 .

FIG. 6.

Summary of ISC rates and level structure of NV 0. We take the experimentally observed D zero-field splitting value from Ref. 38. Additionally, we expect that | 4 A 2 spin relaxation time is close to that of | 3 A 2 of NV because both multiplets are similar orbital singlets.

FIG. 6.

Summary of ISC rates and level structure of NV 0. We take the experimentally observed D zero-field splitting value from Ref. 38. Additionally, we expect that | 4 A 2 spin relaxation time is close to that of | 3 A 2 of NV because both multiplets are similar orbital singlets.

Close modal

The NV defect is often observed under illumination. The most typical photoexcitation of NV defect occurs with green laser (e.g., at 532 nm in wavelength or 2.33 eV in energy), which can lead to charge switching between the negatively charged and neutral NV defects.

Understanding of photoionization of NV defect is difficult because of the possible non-linear optical effects such as successive two-photon absorption processes via real excited states and the presence of long-living | 1 E in NV and | 4 A 2 in NV 0. It is very challenging to calculate the photoionization threshold energies and rates at ab initio level for two reasons: (i) the singlet states in NV cannot be accurately calculated by DFT and (ii) simple ionization and Auger-ionization might compete where the latter is a two-body interaction, an inherently complex task for ab initio methods. Issue (i) can be principally solved by embedding wavefunction techniques that was applied to the ionization of | 1 E of NV .86 We note that the photoionization cross section in that study was calculated within the Γ-point in a 512-atom supercell, which is not fully convergent because the continuous band region at elevated energies above the ionization threshold energy requires much higher k-point sampling87 and band unfolding methods.42 Nevertheless, the calculated photoionization threshold energy at 2.2 eV implies that direct ionization, NV + h ν NV 0 + e C B M , occurs under green illumination. Issue (ii) was dealt with representing the appropriate exciton wavefunctions by combination of the respective Kohn–Sham wavefunctions in the Coulomb-integrals associated with the Auger-ionization rates.88 The method was applied to the | 3 E of NV ionization to the ground state of NV 0 and resulted in about 1 ns inverse rate.88 This result was later criticized with the argument that the electron density in the applied 512-atom supercell was very high so the Auger-ionization rate was orders of magnitude overestimated, which should depend on the carriers density.42 On the other hand, no other ab initio data have been yet published about the Auger-ionization rates of the NV defect to justify this argument.

In this study, we proved that even the ground state of NV 0 is highly correlated that was not considered in the previous studies in the context of photoionization cross section. Because of the high complexity of the multiplet states of NV 0 and numerical challenges of calculating direct and Auger-ionization rates, we rather show the estimated photoionization threshold energies for the NV + h ν NV 0 + e C B M and NV 0 + h ν NV + h V B M + processes (see Fig. 7). Furthermore, we use the Slater–Condon principles to identify the feasible ionization processes as alluded in Refs. 42 and 89 for NV + h ν NV 0 + e C B M and NV 0 + h ν NV + h V B M +, respectively. Briefly, the Slater–Condon selection rule implies that direct ionization and Auger-ionization only occurs for changing maximum one and two single particle wavefunctions between the initial and final multiplet states, respectively.

FIG. 7.

Summary of photoionization processes of NV 0. We round the related transition rates, thus we depict the lifetimes as rough ps, ns, and μs estimates.

FIG. 7.

Summary of photoionization processes of NV 0. We round the related transition rates, thus we depict the lifetimes as rough ps, ns, and μs estimates.

Close modal

We start with the simplest photoionization scheme, i.e., the single photon absorption from the initial ground state to the final ground state. This yields | 3 A 2 + 2.6 eV | 2 E and | 2 E + 2.8 eV | 3 A 2 where the first process was observed at room temperature90 for which the T = 0 K DFT calculation predicted 2.7 eV.48,87 These energies fall in the blue wavelength region. We note that the ground state | 2 E is a correlated multiplet [see Eq. (10)] of 71% | 2 E character, so the direct ionization process is partially quenched as the Slater–Condon selection rule prohibits the direct ionization toward | 2 E with 28% contribution to the | 2 E ground state of NV 0 that we depict by dashed arrows in Fig. 7. On the other hand, the shelving state of NV , | 1 E , can be ionized by green light as | 1 E + 2.2 eV | 2 E where no quenching occurs due to the correlated nature of | 2 E . We find that the shelving state of NV 0, | 4 A 2 behaves very similarly as | 4 A 2 + 2.3 eV | 3 A 2 , which is fully allowed by the direct ionization process. Finally, we note that ionization from | 1 A 1 is unlikely due to its very short lifetime of 100 ps.91 

As was also pointed out in Ref. 42, intersystem crossing is strongly hindered between | 2 A 2 and | 4 A 2 because of selection rules and a large energy spacing between these levels so occupation of the shelving | 4 A 2 could occur via photoionization from NV . Direct ionization from | 3 A 2 yields an ultraviolet threshold energy for this process, but then it would immediately convert it back to NV . However, it was found in an EPR measurement38 that green illumination for NV ensemble diamond samples can effectively pump the system into | 4 A 2 . This may be explained two-photon absorption from | 3 A 2 to | 4 A 2 via | 3 E . In this process, the first photon is absorbed to go from | 3 A 2 to | 3 E (usual neutral excitation) and the second photon is quickly absorbed before it decays within NV . In Ref. 42, they argued that the Slater–Condon rules dictate that direct ionization occurs only toward | 4 A 2 with a threshold energy at 1.2 eV, which explains the enhancement of ionization at this energy in a two-color excitation experiment.92 Furthermore, it could also explain a spin polarization of the | 4 A 2 toward m s = ± 1 2 over m s = ± 3 2.42 Indeed, this phenomenon can be observed in Fig. 1(b) of Ref. 38 where the 3 2 1 2 transition (1st spectrum) exhibits induced emission (positive Lorentzian), whereas the + 1 2 + 3 2 transition (4th spectrum) can be identified as absorption (negative Lorentzian). Our analysis shows that the lifetime of | 4 A 2 can be indeed μs long in dark once the state has been occupied. However, green illumination can photoionize it back to the ground state of NV according to our results. Thus, this is a rather complex process where the average lifetime of | 4 A 2 depends on many factors, e.g., the laser power and the presence of other defects that may induce additional carriers (electrons and holes) upon illumination.

We analyze further the two-photon ionization processes. Our calculations on NV 0 revealed that direct ionization channel exists between | 3 E and | 2 E because of the 28% | 2 E character in the latter. Thus, an additional direct ionization process occurs from 0.7 eV from | 3 E beside the already considered | 3 E + 1.2 eV | 4 A 2 process in Ref. 42. The Auger-ionization process from | 3 E is fully active toward | 2 E . Thus, the sum of direct ionization and Auger-ionization rates will give the total ionization rate. Next, we consider the two-photon ionization processes via | 2 A 2 . The | 2 A 2 + 0.6 eV | 3 A 2 process is very similar to the reverse ionization process where this process is partially quenched by the electron–phonon related mixture of | 2 E into | 2 A 2 similar to that of the quenched radiative lifetime of | 2 A 2 (see Sec. IV E). Additionally, the | 2 A 2 + 1.0 eV | 1 E process is also forbidden at linear order at the Slater–Condon level. However, it becomes partially allowed due the electron–phonon related mixture of | 2 E or the correlated nature of | 1 E (see Ref. 82). We note again that the Auger-ionization process is fully allowed between these two states. Auger-ionization also can connect the | 2 A 2 and | 2 A 1 multiplets, and direct ionization is only allowed in the second order since it requires | 2 E character to be mixed in that is only possible through simultaneous electron–phonon interaction and electronic correlation effects within | 2 A 2 . We argued in a previous study89 that the Auger-ionization rate can be substantial because of the presence of a deep resonant defect state of NV 0 in the valence band, thus we argue that | 2 A 2 + 2.2 eV | 1 A 1 indeed occurs.89 Finally, ionization may occur from | 2 A 2 toward | 3 E upon blue illumination starting at 2.6 eV, thus this route is blocked at the typical green illumination of the NV defect.

The electronic structure of NV 0 was determined by DFT HSE06 calculations where the correlation of orbitals is taken into account with the guide of the group theory, in order to induce correction in the total energy and the geometry in the excited states. We find that the electron–phonon coupling is very strong and it affects both the fine structure of the ground state as well as the radiative lifetime of the excited state. In particular, we find that the 2 E and 2 A 2 excited states are intertwined by symmetry breaking phonons, and the optical properties can be only explained by invoking strong electron–phonon coupling between the doublet excited states. We also analyzed the lifetime of the 4 A 2 metastable state and provided a detailed analysis on the intersystem crossing between the metastable quartet state and the doublet ground state. Finally, the calculated electronic structure of NV 0 yields new insights into the photoionization processes of the NV defect.

We thank L. Razinkovas, R. Ulbricht, and C. Linderälv for fruitful discussions. Support by the National Research, Development and Innovation Office (NKFIH) as well as by the Ministry of Culture and Innovation within the Quantum Information National Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) is much appreciated. A.G. acknowledges the high-performance computational resources provided by KIFÜ (Governmental Agency for IT Development) Institute of Hungary, the European Commission for the projects QuMicro (Grant No. 101046911) and SPINUS (Grant No. 101135699), and the QuantERA II project Maestro (NKFIH Grant No. 2019-2.1.7-ERA-NET-2022-00045). G.T. was supported by the János Bolyai Research Scholarship of the Hungarian Academy of Sciences.

The authors have no conflicts to disclose.

Gergő Thiering: Conceptualization (equal); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). Adam Gali: Conceptualization (equal); Funding acquisition (lead); Investigation (equal); Project administration (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

1.
L.
du Preez
, Ph.D. thesis, University of Witwatersrand, Johannesburg, 1965.
2.
M. W.
Doherty
,
N. B.
Manson
,
P.
Delaney
,
F.
Jelezko
,
J.
Wrachtrup
, and
L. C.
Hollenberg
,
Phys. Rep.
528
,
1
(
2013
).
3.
A.
Gruber
,
A.
Drabenstedt
,
C.
Tietz
,
L.
Fleury
,
J.
Wrachtrup
, and
C. V.
Borczyskowski
,
Science
276
,
2012
(
1997
).
4.
F.
Jelezko
,
T.
Gaebel
,
I.
Popa
,
A.
Gruber
, and
J.
Wrachtrup
,
Phys. Rev. Lett.
92
,
076401
(
2004
).
5.
T. D.
Ladd
,
F.
Jelezko
,
R.
Laflamme
,
Y.
Nakamura
,
C.
Monroe
, and
J. L.
O’Brien
,
Nature
464
,
45
(
2010
).
6.
D. D.
Awschalom
,
L. C.
Bassett
,
A. S.
Dzurak
,
E. L.
Hu
, and
J. R.
Petta
,
Science
339
,
1174
(
2013
).
8.
M.
Pompili
,
S. L. N.
Hermans
,
S.
Baier
,
H. K. C.
Beukers
,
P. C.
Humphreys
,
R. N.
Schouten
,
R. F. L.
Vermeulen
,
M. J.
Tiggelman
,
L.
dos Santos Martins
,
B.
Dirkse
,
S.
Wehner
, and
R.
Hanson
,
Science
372
,
259
(
2021
).
9.
V.
Acosta
and
P.
Hemmer
,
MRS Bull.
38
,
127
(
2013
).
10.
K.
Jensen
,
N.
Leefer
,
A.
Jarmola
,
Y.
Dumeige
,
V.
Acosta
,
P.
Kehayias
,
B.
Patton
, and
D.
Budker
,
Phys. Rev. Lett.
112
,
160802
(
2014
).
11.
S. I.
Bogdanov
,
M. Y.
Shalaginov
,
A. S.
Lagutchev
,
C.-C.
Chiang
,
D.
Shah
,
A. S.
Baburin
,
I. A.
Ryzhikov
,
I. A.
Rodionov
,
A. V.
Kildishev
,
A.
Boltasseva
, and
V. M.
Shalaev
,
Nano Lett.
18
,
4837
(
2018
).
12.
M.
Gulka
,
D.
Wirtitsch
,
V.
Ivády
,
J.
Vodnik
,
J.
Hruby
,
G.
Magchiels
,
E.
Bourgeois
,
A.
Gali
,
M.
Trupke
, and
M.
Nesladek
,
Nat. Commun.
12
,
4421
(
2021
).
14.
M. S.
Barson
,
E.
Krausz
,
N. B.
Manson
, and
M. W.
Doherty
,
Nanophotonics
8
,
1985
(
2019
).
15.
Z.
Su
,
Z.
Ren
,
Y.
Bao
,
X.
Lao
,
J.
Zhang
,
J.
Zhang
,
D.
Zhu
,
Y.
Lu
,
Y.
Hao
, and
S.
Xu
,
J. Mater. Chem. C
7
,
8086
(
2019
).
16.
S.
Baier
,
C. E.
Bradley
,
T.
Middelburg
,
V. V.
Dobrovitski
,
T. H.
Taminiau
, and
R.
Hanson
,
Phys. Rev. Lett.
125
,
193601
(
2020
).
17.
Y.-F.
Gao
,
J.-M.
Lai
,
Y.-J.
Sun
,
X.-L.
Liu
,
C.-N.
Lin
,
P.-H.
Tan
,
C.-X.
Shan
, and
J.
Zhang
,
ACS Photonics
9
,
1605
1613
(
2022
).
18.
C. D.
Clark
and
C. A.
Norris
,
J. Phys. C: Solid State Phys.
4
,
2223
(
1971
).
19.
G.
Davies
and
M. F.
Hamer
,
Proc. R. Soc. A: Math., Phys. Eng. Sci.
348
,
285
(
1976
).
20.
P.
Kehayias
,
M. W.
Doherty
,
D.
English
,
R.
Fischer
,
A.
Jarmola
,
K.
Jensen
,
N.
Leefer
,
P.
Hemmer
,
N. B.
Manson
, and
D.
Budker
,
Phys. Rev. B
88
,
165202
(
2013
).
21.
G.
Davies
,
J. Phys. C: Solid State Phys.
12
,
2551
(
1979
).
22.
N.
Manson
and
J.
Harrison
,
Diamond Relat. Mater.
14
,
1705
(
2005
).
23.
A.
Alkauskas
,
B. B.
Buckley
,
D. D.
Awschalom
, and
C. G. V.
de Walle
,
New J. Phys.
16
,
073026
(
2014
).
24.
G.
Thiering
and
A.
Gali
,
Phys. Rev. B
96
,
081115
(
2017
).
25.
L.
Razinkovas
,
M. W.
Doherty
,
N. B.
Manson
,
C. G.
Van de Walle
, and
A.
Alkauskas
,
Phys. Rev. B
104
,
045303
(
2021
).
26.
See Fig. 5.81. in Ref. 93 where the temperature dependent Debye–Waller (DW) factor is approximately DW 0.095 at T 77 K that translates into S = ln DW 2.35 HR factor. Original (unpublished) work is in Ref. 94.
27.
K.
Wang
,
J. W.
Steeds
,
Z.
Li
, and
Y.
Tian
,
Microsc. Microanal.
22
,
108
(
2016
).
28.
N.
Mizuochi
,
T.
Makino
,
H.
Kato
,
D.
Takeuchi
,
M.
Ogura
,
H.
Okushi
,
M.
Nothaft
,
P.
Neumann
,
A.
Gali
,
F.
Jelezko
,
J.
Wrachtrup
, and
S.
Yamasaki
,
Nat. Photonics
6
,
299
(
2012
).
29.
G.
Liaugaudas
,
G.
Davies
,
K.
Suhling
,
R. U. A.
Khan
, and
D. J. F.
Evans
,
J. Phys.: Condens. Matter
24
,
435503
(
2012
).
30.
D.
Gatto Monticone
,
F.
Quercioli
,
R.
Mercatelli
,
S.
Soria
,
S.
Borini
,
T.
Poli
,
M.
Vannoni
,
E.
Vittone
, and
P.
Olivero
,
Phys. Rev. B
88
,
155201
(
2013
).
31.
E.
Fraczek
,
V. G.
Savitski
,
M.
Dale
,
B. G.
Breeze
,
P.
Diggle
,
M.
Markham
,
A.
Bennett
,
H.
Dhillon
,
M. E.
Newton
, and
A. J.
Kemp
,
Opt. Mater. Express
7
,
2571
(
2017
).
32.
P.
Sáenz de Santa María Modroño
and
G.
Jacopin
, “NV centers in diamond as a CL temperature probe,” poster presented at Hasselt Diamond Workshop 2024—SBDD XXVIII, 2024.
33.
N. B.
Manson
,
M.
Hedges
,
M. S. J.
Barson
,
R.
Ahlefeldt
,
M. W.
Doherty
,
H.
Abe
,
T.
Ohshima
, and
M. J.
Sellars
,
New J. Phys.
20
,
113037
(
2018
).
34.
N. B.
Manson
,
K.
Beha
,
A.
Batalov
,
L. J.
Rogers
,
M. W.
Doherty
,
R.
Bratschitsch
, and
A.
Leitenstorfer
,
Phys. Rev. B
87
,
155209
(
2013
).
35.
A.
Batalov
,
V.
Jacques
,
F.
Kaiser
,
P.
Siyushev
,
P.
Neumann
,
L. J.
Rogers
,
R. L.
McMurtrie
,
N. B.
Manson
,
F.
Jelezko
, and
J.
Wrachtrup
,
Phys. Rev. Lett.
102
,
195506
(
2009
).
36.
L. C.
Bassett
,
F. J.
Heremans
,
D. J.
Christle
,
C. G.
Yale
,
G.
Burkard
,
B. B.
Buckley
, and
D. D.
Awschalom
,
Science
345
,
1333
(
2014
).
37.
Note that 2 λ L ^ z S ^ z spin–orbit coupling operators were defined in Refs. 14 and 16, whereas we define spin–orbit coupling as λ l ^ z s ^ z = λ 2 ( | e + e + | + | e e | | e + e + | | e e | ).
38.
S.
Felton
,
A. M.
Edmonds
,
M. E.
Newton
,
P. M.
Martineau
,
D.
Fisher
, and
D. J.
Twitchen
,
Phys. Rev. B
77
,
081201(R)
(
2008
).
39.
A.
Norambuena
,
A.
Jimenez
,
C.
Becher
, and
J. R.
Maze
,
New J. Phys.
22
,
073068
(
2020
).
40.
A.
Ranjbar
,
M.
Babamoradi
,
M. H.
Saani
,
M. A.
Vesaghi
,
K.
Esfarjani
, and
Y.
Kawazoe
,
Phys. Rev. B
84
,
165212
(
2011
).
41.
W.
Pfäffle
,
D.
Antonov
,
J.
Wrachtrup
, and
G.
Bester
,
Phys. Rev. B
104
,
104105
(
2021
).
42.
L.
Razinkovas
,
M.
Maciaszek
,
F.
Reinhard
,
M. W.
Doherty
, and
A.
Alkauskas
,
Phys. Rev. B
104
,
235301
(
2021
).
43.
G.
Kresse
and
J.
Furthmüller
,
Phys. Rev. B
54
,
11169
(
1996
).
44.
45.
O.
Bengone
,
M.
Alouani
,
P.
Blöchl
, and
J.
Hugel
,
Phys. Rev. B
62
,
16392
(
2000
).
46.
J.
Heyd
,
G. E.
Scuseria
, and
M.
Ernzerhof
,
J. Chem. Phys.
118
,
8207
(
2003
).
47.
A. V.
Krukau
,
O. A.
Vydrov
,
A. F.
Izmaylov
, and
G. E.
Scuseria
,
J. Chem. Phys.
125
,
224106
(
2006
).
48.
P.
Deák
,
B.
Aradi
,
T.
Frauenheim
,
E.
Janzén
, and
A.
Gali
,
Phys. Rev. B
81
,
153203
(
2010
).
49.
A.
Alkauskas
,
M. D.
McCluskey
, and
C. G. V.
de Walle
,
J. Appl. Phys.
119
,
181101
(
2016
).
51.
G.
Thiering
and
A.
Gali
,
Phys. Rev. X
8
,
021063
(
2018
).
52.
G.
Thiering
and
A.
Gali
,
Phys. Rev. Res.
3
,
043052
(
2021
).
53.
K.
Huang
and
A.
Rhys
,
Proc. R. Soc. A: Math., Phys. Eng. Sci.
204
,
406
(
1950
).
54.
See supplementary material in Ref. 95.
55.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
56.
S.
Steiner
,
S.
Khmelevskyi
,
M.
Marsmann
, and
G.
Kresse
,
Phys. Rev. B
93
,
224425
(
2016
).
57.
V.
Ivády
,
I. A.
Abrikosov
, and
A.
Gali
,
npj Comput Mater
4
,
76
(
2018
).
58.
A.
Csóré
and
A.
Gali
,
Phys. Rev. B
102
,
241201(R)
(
2020
).
59.
J. P.
Goss
,
R.
Jones
,
S. J.
Breuer
,
P. R.
Briddon
, and
S.
Öberg
,
Phys. Rev. Lett.
77
,
3041
(
1996
).
60.
M.
Łuszczek
,
R.
Laskowski
, and
P.
Horodecki
,
Phys. B: Condens. Matter
348
,
292
(
2004
).
61.
J. A.
Larsson
and
P.
Delaney
,
Phys. Rev. B
77
,
165201
(
2008
).
62.
C.-K.
Lin
,
Y.-H.
Wang
,
H.-C.
Chang
,
M.
Hayashi
, and
S. H.
Lin
,
J. Chem. Phys.
129
,
124714
(
2008
).
63.
A.
Gali
,
M.
Fyta
, and
E.
Kaxiras
,
Phys. Rev. B
77
,
155206
(
2008
).
64.
P.
Delaney
,
J. C.
Greer
, and
J. A.
Larsson
,
Nano Lett.
10
,
610
(
2010
).
65.
Y.
Ma
,
M.
Rohlfing
, and
A.
Gali
,
Phys. Rev. B
81
,
041204
(
2010
).
66.
S.
Choi
,
M.
Jain
, and
S. G.
Louie
,
Phys. Rev. B
86
,
041202
(
2012
).
67.
J. R.
Maze
,
A.
Gali
,
E.
Togan
,
Y.
Chu
,
A.
Trifonov
,
E.
Kaxiras
, and
M. D.
Lukin
,
New J. Phys.
13
,
025025
(
2011
).
68.
M. W.
Doherty
,
N. B.
Manson
,
P.
Delaney
, and
L. C. L.
Hollenberg
,
New J. Phys.
13
,
025019
(
2011
).
69.
U.
von Barth
,
Phys. Rev. A
20
,
1693
1703
(
1979
).
70.
G.
Thiering
and
A.
Gali
,
Phys. Rev. B
94
,
125202
(
2016
).
71.
I.
Bersuker
,
The Jahn-Teller Effect
(
Cambridge University Press
,
2006
).
72.
I.
Bersuker
and
V.
Polinger
,
Vibronic Interactions in Molecules and Crystals
(
Springer Science & Business Media
,
2012
), Vol. 49.
73.
J.
Zhang
,
C.-Z.
Wang
,
Z.
Zhu
,
Q. H.
Liu
, and
K.-M.
Ho
,
Phys. Rev. B
97
,
165204
(
2018
).
74.
Spherical harmonics are not compatible with the | e ± = [ | e x ± i | e y ] / 2 notation that we use. Instead, there is an extra negative sign: Y 1 ± 1 ( θ , ϕ ) = F ( θ ) e ± i ϕ = F ( θ ) [ cos ( ϕ ) ± i sin ( ϕ ) ] [ | e x ± i | e y ] on which basis the angular momentum ladder operators are defined as l ^ ± | Y 1 0 = 2 | Y 1 ± .
75.
G.
Thiering
and
A.
Gali
, “Nuclear spin relaxation in solid state defect quantum bits via electron-phonon coupling in their optical excited state,” arXiv:2402.19418[quant-ph] (2024).
76.
M.
Mohseni
,
L.
Razinkovas
,
V.
Žalandauskas
,
G.
Thiering
, and
A.
Gali
, “Magneto-optical properties of group-iv-vacancy centers in diamond upon hydrostatic pressure” arXiv: 2408.10407 (2024).
77.
M.
Gajdoš
,
K.
Hummer
,
G.
Kresse
,
J.
Furthmüller
, and
F.
Bechstedt
,
Phys. Rev. B
73
,
045112
(
2006
).
78.
V.
Weisskopf
and
E.
Wigner
,
Z. Phys.
63
,
54
(
1930
).
79.
M. L.
Goldman
,
M. W.
Doherty
,
A.
Sipahigil
,
N. Y.
Yao
,
S. D.
Bennett
,
N. B.
Manson
,
A.
Kubanek
, and
M. D.
Lukin
,
Phys. Rev. B
91
,
165201
(
2015
).
80.
D.
Braukmann
,
E. R.
Glaser
,
T. A.
Kennedy
,
M.
Bayer
, and
J.
Debus
,
Phys. Rev. B
97
,
195448
(
2018
).
81.
M.
Goldman
,
A.
Sipahigil
,
M.
Doherty
,
N.
Yao
,
S.
Bennett
,
M.
Markham
,
D.
Twitchen
,
N.
Manson
,
A.
Kubanek
, and
M.
Lukin
,
Phys. Rev. Lett.
114
,
145502
(
2015
).
82.
G.
Thiering
and
A.
Gali
,
Phys. Rev. B
98
,
085207
(
2018
).
83.
C. M.
Marian
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
187
(
2011
).
84.
K. D.
Jahnke
,
A.
Sipahigil
,
J. M.
Binder
,
M. W.
Doherty
,
M.
Metsch
,
L. J.
Rogers
,
N. B.
Manson
,
M. D.
Lukin
, and
F.
Jelezko
,
New J. Phys.
17
,
043011
(
2015
).
85.
Extrapolating γ in Eq. (5) that of Ref. 84 at T = 0 K and Δ = 47 meV leads to orbital relaxation time in the fs regime.
86.
M.
Bockstedte
,
F.
Schütz
,
T.
Garratt
,
V.
Ivády
, and
A.
Gali
,
npj Quantum Mater.
3
,
31
(
2018
).
87.
E.
Londero
,
E.
Bourgeois
,
M.
Nesladek
, and
A.
Gali
,
Phys. Rev. B
97
,
241202(R)
(
2018
).
88.
P.
Siyushev
,
H.
Pinto
,
M.
Vörös
,
A.
Gali
,
F.
Jelezko
, and
J.
Wrachtrup
,
Phys. Rev. Lett.
110
,
167402
(
2013
).
89.
D.
Wirtitsch
,
G.
Wachter
,
S.
Reisenbauer
,
M.
Gulka
,
V.
Ivády
,
F.
Jelezko
,
A.
Gali
,
M.
Nesladek
, and
M.
Trupke
,
Phys. Rev. Res.
5
,
013014
(
2023
).
90.
N.
Aslam
,
G.
Waldherr
,
P.
Neumann
,
F.
Jelezko
, and
J.
Wrachtrup
,
New J. Phys.
15
,
013064
(
2013
).
91.
R.
Ulbricht
and
Z.-H.
Loh
,
Phys. Rev. B
98
,
094309
(
2018
).
92.
L.
Hacquebard
and
L.
Childress
,
Phys. Rev. A
97
,
063408
(
2018
).
93.
A. M.
Zaitsev
,
Optical Properties of Diamond: A Data Handbook
(
Springer Science & Business Media
,
2013
).
94.
V. S.
Varichenko
, “Luminescence of diamond implanted with high energy ions,” Ph.D. thesis (University of Minsk, 1986).
95.
L. J.
Rogers
,
S.
Armstrong
,
M. J.
Sellars
, and
N. B.
Manson
,
New J. Phys.
10
,
103024
(
2008
).