Ab initio theory of the nitrogen-vacancy center in diamond

1.1 NV center: a brief overview

Point defects may introduce levels in the fundamental bandgap of semiconductors or insulators that radically change the optical and magnetic properties of the host material. In particular, these point defects could be paramagnetic, i.e. the electron spin is greater than zero. A primary example of such a point defect is the NV center in diamond [1]. This defect consists of a nitrogen atom substituting a carbon atom near a missing carbon atom in diamond crystal, i.e. vacancy of diamond. It can accept an electron from the environment and can be negatively charged (see Figure 1A). In this charge state, multiple levels appear in the fundamental bandgap of diamond occupied by four electrons. The corresponding electron configurations constitute an S=1 ground state and an optically active S=1 excited state with below bandgap (5.4 eV) excitation energy [2], [3], [4]. As a consequence, the defect has single-photon absorption and emission spectra in the visible region in the transparent diamond host. Thus, the NV center is also called color center in diamond. The absorption and emission spectra are broad even at cryogenic temperatures caused by the coupling of phonons to the optical transitions. The no-phonon-line or zero-phonon-line (ZPL) optical transition appears at 637 nm (1.945 eV) [5]. In particular, the contribution of the ZPL emission to the total emission, i.e. the Debye-Waller (DW) factor, is about 0.03, with a relatively large Stokes shift of about 0.45 eV. This makes possible to excite the NV center by green light (typically 532 nm excitation wavelength) and to detect the emitted photons in the near-infrared region of 700 to 900 nm. The S=1 spin in the ground state can be measured by conventional electron spin resonance (ESR) techniques [6]. It is important to note that this defect always is a two-spin system as the nitrogen isotopes either have I=1 (¹⁴N) or I=1/2 (¹⁵N) nuclear spin. As a consequence, the hyperfine interaction between the electron spin and the nitrogen nuclear spin always occurs and nuclear quadrupole coupling appears for the largely abundant ¹⁴N isotope. In addition, the hyperfine signature of ¹³C I=1/2 nuclear spins is also observable (natural abundance of 1.1%) in the corresponding ESR spectrum [3], [7], [8]. The spin Hamiltonian of the NV center in the ground state can be written as

Figure 1:

(A) NV defect in diamond. Vacancy is depicted as a circle in the middle of the diamond cage (gray balls are carbon atoms). An electron is provided by the diamond environment, e.g. N-donor, which creates the negatively charged NV defect, i.e. NV(–), briefly NV center. Three carbon atoms and the nitrogen atom have dangling bonds pointing toward the vacancies. The defect has a threefold rotation axis (C₃) about the ⟨111⟩ axis of diamond. The eigenstates of the NV center can be labeled by C_3v symmetry irreducible representations. (B) Electronic structure and decay processes. The double group representations for the triplet states are also depicted with their corresponding spin eigenstates. Red arrows represent emission (weak between singlets). Curved arrows show the ISC processes that are mediated by phonons and the parallel or perpendicular components of spin-orbit interaction (H^SO), where the dominant components are shown with solid lines. Dark gray arrows show the absorption that may lead to the ionization of the NV center. Note that hyperfine interaction occurs between the electron spin and the nitrogen nuclear spin (not depicted).

where the gyromagnetic factors of the electron spin were first observed as isotropic g_xx=g_yy=g_zz=2.0028±0.0003 (Refs. [7], [9]), whereas slight anisotropy, g_xx=g_yy=2.0029(2) and g_zz=2.0031(2), was reported in Ref. [8], which results in the corresponding Zeeman splitting upon external magnetic field B. D is the so-called zero-field splitting (ZFS) tensor, where “zero field” refers to zero magnetic field. In this particular system, the eigenstates of this part of the Hamiltonian can be given as D(Sz2−2/3), where S_z={−1; 0; +1} is the eigenvalue of the electron spin and the D ZFS constant is about 2.87 GHz [7]. This term predominantly arises due to the dipolar electron-spin-electron-spin interaction and results in a D energy gap between ms=0 and ±1 levels in the ground state at zero magnetic field when the quadrupole and hyperfine terms are neglected. C_Q=−4.945 MHz (Ref. [10]) and A_zz,N=−2.14±0.07 MHz; A_xx,N=A_yy,N=A_⊥,N=−2.70±0.07 MHz are the quadrupole strength and hyperfine principal values of ¹⁴N [8], respectively, which are always present for each individual NV center. A hyperfine tensor is nearly isotropic for the given nitrogen isotope but can be highly anisotropic for the proximate ¹³C nuclear spins [3], [7], [8], [11], which generally reduce the symmetry of the spin Hamiltonian. For distant ¹³C nuclear spins, the simple point dipole-point dipole approximation for the hyperfine interaction between the electron spin and nuclear spins is valid; however, the Fermi contact term, i.e. the localization of the electron spin density at the site of nuclear spins, can be significant for the proximate ¹³C nuclear spins [3], [8], [11]. When the magnetic field approaches the ZFS, then the nuclear Zeeman terms of the nitrogen nuclear spin and ¹³C spins with the g_N and g_C nuclear gyromagnetic ratios, respectively, become an important effect [12].

A key property of the NV center is that the ESR is correlated with the spin-selective fluorescence intensity [13] and two-photon ionization probability [14], where the latter occurs by subsequent absorption of two photons via the ³E excited state. If the concentration of NV centers is low in diamond and the photoexcitation source is focused into a single NV center, then the ESR of the single NV center can be detected either optically [13] or electrically [15] that are called optically detected magnetic resonance (ODMR) and electrically detected magnetic resonance (EDMR), respectively. As the latter was achieved by measuring the photocurrent, it is also called photocurrent-based detected magnetic resonance (PDMR). In ODMR measurements, it was found that the fluorescence intensity between the ms=0 levels of the excited and ground states is about 30% stronger than that between the ms=±1 levels of the excited and ground states. This is also called the ODMR readout contrast. In other words, this readout process converts the ESR frequency used in electron paramagnetic resonance techniques in the microelectronvolt region into the frequency of optical photons in the electronvolt region. In addition, the number of optical photons emitted in the ODMR process is several orders of magnitude larger per center than the single microwave photon absorbed in the ESR process, which increases the detection sensitivity by a large amount. Another consequence of the ODMR process is that the optical pumping of the NV center results in almost 100% population of the ms=0 level in the ground state [6], [7], [16]. It has been recently found that a similar PDMR readout contrast is achievable for single NV centers [15]. In the PDMR readout, the readout process ends at the neutral NV defect [17]. Further optical pumping turns the neutral NV defect into the negatively charged NV defect, i.e. the NV center [18], [19], [20].

The readout and electron spinpolarization processes inherently contain flipping of the electron spin during the decay from the excited state to the ground state, which is highly selective to the ms=±1 spin state of the triplet excited state. Selection rules in the decay processes can be efficiently analyzed by group theory. The NV center in diamond exhibits C_3v symmetry with a symmetry axis lying along the ⟨111⟩ directions of diamond [5]. By combining the defect-molecule diagram and group theory [21], [22], [23], one can predict the character of the triplet ³A₂ ground and ³E excited states as well as the dark singlet states also as a function of external perturbations [24], [25]. Spin-orbit interaction may connect the triplet states and singlet states resulting in nonradiative decay and spin flipping. This type of nonradiative decay is called intersystem crossing (ISC). Understanding ISC is the key for ODMR and PDMR readout.

The connection between quantum information science and NV centers is intimately bound to the ODMR and PDMR processes and readout. These readout processes make it possible to coherently manipulate single electron spins in solids by microwave fields and optical excitation [16], [26], where two states of the electron spin of the NV center realize a quantum bit, which can be read out and initialized by optical means. These readout mechanisms operate at room temperature [13], [15] and even at higher temperatures [27], and the measurement scheme can be pushed up to 1000°C with pulsed protocols [28]. The coherence time of the NV center’s electron spin in diamond with natural abundant ¹³C isotopes can reach ≈600 μs [16], [29] and up to 2 ms in ¹²C-enriched diamonds [30] as obtained by Hahn-echo measurements, even at room temperature. By control-NOT operation [31], [32], the quantum bit information can be written from the electron spin to the nuclear spin. Single-shot readout of the quantum bit was demonstrated [33] as well as the coding of quantum information from the electron spin [34] or from the nuclear spin [35] to the polarization of the emitted photon or the photon emission itself at a given frequency [33], which realize spin-to-photon interfaces at cryogenic temperature. This was used to transmit quantum information over 1 km distance [36]. One can conclude that the NV center in diamond satisfies DiVincenzo’s criteria of quantum information processing [37]. In the following, we list these criteria that we amend with the NV properties: (1) scalable physical system with well-characterized quantum bits, i.e. the electron and nuclear spins of a single NV center; (2) ability to initialize the state of the quantum bits to a simple fiducial state, such as initialization of the ms=0 state by optical means; (3) long relevant decoherence times, much longer than the gate operation time, i.e. millisecond coherence time of the NV center; (4) a “universal” set of quantum gates, e.g. control-NOT operation, as demonstrated for the NV center; (5) a quantum bit-specific measurement capability, i.e. ODMR or PDMR readout; (6) ability to interconvert stationary and flying qubits, such as spin-to-photon interface of the NV center; and (7) ability to faithfully transmit flying qubits between specified locations as demonstrated between remote NV centers. Scalability is still an issue because identical NV centers require Stark shift tuning of the levels in the excited state [38], [39], [40], and despite the technological efforts on creating arrays or clusters of NV centers in diamond [41], [42], [43], [44], [45], a robust coherent coupling mechanism between multiple NV centers has not yet been demonstrated for more than three NV centers [45].

The NV center in diamond has a favorable electron spin coherence time, but this coherence time and the (spin) levels are relatively sensitive to the environment. This is undesirable for building up a quantum computer from the NV centers but can be harnessed in the measurement of magnetic [29], [46], [47], electric [48], and strain [49], [50], [51], [52] fields and temperature [53], [54], [55] at the nanoscale. The room temperature operation and favorable coherence time of the NV center paves the way toward nuclear magnetic resonance (NMR) of single molecules at ambient conditions [56], [57], [58], [59], [60], [61], [62], [63], [64]. Furthermore, the possibility of spinpolarization transfer from the electron spin toward the nuclear spins [11], [65], [66], [67] can be employed to hyperpolarize diamond particles [68], [69], [70], [71], [72], [73], [74] or external species [75], [76] attached to the diamond surface. In these quantum sensing and related applications, the NV center should reside close to the surface of diamond.

Theory proposed [77], [78] that near-surface NV centers can be also used as a resource to carry out quantum simulations on frustrated quantum magnetism or quantum spin liquid if the surface of diamond is decorated by species with nuclear spins, e.g. I=1/2 fluorine [77] or I=1 ¹⁴N [78], respectively.

NV centers can be observed in natural diamonds, but controlled preparation is required for the aforementioned quantum technology applications. The intentional production of NV defects often starts with a high-quality diamond that was grown by chemical vapor deposition (CVD) with minimal contamination. In the next stage, nitrogen ions are implanted into diamond. After implantation, the diamond sample is annealed to remove the damage created by implantation and facilitates the formation of the NV center in diamond [79], [80], [81], [82], [83]. The negative charge state of the NV defect is most likely provided by substitutional nitrogen (N_s) donor defects that are called P1 center named after its ESR signature [84], [85], [86]. The shallow NV centers near the surface of diamond are dominantly created by this implantation technique [82], where the depth of the defects can be controlled by the energy of the bombarding ions. It was found that the coherence time and often the photostability of these shallow NV centers are compromised [82], which might be either related to the quality of the diamond surface or the quality of diamond crystal around the shallow NV centers. Indeed, the nitrogen delta-doping technique of diamond growth for creating shallow NV centers improved the coherence time [87].

In some sensing applications, large ensembles of NV centers are needed. In that case, the starting material could be a heavily nitrogen-doped diamond and then vacancies can be formed by different irradiation techniques (such as electron [88], [89] or neutron irradiation [17] or implantation with inert ions [90], [91]) in that sample. Finally, annealing is applied to remove the irradiation damages and facilitate the formation of NV defects [92]. It was found that the optimization of the annealing stages is very important to achieve good coherence properties of the resulting NV centers [93].

This brief overview could lead to the impression that the NV center in diamond is probably the most studied point defect in the experiments. Still not all the signatures are well interpreted and understood purely from experimental spectra and simple models. The ab initio theory could significantly contribute to understanding the formation, photoexcitation, and photoionization processes as well as the ISC processes of NV centers. In particular, results from ab initio calculations could guide the idea of improving the coherence properties of NV centers after nitrogen implantation [94] and the notion [14], [20] and optimization [95], [96], [97] of PDMR readout of the NV center in diamond. We provide a recipe below how to achieve the full ab initio description of a solid-state defect quantum bit on the exemplary NV center in diamond.

1.2 Ab initio description of the NV center in diamond: a desideratum

Two major goals can be identified in the ab initio description of the NV center in diamond: (i) creation of the NV center and its interaction with the other point defects and diamond surfaces and (ii) determination of ionization energies and magneto-optical properties with the corresponding radiative and nonradiative rates, also as a function of external perturbations such as magnetic and electric fields, strain, and temperature.

Issue (i) is a very common target in point defect studies, where the formation energy or enthalpy of the point defects should be determined and the surface of the host should be studied as a function of the environment, which may influence the surface morphology and termination of the host material. However, issue (ii) requires in-depth investigation of the behavior of a point defect, much deeper than usual in the community of researchers working on point defects in solids, as the thermodynamic properties of point defects are often the main target in the vast majority of such investigations. However, vacancy formation is not necessarily a quasiequilibrium process. ODMR and PDMR signals of the NV center arise from photoexcitation, which is out of thermal equilibrium of the electrons. In addition, the magnetic properties play a central role in solid-state defect quantum bits that should be described in detail in both the ground and excited states.

The list of properties for the full ab initio description of the NV center can be recognized in Figure 1B. In the ground state, the ZFS due to dipolar electron spin-electron spin interaction should be computed. The hyperfine interaction of the electron spin with the nuclear spins is also very important in the entanglement schemes [16], quantum memory [35], [98], [99], quantum error correction [100], and hyperpolarization [65] schemes as well as in understanding the decoherence of the NV center’s electron spin. For instance, by controlling the proximate nuclear spin states around the NV center, the coherence time of the electron spin has been pushed beyond 1 s [101]. The calculation of the nitrogen quadrupole constant can be important in understanding the so-called “dark state” NMR of ¹⁴N [102], which can be used to identify the charge state of the NV defect [10].

In the ³E excited state, the spin levels are heavily temperature dependent and show complex features as observed in the photoluminescence (PL) excitation (PLE) spectrum [103]. The ZFS occurs due to the dipolar spin-spin interaction, but also spin-orbit interaction takes place. As a consequence, the optically active spin triplet state should be calculated with the corresponding electron dipolar spin-spin interaction and spin-orbit interaction.

In the ODMR readout process, the calculation of ISC rate requires the spin-orbit matrix element between the triplet states and the corresponding singlet states (see Figure 1B). This assumes that the optically inactive or dark singlet excited states and levels should be calculated too.

In the PDMR readout process, the probability of photoionization either directly from the ground state or via the real excited state is a key issue. Photoionization could be direct or via Auger process; thus, the photoionization rate of both processes should be calculated. These nonradiative rates compete with the radiative rate from the ³E excited state to the ³A₂ ground state, so the latter should be also determined for the sake of a complete description.

Phonons are involved in all of these processes; thus, electron-phonon coupling should be calculated to understand the DW factor of the PL spectrum and the ISC processes between the triplet and singlet states [104], [105].

One should recognize that all the interactions between electron orbitals, electron and nuclear spins, phonons, and external fields should be considered and calculated ab initio for a complete description of the operation of solid-state defect quantum bits in a realistic solid and environment. At least, two main challenges can be identified: (a) calculation of the excited states for sufficiently large models of the NV center and (b) treatment of the electron-phonon coupling in the radiative and nonradiative processes. In the following, recent efforts along these directions are presented together with the other developments and results in the field.

1.3 Contents of the review paper

The rest of the paper is structured as follows. The basic ab initio methods for studying solid-state defect quantum bits are presented in Section 2. We then summarize the recent methodology developments for the description of the NV center in diamond and the corresponding results in Section 3. This section starts with the basic modeling and ground-state properties in thermal equilibrium and then continuous with the treatment of excited states and optical properties, which includes the discussion of the participation of phonons in the optical transition. Next, the computational methods of the magnetic parameters are discussed, which would complete the description of the magneto-optical properties of the NV center in diamond. This is the starting point to determine the radiative and nonradiation rates as described in the next section. Finally, the simulation tools for calculating the various sources of perturbation on the aforementioned properties are presented. An outlook is provided for the next steps toward the full ab initio description of the NV center in diamond in Section 4. Finally, the paper is summarized in Section 5.

3.1 Formation energies and charge transition levels

The formation energy of the defect (Eformq) in the charge state q can be calculated from the total energy of the defect (Etotq) and the corresponding chemical potential of the atoms (μ) constituting the defect and the electron [123], which can be written for the 512-atom diamond supercell with the NV center as

where the chemical potential of the carbon atom can be expressed from the total energy of the perfect diamond lattice [E_tot(C512)], whereas the chemical potential of the nitrogen atom (μ_N) depends on the growth conditions. It can be set to the half of the total energy of the nitrogen molecule as a reference [124]. E_F is the Fermi energy referenced to the valence band top E_V. Ecorrq is the correction of the total energy for charged supercells, because the charged supercells are neutralized by a jellium background in the plane-wave supercell calculations that can all interact with their periodic images [125] with a resulting shift in the total energy. For the medium localization of defect states, such as the NV center in diamond, Lany-Zunger correction [126] and Freysoldt correction [127] yield equivalent and relatively accurate results in sufficiently large supercells [128]. The formation energy of the defect provides at least two very important quantities: (i) the concentration of the defect (N_i) in thermal equilibrium can be calculated as

where Ni0 is the density of i sites in the perfect lattice, k_B is the Boltzmann constant, and T is the temperature in Kelvin, and (ii) the adiabatic ionization energies or occupation level of the defect between q and q+1 charge states can be calculated as

which gives the position of the Fermi level with respect to E_V, where the concentrations of the defects in charge states q and q+1 are equal, E(q|q+1)=Etotq+1+Ecorrq+1−Etotq−Ecorrq. Although the controlled formation of the NV center is not a thermal equilibrium process, still formation energies and the corresponding concentration of defects can provide information about their abundance, in particular, when the occurrence of various defects in diamond is studied and compared to each other. Regarding the adiabatic ionization energies, it is highly important to apply such a functional that is able to reproduce the experimental bandgap; otherwise, the calculated ionization energies with respect to the band edges are not comparable to the experimental data. As explained above, LDA or PBE DFT predicts about 1.2 eV lower bandgap than the experimental one. In contrast, it has been found by applying intensive tests on defects group IV semiconductors and diamond [129], the so-called HSE06 screened and range separated hybrid functional [130], i.e. the mixture of screened Fock exchange and PBE exchange functionals, are able to reproduce the experimental bandgaps and the ionization energies are also reproduced within 0.1 eV accuracy. The test provided good results also for the NV center in diamond [124], [129]. This is a huge improvement over the accuracy of LDA and PBE functionals. The calculated HSE06 formation energy of the NV center in diamond is illustrated in Figure 3. These calculations were carried out for N_s. Carbon vacancy, nitrogen di-interstitial (N₂), N₂V (two nitrogen atoms substitute carbon atom near an adjacent vacancy), divacancy (V₂), and NVH defect, where a hydrogen atom saturates one of the carbon dangling bonds in the NV defect. By calculating the formation energy of these defects and applying the charge neutrality condition, the concentrations of defects were determined under quasi-thermal equilibrium conditions at a given temperature that might occur in the CVD process [124]. It was found that the concentration of the NV center would be very low at any nitrogen concentrations because of the favorable formation of either N₂ or NVH. NV centers can be rather created in nitrogen-doped diamond after irradiation, which creates carbon vacancies (V). The formation energy of the defects can be used to calculate the energy balance of defect reactions, i.e. to study defect chemistry in diamond. As an example, at relatively low Fermi-level positions (such as at E_V+2.0 eV, where V is neutral and mobile), the combination of two vacancies results in V₂ with energy gain of 4.2 eV, whereas the energy gain of combining N_s and V to form NV is only 3.3 eV. This means that rather divacancy forms than the NV center when single vacancies diffuse in the diamond crystal upon annealing. However, it is much likely (by ~5 eV) to remove a carbon atom near N_s than that from perfect lattice. These results implied that the formation of NV occurs, when the vacancy is formed very close to N_s. Divacancies are paramagnetic, electrically and optically active defects (see Ref. [124] and references therein), which can be detrimental for the charge state stability and photostability of the NV center [124]. It was proposed that high-temperature annealing at about 1200°C to 1300°C may reduce the concentration of vacancy clusters that can improve the properties of NV centers [124].

Figure 3:

Formation energy of NV center in diamond [124].

The chemical potential of nitrogen was set to the energy of nitrogen atom in the nitrogen molecule at T=0 K. The crossing lines correspond to the adiabatic ionization energies. The Fermi level is aligned to the top of the valence band. We note that these data differ, in particular, for the donor level from Ref. [131], which was obtained in a small 64-atom supercell.

This result could motivate a recent experimental effort to hinder the formation of divacancy or vacancy clusters by engineering the charge state of the single vacancies so that they are not mobile anymore [94]. In particular, in boron-doped diamond, the single vacancies become positively charged or even double positively charged [124], and the formation of vacancy clusters should be significantly reduced. To this end, Fávaro de Oliveira et al. produced a thin boron-doped layer near the region where nitrogen ions were implanted [94]. Indeed, this method resulted in an improved charge stability and coherence properties of NV centers in diamond [94]. We note here that atomistic molecular dynamics simulations using empirical potentials showed in this study [94] that indeed divacancies and larger vacancy aggregates are formed near the NV center after N ion implantation when the vacancies are neutral, and the proximate vacancy aggregates’ spin causes the decoherence of the NV center’s electron spin [94].

3.2 Diamond surface

The shallow implanted NV centers below the diamond surface may suffer from the vacancy cluster formation as explained above, (e.g. Ref. [132]) but the presence of the diamond surface can also pose a problem. Typically, NV centers were implanted into the (100) diamond surface. It was found that shallow NV centers in hydrogen-terminated diamond are not stable. This can be explained by the fact that hydrogenated (100) diamond has an electron affinity of about −1.3 eV [133] as predicted by Shockley [116]. When water absorbs to this diamond surface, it creates a band bending of diamond [134], which shifts down the Fermi level and converts the NV center to neutral NV defect [135], [136], [137]. As a consequence, the stability of NV centers requires diamond surface with positive electron affinity. An additional criterion is that no surface-related level should appear in the fundamental bandgap of diamond. Ab initio DFT supercell calculations can predict the electron affinity of surfaces of solids by comparing the calculated conduction band minimum (E_C) with respect to the vacuum level. In the supercell modeling of surfaces, a slab model is the only option, e.g. diamond slab with (100) bottom and top surfaces [138]. It is often desirable to choose the same type of termination at the bottom and top of the slab to avoid an artificial polarization across the slab, and the size of the vacuum region should be sufficiently large for convergent electron affinity calculations (see Figure 4A). Again, the choice of the functional is crucial. In early calculations, LDA or PBE DFT calculations were applied with too low bandgap [140]. To obtain realistic results, the calculated E_V was fixed and E_C was shifted to reproduce the experimental bandgap, i.e. scissor correction was applied to E_C [140]. However, HSE06 functional calculations with accurate band structure calculation showed [120] that E_V should shift down and E_C should shift up by about the same amount in the correction of PBE band edges. Thus, the trends could well produce the LDA calculations in Ref. [140], but accurate results are expected from HSE06 electronic structure calculations [120].

Figure 4:

(A) Schematic diagram of the slab modeling of point defects at the surface: NV is an acceptor, whereas N is a donor. The entire slab is neutral. (B) Proposed smooth oxygenated (100) diamond surface for NV sensor applications. (C) Nitrogen-terminated (100) diamond surface. (D) Nitrogen-terminated (111) diamond surface. (E) Mixed oxygen-terminated (113) diamond surface with an epoxy C–O–C bonding configurations shown also from the side view with typical bondlengths in Ångström unit. Data were taken from Refs. [78], [120], [139].

In the experiments, oxygenation is applied to stabilize the charge state of shallow NV centers in (100) diamond. It was found that alcohol groups (–OH) would not turn the negative electron affinity of hydrogenated diamond surface to positive electron affinity. Rather, ether bridges (C–O–C bonds) are responsible for this issue [120]. However, ether bridges will stiffen the diamond lattice with introducing defect levels in the gap of diamond, and these groups most likely roughen the diamond surface with creating optically and electrically active defects. Kaviani et al. suggested to apply smooth and well-controlled oxygenation techniques that would result in a mixture of C–H, C–OH, and C–O–C bonds at the diamond surface with no stress at the diamond surface but with positive electron affinity and clean bandgap (see Figure 4B). It was also found [120], in agreement with Ref. [140], that smooth fluorine termination provides a large positive electron affinity (100) diamond surface, although an empty band appears just below E_C according to HSE06 DFT [120]. Motivated by experimental efforts [141], nitrogen-terminated diamond surface was also studied from first principles, and positive electron affinity was predicted for larger than 0.5 monolayer of N–N groups rather than C–H groups at the surface (see Figure 4C). Indeed, a recent study has found improved NV properties on nitrogen-terminated diamond (100) surface over oxygenated one, although the improvement was not striking [142]. Most likely, the nitrogen termination of the (100) diamond surface was not perfect in that experiment. HSE06 DFT calculations predicted [78] that nitrogen-terminated (111) diamond surface is more favorable than that on (100) diamond surface, as nitrogen can naturally replace the top C–H layer on (111) diamond without introducing any strain [78]. Thus, it is likely that (111) nitrogen-terminated diamond surface may host NV centers with excellent properties for sensing applications (Figure 4D). This is appealing as preferential alignment of NV centers along the (111) axis can be realized [143], [144], [145], [146], up to 99% [147], [148], in the growth process of diamond, which is desirable for magnetometry and related applications. However, the growth of high-quality (111) diamond at sufficiently high rate is still a big challenge. Alternatively, (113) diamond can be grown with 79% preferential alignment of NV centers at a considerable growth rate [149], [150]. HSE06 DFT calculations have recently obtained a surprising result that nitrogen termination is not preferred for this surface because of introducing surface-related bands into the bandgap of diamond, but rather oxygen termination may result in an excellent environment for hosting NV centers for quantum sensor applications [139]. On the (113) diamond surface, oxygen can form so-called epoxy bonds with the surface carbon atoms that are stable according to the ab initio simulations (Figure 4E). This bonding situation does not frustrate the top carbon layers; thus, it produces clean bandgap and positive electron affinity.

The calculated electron affinities in seeking positive electron affinity are presented in Table 1. We note that we do not provide here a comprehensive list of surface terminators creating negative electron affinities that are not relevant for NV sensor applications.

Previously, we mentioned above that the preferential alignment of NV centers have been reported during the in-growth process of NV centers in the CVD of the (111) and (113) diamond. We note here that a single DFT study predicted [151] that the preferential alignment of NV centers might be achieved in bulk diamond by inducing 2% biaxial strain during the diamond growth at about 970°C temperature; however, this has been not verified in the experiments. The preferential alignment of NV centers in (110) and (111) diamond surfaces was studied by ab initio calculations [152], [153]. It was found that nitrogen prefers to reside on the carbon layer beneath the top carbon layer of the surface [152], and the NV centers will be preferentially aligned at the kink of the diamond terrace during the growth of the (111) diamond under hydrogen-rich environments [153].

Realistic diamond surfaces may contain steps, voids, and other defects that were not considered in previous modeling. An exception is the carboxyl group (double carbon-oxygen bonds) at the oxygenated (100) diamond surface, which creates levels in the bandgap diamond with localized states that can absorb light in the visible and perturbs NV measurements [120]. Thus, carboxyl groups should be eliminated from the oxygenated diamond surface. Recently, very surprising results have been reported on the (100) diamond surface that could have a direct relation to surface-related charge or spin noise [154]. Ab initio modeling tentatively identified stable pairs of sp² C–C bonds in the void defect, i.e. single surface carbon vacancy at the top of diamond surface, which introduces an empty level about 1.4 eV above E_V as found by surface-sensitive analysis techniques [154]. This empty state may be filled by electrons during illumination of the diamond surface that provides a surface charge and spin. This defect is very general and relatively abundant on the (100) diamond surface that might be the source of surface noises felt by near-surface NV centers [154]. The interaction of the NV center with a nearby acceptor defect has been recently modeled in a bulk supercell by DFT calculations [155]. The bulk model was chosen to avoid modeling problems of the slab calculations. The acceptor defect was chosen to be the neutral NV defect for simplicity [155]. A common sense among experimentalists is that the NV center is an atomic-like defect with very localized orbitals. Chou et al. showed in the study of interaction of the NV center with neutral NV acceptor [155] that the NV orbitals spread to 4 nm from the center of the NV defect in the plane of the three-carbon atoms nearest to the vacancy. If the acceptor defect is closer than about 7 nm distance from the NV center with similar wave function extension, then they can directly interact quantum mechanically without any illumination of the NV center, which can lead to a decrease in the spin coherence time of the NV center [155]. The calculated quantum mechanical tunneling rate between the defects could well explain the experimental data in diamond sample with high density of NV defects [156], which highlights the predictive power of DFT methods. Recently, it has been found in experiments [157] that shallow NV centers can ionize in the dark on experimentally relevant timescales, which can be understood as tunneling to a single local electron trap as the mechanism behind this process.

NV centers were also considered close to diamond surfaces using slab models. The modeling of the NV center in slab models faces several problems: (i) the charge correction of the negatively charged defect is not straightforward as the potential of a point charge in a surface is only conditionally convergent; thus, the handling of charged slab supercell is painstaking. (ii) Artificial polarization might appear due to the bottom of the slab (double surface). Lastly, (iii) the finite width of the slab may introduce a quantum confinement effect. Basically, the total energy of the NV center in the slab model may be converged by increasing the lateral size of the slab model [138] and by applying total energy correction techniques that were recently proposed in the literature [128], [158], [159], [160], [161], [162], [163], [164], [165], [166]. The sufficiently large lateral size is necessary to minimize the artificial interaction of the defect with its periodic images too [167] as usual in 3D bulk modeling. The other two problems may be minimized by adding the same surface termination at the bottom and top of the slab and by using a sufficiently large width of the slab [120]. According to intensive tests, about 2.2 nm width produces minute quantum confinement effect [120].

Kaviani et al. invented to do a workaround in the problem of charged slab supercell by replacing it with another but readily solvable problem [120]: a neutral slab model is used for the negatively charged NV defect (NV center) at the expense that another defect enters the slab. Namely, the N_s donor will donate an electron to the neutral NV acceptor defect by creating a pair of the NV center and positively charged N_s. If these defects are placed into the same layer of the slab and the slab has a cubic-like shape, then the dipole-dipole interaction between the periodic images of the defect pairs can be minimized. Rather, the presence of ionized N_s near the NV center may shift the levels or split the degenerate levels of the NV center. Extensive tests showed that if the defects are placed at least 7.5 Å far from each other, then the degeneracy of the corresponding orbitals is maintained and the constant shift in the KS levels can be well monitored and corrected [120]. In a recent study, the pair of the NV center and ionized N_s was analyzed in detail by arriving at the same conclusion that 7.5 Å distance between the two defects suffices to avoid the splitting of the degeneracy of the corresponding orbitals, and the NV center can be well approximated as isolated [168]. These achievements make possible to directly simulate the NV center with correct total energies that is necessary for modeling direct interaction of the NV center with surface species [120].

We mention in this context that direct tunneling between the nitrogen donor and neutral NV defect has been recently discussed in detail in type 1b diamond [169], which study largely extends the original idea from Collins [170] about the electron transfer between these nearby defects. We note that one conclusion from the ab initio study of the interaction of the NV center and neutral NV acceptor defect that the rate of tunneling depends on the actual orientations of the nearby defects because of the direction of the spatial extension of the orbitals [155]. This should hold for the nitrogen donor-NV defect pair too as the neutral nitrogen donor orbitals shows special spatial extension due to the giant Jahn-Teller distortion (e.g. Refs. [124], [171]); thus, accurate determination of the tunneling rate as a function of distance and relative orientation between the two defects requires ab initio calculations.

3.3 Excited states

3.3.1 Electronic solution

The quantum bit operation and readout works via optical excitation; thus, understanding the absorption and decay from the excited state is highly important. To this end, the first task is to calculate the excited states and levels. This is far from trivial. The extension of wave function of the NV center in diamond requires sufficiently large models for accurate calculations [105], [155]. In addition, if diamond band edges are also involved in the photoexcitation process, then the diamond host should be modeled at equal footing with the NV center itself. As a consequence, the highly accurate but extremely costly wave function-based quantum chemistry methods have only limited accuracy, as these methods can be applied on small molecular cluster models [108], [109]. KS DFT, including HSE06 DFT, which can be used in large supercells in practice, has built-in limitations because of the approximations in the KS DFT functionals; in particular, the highly correlated singlet excited states, e.g. the ¹A₁ state, cannot be directly calculated by KS DFT correctly. Highly correlated states may be recognized in the combination of defect-molecule picture and group theory as multideterminant states [2], [23], [24], [25]. Defect states in solids may be viewed as highly correlated electron states in contact with a bath of extended states. In a recent work, the wave function method is embedded into the DFT framework to synthesize the advantages of both methods and follow the character of the system [172]. For the chosen orbitals (typically defect orbitals), the electron-electron interaction is exactly calculated as Coulomb interaction between the electrons, i.e. configurational interaction (CI), whereas the interaction of the other electrons is treated with HSE06 DFT. The interaction between the chosen orbitals and the extended states is treated within random-phase approximation in a way that the corresponding dielectric function calculated with the constraint of extended states (cRPA) [172]. This CI-cRPA method does not contain any fitting parameters, and it is transferable in the sense that the number of orbitals in the CI active space can be systematically increased. This method can be used to analyze the character of the wave functions and the features of the other methods. In particular, the position of the ¹A₁ level is very sensitive to the screening of the electrons from the diamond bands [172]. As a consequence, pure quantum chemistry CI method in tiny diamond cluster models yield too high ¹A₁ level [109]. In a diamond cluster Hubbard model calculation, it was shown that the ¹A₁ state inherits double excitation from the lower a₁ level to the upper e level in the gap [173], which provides an important insight about the nature of the ¹A₁ state, although the limitation of the diamond cluster model resulted in a false position of the ¹A₁ level in this calculation [172]. For this reason, the very popular many-body perturbation method on top of DFT calculations, called GW+BSE (see Ref. [174] and references therein), fails to properly describe ¹A₁ level [175], [176] because BSE can typically describe excited states with the combination of single excitation configurations. The extended Hubbard model calculation fit to the GW calculation of the in-gap defect levels provides accurate singlet and triplet levels of the NV center [172], [176]. We note here that the NV center is exceptional in the sense that the accurate description of the electronic structure of other solid-state quantum bits requires to consider resonance orbitals from the valence band in the CI active space, for which the fitting procedure is ill defined [172]. The final conclusion is that the ¹A₁ level resides about 0.4 eV below the excited state ³E level. We note that the electronic ¹E state is basically stable against distortion (it can be described as a single Slater determinant in an appropriate basis of e_x,y orbitals in Figure 5A [4], [177]), and it becomes slightly unstable against distortion because of the appearance of ¹E′ character in the ¹E state [172], [177]. The predominant simple character of the ¹E state makes the ¹E level almost insensitive to the choice of the computational method that lies about 0.4 eV above the ground-state ³A₂ level. The final level diagram is shown in Figure 5 in which the global energy minimum of the corresponding APES is taken into account (see below). The final picture well accounts for the measured ZPL energy between the singlet states of about 1.19 eV [178], and derivation from the measured nonradiative rates from the ³E level to the ¹A₁ combined with a phenomenological model on the phonon participation and density of states [104] that concluded about 0.4 eV gap between these two levels.

Figure 5:

(A) Single-particle scheme of the electronic structure of the NV center. The spinpolarization between spin-up and spin-down electrons in the KS DFT results in different levels for spin-up and spin-down orbitals. VB and CB refer to valence band and conduction band, respectively. The fundamental bandgap of diamond is 5.4 eV. Green arrow represents the optical transition within a single-particle scheme. (B) Level structure of the NV center with taking the geometry relaxation in the corresponding electronic states into account. The character of the wave functions are depicted. Tilde labels the vibronic or polaronic nature of the state, which is particularly strong for the ³E and ¹E states.

The CI-cRPA method also justified the so-called constraint occupation DFT or ΔSCF method for calculating the ³E excited state [2], [179] of which level is relatively insensitive to the size of the diamond cluster and different computational methods. In the constraint occupation DFT, the electronic configuration is set to build up the Slater determinant of the ³E excited state as shown by the green arrow in Figure 5A. A big advantage of the ΔSCF method is that the forces acting on the ions can be straightforwardly calculated; thus, the APES of the ³E state can be mapped by the KS DFT method. This is of high importance in understanding the nature of the interaction of electrons and phonons in the excited state and its role in the optical spectrum and nonradiative decay.

3.4 Strong electron-phonon coupling: vibronic wave functions

Experiments indicated [180], [181] that the dynamic Jahn-Teller (DJT) effect occurs in the ³E excited state, which is double degenerate and has two components: |E±〉=12(|Ex〉±i|Ey〉). Symmetry breaking e_x,y-type phonons distort the system and couple strongly to the ³E electronic state [105], [182]. This is a so-called E⊗e DJT system for which the APES shows a sombrero shape (two shifted parabolas with minima at distorted C_1h geometries, which crosses at a conical intersection point at the high C_3v symmetry as shown in Figure 6) rather than a simple parabola with the minimum at the high symmetry point [105], [183]. By calculating the APES of the ³E state and applying the E⊗e DJT theory from Bersurker [184], one can setup an electron-phonon Hamiltonian with a single effective phonon mode for which the electron-phonon parameters can be derived from the calculated APES [105], [183]. The limitation of the single effective phonon mode is briefly discussed in Outlook Section 4.

Figure 6:

(A) Schematic diagram of APES of typical E⊗e linear DJT system along the appropriate configuration coordinate. E_JT is the Jahn-Teller energy that separates the total energy at the high symmetry geometry (conical intersection point) and the lowest energy in APES. (B) Calculated HSE06 DFT APES of the NV center in diamond in the electronic ³E excited state. The barrier energy (δ_JT) among the three minima with C_1h symmetry breaks the axial symmetry of DJT that can be handled by the so-called quadratic Jahn-Teller Hamiltonian (see Ref. [105]). Q_x and Q_y are the configuration coordinates associated with the effective double degenerate e_x and e_y vibration mode of=77.6 meV, respectively.

DFT HSE06 calculations yielded about 42 meV Jahn-Teller energy and about 9 meV barrier energy between the global minima of APES [105] that might be superior over the DFT PBE results, yielding smaller values [183]. Because of the finite barrier energy, the DJT should be solved in quadratic Jahn-Teller approximation. The exact solution can be expanded into series as

where the expansion is convergent with maximum of four oscillator quanta (n+m≤4), |n, m⟩ is the occupation representation of e_x and e_y vibrations, respectively, and c_nm, d_nm coefficients are obtained from the solution of the E⊗e DJT electron-phonon Hamiltonian [105]. We note that other systems may require a larger number of oscillator quanta (e.g. the product DJT system of neutral silicon-vacancy qubit in Ref. [185]). The tilde over Ψ labels the vibronic or polaronic nature of the wave function: the electron-phonon wave function cannot be factorized into electronic and phonon wave functions as assumed in BOA but will be a linear combination of such wave functions. Thus, one has to go beyond BOA in DJT systems, i.e. in strongly coupled electron-phonon systems. The vibronic spectrum starts with an E˜ ground state and then is followed by A˜1,A˜2, and E˜ vibronic levels by 39, 57, and 90 meV, respectively. We note that the A˜1 and A˜2 levels split due to the quadratic Jahn-Teller interaction. The vibronic solution will have a serious consequence in the effective spin-orbit splitting between the ³E spin levels and the ISC processes toward the ¹A₁ state, and it also manifests in the PL spectrum. We note that the singlet counterpart of |³E⟩, |¹E′⟩ experiences the same type of DJT effect as |³E⟩.

Phonons may couple nondegenerate states too, which is often called the pseudo-Jahn-Teller (PJT) effect. This occurs for the electronic |¹A₁⟩ and |¹E⟩, although the energy spacing between the two is more than 1 eV [177]. The dynamic coupling between these two states may be rationalized by invoking the symmetry breaking e_x,y phonon modes, which distort |¹E⟩ to |¹A′⟩ and |¹A″⟩, and |¹A₁⟩ becomes |¹A′⟩, so the |¹A′⟩ component of the distorted |¹E⟩ can couple to |¹A₁⟩ by Coulomb interaction. According to the DFT+CI-cRPA method, the contribution of |¹A₁⟩ into |¹E⟩ is about 2% in the distorted geometry [172]. At high C_3v symmetry, the contribution of |¹E′⟩ to |¹E⟩ is about 10%, which is an electron-electron correlation effect. We call the resulting multideterminant singlet state |1E˜〉. As a consequence of this fact, |¹E′⟩ brings the DJT effect into |¹E⟩, which is damped by the contribution factor. The full electron-phonon Hamiltonian contains the PJT and damped DJT effect that looks like a superlinear DJT Hamiltonian; although the PJT and DJT Hamiltonian is written in linear approximation for the sake of simplicity, the final equation is not linear because PJT and DJT distorts the system in a different manner [177]. The solution of PJT+DJT Hamiltonian is a set of polaronic wave functions. The combined |1​A˜1〉 ⊕|1​E˜±〉 states may transform as E, A₁, and A₂. The A˜2 vibronic states do not play a significant role; thus, we only show the expressions for the 1​E˜± and 1​A˜1 vibronic states as follows:

that govern the shape of the phonon sideband in the optical spectra of singlets. We label the symmetry adapted vibrational wavefunctions, e.g. |χ₁(A₁)⟩=|00⟩, or in general by |χ_i(…)⟩. The g_i coefficients are generally tiny and can be ignored. In contrast, the nonzero d_i and c′i (f_i and d′i) coefficients drive the ISC process, and they are also responsible for the shape of the PL spectrum of the singlets. The resultant electron-phonon spectra are depicted in Figure 7 and the values of the coefficients are listed in Ref. [177]. The lowest-energy dark A₁ excited vibronic level appears above the ground state |1​E˜〉 because of PJT effect and plays a crucial role in the temperature-dependent lifetime of |1​E˜〉.

Figure 7:

Vibronic levels of |1​E˜〉 and |1​A˜1〉.

The selection rules for the PL spectrum are indicated. Here, the ZPL energy of 1.19 eV between |1​E˜〉 and |1​A˜1〉 is not scaled for the sake of clarity. The calculated effective phonon frequency of the P_JT+D_JT Hamiltonian is=66 meV. The vibronic levels of |1​E˜〉 do not follow the ladder structure at all, whereas the vibronic levels of |1​A˜1〉 do show the ladder structure but significantly larger (92 meV) than ħω_e. These results are taken from Ref. [177]. The optical transition activated by uniaxial stress is observed in Ref. [186].

3.5 Optical properties

Similar to molecules, vibrations or phonons can contribute to optical transitions of point defects in solids. The optical transitions of molecules are often modeled by the Franck-Condon theory, which applies BOA. One can sketch the APES in the electronic ground and excited states together with the phonon states and levels as shown in Figure 8. The most simple case is that the APES can be described as a parabola in both the electronic ground and excited states, so the Franck-Condon theory perfectly works. Principally, the two APES parabolas may have different tangents or effective vibration modes, but they may be similar in the electronic ground and excited states. The Huang-Rhys theory assumes that two APES parabolas are identical, except for a shift in the minimum of the parabolas. In this case, the derivation of the phonon lineshape in the optical excitation spectrum is simplified. In the experiments, the DW factor can be measured, which has a direct relation to the Huang-Rhys factor S=−ln DW, which measures how many effective phonons participate in the optical transition on average. Generally, the larger the distance between the two APES parabolas, i.e. the larger movement of ions upon optical excitation, the more likely to incorporate phonons in the optical transition. The NV center has a broad absorption and emission spectrum (see Figure 8) even at low temperatures, and DW≈0.03 at cryogenic temperature, which means that about 3.5 effective phonons participate in the radiative decay.

Figure 8:

(Top) Total energy against configuration coordinate sketches the APES. The configuration coordinate is related to the position of ions. ħω_ν is the energy quantum of the (effective) ν phonon mode in quasiharmonic approximation, i.e. parabolic APES. In the Huang-Rhys approximation, the APES of excited state (top parabola) and that of the ground state (bottom parabola) are the same but the minima shifted by ΔQ_ν. The wave functions of the phonons are shown, and the overlap of the phonon wave functions of the two APES will appear as phonon sideband in the optical spectra. (Bottom) Absorption and emission spectra of NV center in diamond at low temperatures. The phonon modes coupled strongly to the optical transition can be identified as bumps in the phonon sideband of the corresponding spectrum. ZPL is the ZPL optical transition. Note that the two spectra are not identical in the phonon sideband as features at about −0.15 eV from ZPL are visible in the luminescence spectrum that are missing at +0.15 eV from ZPL in the absorption spectrum. These features are explained by the vibronic nature of the ³E state (see Refs. [105], [187]).

The construction of the absorption and emission spectra in the Huang-Rhys approximation requires the calculation of APES in both the electronic ground and excited states as well as the vibration modes of the point defect in the electronic ground state. DFT calculations showed that quasilocal vibration modes appear at about 65 meV in which the carbon and nitrogen atoms vibrate the most around the vacancy [182], [187], and these modes appear in the absorption and PL spectra as bump features (see Figure 8). DFT ΔSCF method can be used to calculate the APES in the electronic excited state. By forcing C_3v symmetry in the electronic excited state, Alkauskas et al. calculated the PL spectrum of the NV center in diamond [188]. They developed an embedding technique to calculate the contribution of long-wavelength phonons [188], and they obtained a good agreement with the experimental data, except for a small feature in the phonon sideband of the PL spectrum (see Figure 8) as pointed out in Ref. [105]. It was shown above that DJT occurs in the ³E state, which dynamically distorts the symmetry. By taking one of the parabolas in the sombrero APES, i.e. the distorted C_1h geometry, the PL spectrum can be calculated within the Huang-Rhys theory, and all the features in the phonon sideband of the PL spectrum were well reproduced [105]. The ab initio theory revealed that about 10% of the phonon sideband emission is associated with the e_x,y phonon modes that are responsible for the dynamic distortion in the ³E state [105]. This proves that the DJT feature indeed appears in the PL spectrum of the NV center. Furthermore, the DJT nature of the ³E state was considered to be responsible for the anomalous temperature dependence of the ZPL width [183].

The aforementioned asymmetry in the phonon sideband in the PL and absorption spectra is attributed to the DJT nature of the excited state. We note that a characteristic double peak in the absorption spectrum was associated with a Jahn-Teller feature [5], but the ab initio theory (see Ref. [188] and private communication with Audrius Alkauskas) produces this double peak by the contribution of A₁ phonons. The reason behind this asymmetry may be identified as the consequence of the dynamics of ions in DJT systems. By taking the low temperature limit for the sake of simplicity, the electron occupies the lowest-energy E˜ vibronic level in the ³E excited state upon illumination. This state may be viewed as the electron continuously tunnels from one of the global C_1h minima to the other. As the spontaneous emission of the photon from this excited state is an instantaneous process, the electron just stays in one of these global C_1h minima at any time when the radiative decay starts. The ³A₂ ground state has a high symmetry; thus, the phonon sideband of the PL process should be calculated as if the excited system was frozen in the distorted C_1h symmetry and it arrives to the high C_3v symmetry in the electronic ground state. As there is a continuum of E phonons that distort the symmetry from C_3v to C_1h in the ³E state, the Huang-Rhys theory therefore can be used to calculate the shape of the PL spectrum. In contrast, the absorption process is different. The rate of phonon absorption may go with the radiative lifetime of the excited state, which is associated with the optical transition dipole moment. This rate is Γ_rad≈13 MHz for the NV center in diamond (e.g. Ref. [104] and references therein). However, the tunneling rate of the electron in the ³E excited state is much faster. Bersurker [184] analyzed this for the lowest-energy solutions of E⊗e DJT system, and the Γ_tunnel goes as Γ_tunnel∝ΔE/h, where ΔE is the energy gap between the first excited and ground state vibronic levels and h is the Planck constant. This formula yields Γ_tunnel≈100 GHz. The exact solution is 112.6 GHz for the ³E state in the NV center, where the corresponding equations are Eqs. (7a) to (7c) in Ref. [189] (see also references therein). As a consequence, Γ_tunnel≫Γ_rad; therefore, the phonon sideband of the absorption spectrum should be calculated as the combination of the A₁ phonons in the Huang-Rhys approximation (predominant contribution) plus direct calculation from the ground-state vibration function toward the high-symmetry polaronic solution caused by the E phonons. This has not yet been published, to our knowledge, for the optical transition between the triplets, but a similar spectrum was published for the optical transition between the singlets [177], which will be discussed below.

The excited states and excitation spectrum of the triplets can be basically calculated by time-dependent DFT (TD-DFT), too [190], [191]. The usual linear response approximation was applied in the TD-DFT calculations of the NV center [191], [192]. We note that the accurate calculation of the spectrum requires the proper choice of the DFT functional. In 1.4 nm diamond molecular cluster model, the PBE0 [193] hybrid density functional provides quantitatively good results within TD-DFT [192]. The quantum mechanical forces acting on the ions in the electronic excited state can be calculated within the TD-DFT formalism. The experimental Stokes shift of the NV center of about 0.45 eV could be well reproduced by PBE0 TD-DFT calculations (see Supplementary Materials in Ref. [194]) despite the limitations of a molecular cluster model as explained above. This shows the predictive power of the TD-DFT method. TD-DFT calculations can be combined with molecular dynamics simulation to monitor the evolution of the ionic motions in the electronic excited state. It was found that phonons move the system from the phonon excited state to the zero point energy (lowest-energy phonon level) after photoexcitation from the electronic ground state within 50 fs, which agrees well with the observed decay time from pump-probe PL measurements [195]. The reason behind the peculiar ultrafast motion is tentatively associated with the DJT nature of the ³E excited state. The observed femtosecond electronic depolarization dynamics of the NV center [181] was associated with the nonadiabatic transitions and phonon-induced electronic dephasing between the two components of the ³E state based on ab initio molecular dynamics simulations.

We now discuss the PL spectrum of the singlets, in particular, the phonon sideband. The starting point is that the final state, |1​E˜〉, is very far from the quasiharmonic vibration spectrum; thus, the usual Franck-Condon approximation does not hold. In addition, Γ_tunnel≈31 GHz in the vibronic ground state of |1​E˜〉, which is several orders of magnitude faster than the inverse radiative lifetime of 1​A˜1 [172]. Therefore, the optical transition dipole moment should be calculated directly between the polaronic states, which can be written as [177]

It was found from direct calculation of the intensities in Eq. (7) that the optical transition to the first vibronic A₁ state of the 1​E˜ state is not allowed. However, there is a significant optical transition dipole toward the split E vibronic states about 45 meV. After switching off the small DJT effect in the electron-phonon Hamiltonian, only a single E mode appears with a smaller optical transition dipole moment. This clearly demonstrates that the small DJT effect does play an important role in understanding the optical features of the singlet states [177]. The simulated PL spectrum from ab initio wavefunctions is shown in Figure 9 that can be directly compared to the low-temperature experimental PL spectrum [178]. Clearly, the broad feature with the maximum intensity at ≈43 meV can be reproduced (red curve). It was found that the broad feature consists of two vibronic excited levels (see red ink text in Figure 7). The experimental intensity and the shape of this broad feature can be well reproduced by invoking our electron-phonon Hamiltonian (red curve). This theory is further supported by an uniaxial stress experiment on the PL spectrum, which showed the existence of a forbidden state at ≈14 meV [186]. This can be naturally explained by the calculated A₁ vibronic excited state (see blue ink text in Figure 7). This A₁ state will play an important role in the temperature dependence of the ISC rate where ≈16 meV phonon mode was deduced from the temperature-dependent ISC rate measurements in nonstressed diamond samples [33], which should be identical with the optically forbidden vibronic mode.

Figure 9:

Experimental PL spectrum of the singlets at low (black solid line) and room (dotted black line) temperatures compared to the simulated spectrum from the ab initio solution (red curve from Ref. [177]).

We note that the experimental spectra show a substantial and minor background at room temperature and cryogenic temperature, respectively. The simulation curve does not include background signal. The ZPL energy is now set to zero to easily read out the position of vibration features in the spectrum. 2, 5, and 10 meV Gaussian smearing was used for the linewidth of the ZPL and first and second vibronic emissions, respectively, where the width of the ZPL and vibration bands were read out from the experimental spectrum recorded at cryogenic temperature. This theory does not account for the features at 133 and 221 meV. These features seem to disappear at room temperature PL spectrum, and they may not belong to the NV center.

The absorption spectrum of the singlets should be calculated in a different manner. The PJT effect produces a barrier energy for the damped DJT APES energy surface of the 1​E˜ state; thus, at any event of absorption of the photon, the electron stays one of the distorted global energy minima. The photoexcited electron arrives to a highly symmetric state with equidistant 1​A˜1 vibronic levels that look like a perfect quasiharmonic solution. As a consequence, the absorption spectrum can be simulated from a frozen distorted structure in the ground state toward a highly symmetric structure in the excited state [177] that can be calculated within Huang-Rhys approximation. Finally, the simulated shape of the absorption spectrum is very different from that of the PL spectrum [177], in agreement with the experiments [178], [196]. We note that as the APES of the singlets was not directly calculated in this procedure [177], the sharp feature with energy above the phonon bands of diamond in the absorption spectrum therefore is not reproduced by this method, the feature of which was associated with the nitrogen-carbon local vibrations [196].

3.6 Magnetic properties

We previously discussed the orbitals of the NV center and their interaction with phonons. The group theory analysis was very powerful to identify the electronic structure of the ground and excited state orbitals. The fine structure of the NV center can be further analyzed by invoking the analysis of the C_3v double group, which treats spinor functions. As the NV center has an even number of electrons (holes) the single group rows in the C_3v double group appears as irreducible representation of the spin states [22], [23], [24], [25]. In particular, the ³A₂ ground state splits to A₁ms=0 and E ms=±1 states, whereas the ³E state splits to E_x,y ms=0, E_1,2ms=±1, and A₂ and A₁ms=±1. These splittings are intrinsic to the C_3v symmetry defect, and not the external fields are responsible for their presence. Basically, spin-orbit and dipolar electron spin-electron spin interactions may introduce energy spacing between spin levels. In the ³A₂ ground state, the effective angular momentum (L) of the orbital is zero; therefore, the spin-orbit interaction (λzL^zS^z) is zero, where λ_z is the strength of interaction along the symmetry axis z. The ZFS D between A₁ and E is therefore caused by dipolar electron spin-spin interaction as already stated in Eq. (1). In contrast, L=1 can be envisioned for the ³E state; thus, both spin-orbit and dipolar electron spin-electron spin interactions are present. The spin-orbit coupling and the zz component of the dipolar electron spin-electron spin tensor split the ms=0 and ±1 levels in the ³E state. Further splitting between A₁ and A₂ms=±1 levels is caused by the x²-y² component of the dipolar electron spin-electron spin tensor [24], [25], whereas its xz component can couple the E_x,y ms=0 with E_1,2ms=±1 states [24]. The latter is responsible for the very weak radiative decay from ms=0 toward the ground state via the E_1,2 state (see Section 3.6), which is the base for realizing the so-called Λ-scheme for single-shot readout [33]. These are the basic ZFS parameters in the ground and excited state spin manifolds.

The spin states and levels behave in a strikingly different way as a function of temperature [103], [197]. In the ground state manifold, the ZFS D constant decreases slightly [−74.2(7) kHz/K] but well resolved [197]. In contrast, the fine structure of the ³E state can be well observed at cryogenic temperatures with λ_z=5.33±0.03 GHz [103], [198] and D=1.42 GHz, DA1,A2=Dxx−Dyy= 3.1 GHz [103] with DE1,2,Ex,y=Dxz= 0.2 GHz coupling parameter [38], but only ZFS caused by D_zz remains observable at T=20 K and temperatures above, whereas the gap between A₁ and A₂ as well as the spin-orbit-related splitting (λ_z) completely disappears [103]. Understanding these features requires to calculate the ZFS parameters, i.e. the spin-orbit energy and the dipolar electron spin-electron spin tensor as well as understanding dynamical effects in the ³E state. The latter will be discussed in Section 3.7.3.

The spin Hamiltonian of the spin-spin dipolar interaction may be written as

where r_ij=r_i−r_j. The 3×3 D-tensor can be diagonalized to find the spectrum and spin eigenstates. In the C_3v symmetry, the D_zz component lies in the symmetry axis, which also sets the quantization axis of the electron spin along the symmetry axis of the NV center. In the ground and excited state manifold, the ZFS is D=(3/2)D_zz and the corresponding spin levels can be computed as D(Sz2−S(S+1)3) for S_z={0,±1} with S=1; thus, the components, e.g. zz component, of the tensor should be calculated. The D-tensor is associated with the two-particle spin density matrix, n₂(r₁,r₂), which can be approximated using the Slater determinant of the KS wave functions ϕ of the considered system, so that n₂(r₁, r₂)≈|Φ_ij (r₁, r₂)|², where Φij(r1, r2)=12(ϕi(r1)ϕj(r2)−ϕj(r1)ϕi(r2)) and then

where r_a,b=(r₁−r₂)_a,b and χ_ij is either 1 or −1 for KS i, j states of the same or different spin channels, respectively. Note that in DFT the KS states are not spin restricted as mentioned previously. Consequently, not only the unpaired KS states but also the rest of the occupied states can contribute to the spin density and the ZFS [199]. The first implementation of Eq. (9) in the DFT supercell version with LCAO basis set was reported in Ref. [199]. The first implementation into plane-wave supercell code was done by Ivády et al. [200], which only used pseudo-wave functions in the vasp code. The theory of all-electron PAW version of the D-tensor was reported by Bodrog and Gali [201]. The all-electron version was then implemented by Martijn Marsman into vasp code in 2014. Other implementations also followed later: both the pseudo-density version [202] and the all-electron versions have been recently implemented into Quantum Espresso code [203], [204] or an adaptive finite-element 3D grid code [205].

The calculated D constant in the ground state was within 3% with respect to the low-temperature value at ≈2.87 GHz [200] in Γ-point 512-atom supercell calculation. By increasing the supercell size, this value did not change, but a larger supercell might be needed for similar defect quantum bits [206]. The ΔSCF method can be used to generate the KS states for calculating the components of the D-tensor in the ³E excited state. By following the previous group theory analysis [24], [25], the calculated D=1.61 GHz, DA1,A2=Dxx−Dyy= 1.95 GHz, and DE1,2,Ex,y=Dxz=0.15 GHz do not show such an excellent agreement with the derived parameters from PLE measurements [103], [207], which resulted in 1.42, 3.1, and 0.2 GHz, respectively. It is not yet clear why there is a larger discrepancy with respect to the experimental data in the excited state than that in the ground state. The ΔSCF method should work well for the ³E state, and tests on the size of the supercell do not significantly change the results (Á. Gali, unpublished results).

The spin-orbit interaction Hamiltonian in zero-order approximation can be written as

where V is the nuclear potential energy, m_e is the electron mass, and p^i and S^i are the momentum and spin of electron i. The elements of the orbital operator vector O^=∇iV×p^i can be calculated from the KS orbitals, which defines the effective angular momentum of the NV center. Note that the crystal field of a solid breaks the spherical symmetry of the spin-orbit interaction. In C_3v symmetry, the interaction Hamiltonian can be rewritten as [24]

where λ_⊥ and λ_z are the basal or nonaxial and axial parameters of the interaction, respectively. The former flips the orbital and spin state; thus, it does not play a role in the solution of the double group of C_3v symmetry but can play a role in the ISC by coupling ms=±1 of the triplet to a singlet state [24], [25]. The latter will only shift spin levels in the ³E triplet or may couple the ms=0 of the triplet to a singlet state [24], [25].

The axial spin-orbit interaction strength in the ³E excited state can be obtained both from total energy difference calculations and from the splitting of KS states in ΔSCF noncollinear spin-orbit DFT calculations [105]. It has been shown that the two approaches gave the same result. The latter is more practical because it only calls a single calculation but should be done in the Γ-point of the BZ, as the dispersion and splitting of the defect states in low symmetry k-points may be larger than the axial spin-orbit splitting. To accurately determine small values of the axial spin-orbit interaction, the calculations require high numerical convergence and accuracy as the spin-orbit energy falls in the μeV regime. Finite size effect turned to be crucial for spin-orbit interaction calculations [105] (see Figure 10). Thiering and Gali [105] attributed the observed finite size effect to the overlap of the defect states and used an exponential fit to eliminate the supercell size dependence of λ_z.

Figure 10:

Convergence of the axial spin-orbit coupling parameter λ_z in the excited state of the NV center [105].

Horizontal axes shows the supercell size and the number of carbon atoms in the defect free supercells within Γ-point sampling. DFT PBE calculations prove the exponential decay of the value that can be used to fit the converged value for DFT HSE06 calculations.

The calculated HSE06 λ_z=15.78 GHz is about three times larger than the data derived from PLE measurements at cryogenic temperature [103]. This difference is not a discrepancy but rather the consequence of the DJT nature of the ³E state as stated in Eq. (5). Ham already suggested that the effective spin-orbit energy will be reduced in such systems with a Ham reduction factor p [208]. In the particular ³E state of the NV center, p=∑nm[cnm2−dnm2], which represents the mixture of the E₊ component with the E₋ component of the ³E state that results in the quenching of the effective angular momentum (see the derivation in Ref. [209]). The c_nm and d_nm coefficients were taken from the solution of the E⊗e DJT electron-phonon Hamiltonian, which finally results in p=0.304 with an effective spin-orbit splitting of 4.8 GHz [105].

The ISC between ¹A₁ and ³E states is associated with the perpendicular component of the spin-orbit operator. Unlike the case of λ_z, one has to apply an approximation to calculate this, namely, that the KS wave functions building up the ³E state and the ¹A₁ multiplet do not change. This enables to compute λ_⊥ab initio using the a₁ and e_+,− KS wave functions of the NV(–) in the ³E excited state as 〈a1|H^SO|e+〉/2. The converged intrinsic value is λ_⊥=(2π)56.3 GHz in rad/s unit, which is relatively large, and it is not damped by DJT because the ¹A₁ state is not part of the DJT effect. We note that the nitrogen contribution is minor in λ_z, whereas it is significant in λ_⊥ that explains the large anisotropy between λ_z and λ_⊥.

In the NV center [see Eq. (1)], the ¹⁴N nuclear spin has a quadrupole splitting, C_Q, which can be written as

where h and e are the Planck constant and the charge of the electron and V_zz is the gradient of the electric field changes around the nitrogen nucleus in the z-direction (chosen symmetry axis of the defect). Q_N is the nuclear electric quadrupole moment of ¹⁴N, which scatters between 0.0193 and 0.0208 barn in the literature [210]. V_zz can be calculated ab initio within the PAW framework. In this calculation, the criterion for the wave functions (plane-wave cutoff) and the forces should be set accordingly to achieve convergent V_zz (see the Supplementary Information in Ref. [10]). The calculated C_Q value for the NV center is −5.019±0.19 MHz, where the uncertainty comes from the uncertainty in the experimental data of Q_N. Nevertheless, the ab initio result is very close to the experimental data and could be used to follow the change in the quadrupole moment in NV defect as a function of the charge state [10].

The last remaining terms in Eq. (1) at zero magnetic field are related to the hyperfine interactions. Hyperfine interaction describes the interaction between the electron spin and the nucleus spins that can be generally written as

where A is the hyperfine tensor and I^ is the nuclear spin vector operator. Unlike the case of electron spin-electron spin interaction, where the Pauli exclusion principle does not allow two electrons with the same spin at the same place, the electron spin may be localized at the place of the nucleus spin. This is called the Fermi-contact term, whereas the dipole-dipole term is familiar to what we learned for the electron spin-electron spin interaction. These can be respectively written as

where n_s(r) is the electron spin density, r is the vector between the electron spin and nuclear spin at R, g_n is the nuclear g-factor, and μ_n is the nuclear magneton for a given nucleus n. We note that the nitrogen atom and distant carbon atoms sit in the symmetry axis by creating a hyperfine field that preserves the C_3v symmetry. For those nuclear spins, the hyperfine tensor is diagonal with diagonal element A_xx=A_yy=A_⊥ and A_zz=A_||. These parameters can be expressed by the Fermi-contact term a and a simplified dipolar coupling term b as A_⊥=a−b and A_||=a+2b. The other nuclear spins reduce the symmetry; thus, all the principal values of the hyperfine tensor differ for those nuclear spins. One can recognize from the definition of the Fermi-contact term that the spin density at the nucleus site should be calculated very accurately. Thus, all-electron methods should be applied: either treating the core electrons explicitly (typical for GTO codes [122]), PAW method for plane-wave codes [3], [11], [115], [211] or finite-element numerical method on a grid [205].

The hyperfine constants of the NV center were first determined by plane-wave PAW DFT LDA method [3]. It was found that the spin density is not evenly distributed around the core of the defect but rather elongated in the plane of the three-carbon dangling bonds where the amplitude of the spin density decays like sin x/x function, where x is the distance from the vacancy (Figure 11A and B). This spin density distribution has an analog of a Friedel oscillation with the NV center as an impurity. This spin density distribution is responsible for the fact that direct quantum mechanical interaction of the NV center with acceptor defects is much stronger in the basal plane of the three-carbon dangling bonds of the NV center, i.e. in the (111) plane than that along the ⟨111⟩ direction [155]. Another important consequence of this finding is that the Fermi-contact term cannot be neglected for ¹³C spins in the (111) plane near the NV center, but some proximate ¹³C spins have small hyperfine constants that lie in the node of the spin density oscillation. In this study [3], some ¹³C spin quantum bits were identified that were employed in a previous entanglement measurement between the electron spin and nuclear spins of the NV center [16].

Figure 11:

Calculated spin density in a 512-atom supercell in Ref. [3].

The isosurfaces of blue, green, yellow, orange, and red lobes are −0.014, −0.007, +0.014, +0.126, and +0.478, respectively. (A) ³A₂ state top view of the (111) plane. The threefold rotation symmetry of the C_3v point group can be well identified with symmetrically equivalent ions around the vacancy. The red lobes show the position of the three carbon atoms near the vacancy with the largest localization of the spin density. (B) ³A₂ state side view of the (110) plane. The spin density shows an oscillation in the plane of the three nearest carbon atoms around the vacancy. (C) ³E side view of the (110) plane. The spin density is significantly localized on the nitrogen atom below the nearest three carbon atoms around the vacancy.

The spin density significantly redistributes in the ³E excited state and is significantly localized around the nitrogen atom (Figure 11C). As a consequence, the hyperfine constants on the nitrogen spins increases up to ≈50 MHz, which broadens the resonance condition of external magnetic fields for transferring the polarization of the electron spin toward the nitrogen nuclear spin in the excited state (see Refs. [65], [66], [67]).

This type of calculations was repeated by the PAW DFT PBE method [11], [66], where the agreement between the reported hyperfine constants of ¹⁴N and ¹³C nuclear spins was very good [6], [8]. Later, it was found that this was partially fortuitous agreement [211] because of the PAW implementation with frozen core orbitals [115]. The core polarization of carbon 1s orbital significantly contributes to the Fermi-contact term, in particular, where the spin density is much localized [211]. Generally, core polarization correction is required in the PAW implementation of accurate hyperfine tensor calculation [211], [212]. DFT HSE06 calculations localized the wave functions more around the vacancy than DFT PBE does, but core polarization correction will reduce the Fermi-contact term of ¹³C spins, resulting in an excellent agreement with the experimental data [211]. It was shown [211] that cancellation of errors (DFT PBE underestimation of spin density and neglect of core polarization) does not generally occur for other defect quantum bits; thus, the implementation of core polarization in the PAW formalism is crucial. In a recent study [122], nonflipping proximate ¹³C spins in the NV center going from the ground state to the excited state were identified in an 510-atom molecular cluster model by means of hybrid functional calculations of hyperfine tensors as implemented in an all-electron GTO code that can be used as long-living quantum memories.

We note that orbital hyperfine interaction may contribute to the hyperfine dipole-dipole interaction in defects with an effective angular momentum of the associated orbitals (e.g. Ref. [213]), which was implemented into a Green function-based code and applied to nickel defects in diamond [214] but not to known defect quantum bits.

3.7 Radiative and nonradiative rates: a complete theory on the optical spinpolarization loop

By having the magneto-optical parameters in hand, one can discuss the ab initio description of ODMR and PDMR readout of the NV center in diamond. In the ODMR readout, the interplay between the radiative and ISC decay rates from the ³E excited state plays a crucial role that will be first discussed below. In the PDMR readout, the direct photoionization and Auger recombination processes also enter, which we will describe next.

Radiative rates can be readily calculated from KS DFT by taking the corresponding a₁ and e_x,y KS orbitals of the NV center and calculating the 〈ex,y|er^|a1〉 transition dipole matrix element (e.g. Refs. [20], [105]). Using the Wigner-Weisskopf theory of fluorescence, the radiative lifetime τ can be calculated as [215]

where n_r is the refractive index of diamond, ħω is the excitation energy, c is the speed of light, and er is the optical transition dipole moment. One can realize that the ez component will be zero, and only the perpendicular components will give nonzero results as was also found in the experiments [216]. The polarization property of the absorbed and emitted photons is naturally obtained by KS DFT calculations. Finally, Eq. (15) results in about 8 ns radiative lifetime, which is relatively close to the experimental one (~12 ns that corresponds to about (2π)13 MHz rate; see Refs. [104], [217]). For the triplet optical transition, this is a very good approximation because of the relatively simple character of the electronic ground and excited states. The calculation of the optical transition between the singlet states has just become feasible by the DFT+CI-cRPA method [172], where the multiplet wave functions are considered in the optical transition dipole moment calculation. The calculated radiative lifetime of the ¹A₁ state is 1878 ns, which is two orders of magnitude longer than that for the triplet [172]. This ab initio result sheds light on the experimental fact that not just the competition with nonradiative decay and the low fractional population of the singlets but the inherently tiny optical transition dipole moment is responsible for the very weak infrared PL signal of the NV center [178]. This agrees well with the conclusion of absorption measurements between the singlets [196]. The very recent pump-probe PL measurement has observed 100 ps lifetime of the ¹A₁ state [218]. The combination of this recent experimental data with the ab initio radiative rate [172] and observed absorption [196] indicates that the nonradiative decay from the ¹A₁ state to the ¹E is much faster than the radiative one. This should be verified by future ab initio calculations of the nonradiative rate between the singlets.

We note that these calculations inherently assume that the Franck-Condon principles hold for the optical transitions, i.e. the optical transition dipole moment is independent from the ionic movements. This is a valid approximation for the optical transition between the triplets of the NV center in diamond, but it does not necessarily hold for all types of vibrations and solid-state defect quantum bits, in which the Herzberg-Teller interaction is considerable [219]. We showed above that the fluorescence of the singlets of the NV center in diamond does not follow the Franck-Condon principles either, and one should use polaronic wave functions in the calculation of the optical spectrum.

The nonradiative rates are associated with the spin-orbit interaction between the triplet and singlet states mediated by phonons, i.e. ISC processes. Thus, the ISC rates from the ³E spin states to the ¹A₁ state and from the ¹E state to the ³A₂ spin states should be calculated. In the Franck-Condon type of approximation, the ISC rate from the A₁ state of the ³E state may be calculated [104] as

where λ_⊥ is given in rad/s unit, Δ is the energy spacing between |A₁⟩ and |¹A₁⟩, and δ is the Dirac delta function. |χ₀⟩ is the ground vibrational level of |A₁⟩, and |χ′νn〉 is the vibrational level of |¹A₁⟩ with energies ν_n above that of |¹A₁⟩. The group theory implies that ISC can only occur between |A₁⟩ and |¹A₁⟩ (see Figure 1B). However, three different ISC rates were observed in the experiments at cryogenic temperature in which the ratio of two rates is about 0.52±0.07 and the third rate was much smaller than the other two [104]. We already showed above that the ³E state is not a pure electronic state, but the e_x,y phonons couple the components of the ³E state. By solving the E⊗e electron-phonon Hamiltonian and analyzing the solution by the group theory, one finds A˜1,A˜2, and E˜1,2 vibronic states as eigenstates, where the contribution of the electronic |A₁⟩ in | A˜1〉,| E˜1,2〉, and | A˜2〉 is c_i, d_i, and f_i, respectively, where the running variable i groups the phonon wave functions for given representations of phonon wave functions and n_i is the sum of quantum numbers of e_x and e_y phonons [105]. Thiering and Gali demonstrated that the DJT nature of |³E⟩ naturally explains the three observed rates that can be written with modifying Eq. (16) as

where |χ′νn〉 is now restricted to the a₁ vibrational levels of |¹A₁⟩ [105]. The ab initio solution yielded f₁<0.001, which explains the very low ΓA2 at cryogenic temperature. In contrast, c₁=0.578 and d₁=0.331 are comparable. As a consequence, ΓA1/ΓE12≈0.5 for Δ=0.4 eV, which agrees well with the experimental findings [104], [105]. In contrast, the absolute values of the ΓA1 and other rates were significantly larger than the experimental data by taking λ_⊥=56 GHz from first-principles spin-orbit matrix element from a₁ and e_x,y KS wave functions. It is suspected that KS wave functions building up the ³E and ¹A₁ states may change, so it could be not a good approximation to calculate λ_⊥ from KS wave functions fixed at ground-state electronic configuration.

The next step is to calculate the ISC rate from the ¹E state to the ³A₂ spin states in the lower branch that will also determine the lifetime of the ¹E state. It is known from experiments that the lifetime of |¹E〉 is temperature dependent and varies between about 371 and 165 ns going from cryogenic temperature to room temperature [220]. The temperature dependence could be well understood by a stimulated phonon emission process with an energy of 16.6±0.9 meV [220]. On the other hand, the ODMR contrast does not seem to vary significantly as a function of temperature. From PL decay measurements, it was derived that the decay rates toward ms=0 and ±1 favors ms=0 only by about 20%, which was a surprising result but does not contradict with the measured ODMR contrast because of the high spin selectivity in the upper branch [220]. Indeed, the nature of the ¹E state is complex: it is a polaronic 1​E˜ state in which the character of ¹A₁ and ¹E′ states appear. The former and latter links 1​E˜ to ms=0 and ±1 spin states of ³A₂, respectively, with the corresponding ISC rates of Γ_z and Γ_⊥=Γ±+Γ_∓, respectively. These rates may be described as [177]

where the summation over all vibration wave functions of ³A₂ collapses to the number of |χi(e±)〉 vibration modes in the phonon overlap integral. Here, d_i coefficient is responsible for the contribution of the electronic ¹A₁ state in |1​E˜〉 that is linked to ms=0 of ³A₂ by λ_z. ∑ is the energy spacing between ³A₂ and 1​E˜ states of about 0.4 eV, and (1−C²)≈0.1 is the contribution of the ¹E′ state in the 1​E˜ state. The other rates can be expressed as [177]

and

It was again found that λ_⊥/λ₂=1.2 provides consistent absolute total ISC rate (TE−1=Γz+Γ⊥) of about 2.7 MHz with ∑≈0.4 energy gap, and the λ_⊥ derived from ab initio spin-orbit matrix element between KS a₁ and e_x,y states is largely overestimated. The calculated Γ_z/Γ_⊥ is about 5.5, which means that the ab initio theory predicts [177] about 84% selectivity toward ms=0 state over ms=±1 state in the ISC process in the lower branch of the optical spinpolarization loop. In a recent experimental study [221], the refined modeling of the signals resulted in Γ_z/Γ_⊥=4±0.5, which agrees well with the ab initio theory. This means that relatively high spin selection appears in the ISC process between the 1​E˜ state and ³A₂ ground spin states. Finally, the temperature dependence of the lifetime T_E of the 1​E˜ state can be calculated using the calculated vibronic levels of 1​E˜ and the ISC theory in Eqs. (20) to (22) (see Figure 12). The Boltzmann occupation of vibronic levels was used at the given temperature to compute ISC rates in Ref. [177]. A very good agreement was found with the experimental data [33] as the calculated lifetime is reduced from 370 ns at cryogenic temperatures down to 171 ns at room temperature to be compared to 371±6 and 165±10 ns, respectively. The ab initio theory calculations showed quantitatively that the vibronic state associated with the optically forbidden phonon feature at ≈14 meV in the PL spectrum plays a key role in the temperature dependence of the ISC rates [177] where this connection was previously hinted in Ref. [222]. The calculated Γ_z/Γ_⊥ is only reduced by ~5% going from cryogenic temperature to room temperature in the simulations [177], which means that the spinpolarization efficiency per single optical cycle does not degrade significantly as a function of temperature. These results demonstrate that the ab initio theory developed by Thiering and Gali can account for the intricate details of the ISC processes in the NV center in diamond and reproduce the basic experimental data.

Figure 12:

The calculated lifetime [177] of the singlet shelving state is plotted as a function of the temperature with the observed lifetimes for two single NV centers (dot and triangle data points with uncertainties) taken from Ref. [33].

We emphasize that Eqs. (17) to (22) complete the theory of optical spinpolarization loop of the NV center in diamond where the corresponding matrix elements and coupling parameters could be mostly determined from first-principles calculations [105], [172], [177]. The remaining issues are discussed in Section 4.

Next, we discuss the photoionization of the NV center and related defects that are important in the PDMR readout of the NV center’s spin state. The idea of PDMR readout is based on the standard ODMR measurements of a single NV center. The phonon sideband photoexcitation of the NV center (typically 532 nm wavelength) with high laser power can lead to the unintentional ionization of the defect via the absorption of a second photon before decay from the ³E state toward the ground state. The ab initio theory predicts that ≈2.7 eV requires to ionize NV center from the ground state to neutral NV defect at cryogenic temperature. The absorption of the second photon at the ³E state likely occurs at the zero-point energy as the phonon relaxation in the ³E state is extremely fast [195], which is 1.945 eV above the ground state (ZPL energy of the NV center). Thus, the absorption of the second photon would promote the electron to ≈4.3 eV above the ground-state level of the NV center, which is much larger than the ionization threshold. The two-photon absorption and the corresponding ionization of the NV center will convert it to neutral NV [18]. The ab initio theory proved that this two-photon ionization is an Auger recombination process [20]. It is well known in semiconductor physics that two types of ionization of defects may occur: (i) direct ionization or (ii) Auger recombination. In Auger recombination, multiple defect levels and states may play a role. In the NV center, the direct ionization from the e_x,y orbital in the ³E state competes with the Auger process where the high-energy electron promoted to the conduction band will fall back to the hole left in the a₁ state, and the energy gain is used to kick out an electron from the e_x,y level [20]. The Auger recombination rate, i.e. the inverse process of impact ionization, can be written as a Coulomb interaction between the initial exciton state and all the possible final exciton states. The initial state i is an excitation from orbital j to a with spin ↓; here, j is the hole on the a₁ level and a is a state close to the conduction band edge, which can be simplified as ϕ_j↓ϕ_a↓ electron-hole KS wave functions. The energy of this exciton is approximated as ε_a↓−ε_j↓, where ε_rσ is the KS level of state rσ. The final exciton state is approximated similarly as ϕ_k↑ϕ_b↑ electron-hole KS wave functions, where k denotes the hole on the double degenerate e_x,y orbital, b is the electron above the conduction band minimum, and the energy of the exciton is ε_b↑−ε_k↑. To obtain the rate of transition, we need to sum over all possible final states: the degenerate e_x,y level (index k) and the electron in the conduction band minimum (index b). The final formula may be given as (see Supplementary Materials in Ref. [20] and also Ref. [223])

We obtain a single Γ by averaging over initial states:

where ħω is the sum of two energy pulses (two-photon absorption) and V can be generally expressed as

where ϕ_s,r,u,t are the corresponding KS wave functions. Vörös et al. [223] implemented the computation of this rate into a plane-wave supercell code [203], which was applied for the NV center in diamond, and ε_r was simply estimated from the refraction index of diamond in this particular calculation [20]. This could be the first ab initio calculation of the Auger recombination rate for point defects in solids. The direct ionization rate can be calculated like the usual absorption rate with optical transition dipole operator just the final state is at the conduction band minimum. The calculated Auger recombination rate was 800 ps, whereas the direct ionization rate was 0.5 μs. A very important consequence of this finding is that the two-photon absorption of the NV center results in the “ground state” of neutral NV defect; thus, no PL signal of neutral NV defect is expected. The typical 532 nm excitation wavelength can excite the neutral NV defect for which the same type of two-photon absorption occurs; finally, it ends at the ground state of the NV center [20]. It was found that the two-photon absorption of the neutral NV defect is very effective at the excitation wavelength of its ZPL wavelength at 575 nm (yellow light). One can conclude that the full cycle of photoexcitation and ionization of the NV center will lead to emission of an electron and a hole, and the final state is the ground state of the NV center.

As the ejection of the electron goes through the ³E excited state in which the spin states have different lifetimes, the rate of ionization (absorption probability of the second photon) should be also spin dependent. This idea was tested on a diamond device where electrodes were built on the diamond structure for the observation of the photocurrent upon the photoionization of the NV centers [14]. They observed a contrast in the photocurrent as they hit the microwave resonance of the electron spin in the ground state of ensemble NV centers [14], [96], i.e. PDMR readout. The full ab initio description of a PDMR cycle would require to calculate the Auger recombination rate for the neutral NV defect, which has not yet been reported. In contrast, the emission of the electron from the NV center and its spin selectivity is now fully understood using the ab initio theory. CI-cRPA calculations showed little or no photoionization from the 1​E˜ state toward the conduction band edge by 532 nm green excitation [172]. Thus, all the important rates are in hand to calculate the PDMR contrast, as the rates for ODMR readout can be used in combination with the Auger recombination rate.

In the first PDMR measurements [14], the PDMR contrast was relatively low compared to ODMR contrast because of the background current coming from mostly the P1 center in diamond. Therefore, it is highly important to understand the photoionization of the NV center and other parasitic defects to optimize the PDMR readout technique. Recent ab initio calculations have found that green light excitation much favors the ionization of the P1 center over that of the NV center; however, the relative ionization rates can be pushed toward increasing that of the NV center if blue light excitation is applied [95]. In these calculations, the valence and conduction bands were included in the computation of the absorption rates (not just the band edges); therefore, the BZ sampling should go up to 6×6×6 MP k-point set for convergent calculations in a 512-atom supercell [95]. Indeed, dual-beam photoionization of the diamond sample with an ensemble of NV centers resulted in an 3× enhanced PDMR contrast over the single-beam photoionization scheme, where blue light was employed for direct photoionization and green light was applied to induce spin-selective nonradiative decay from the ³E excited state [95]. This again proves the predictive power of ab initio calculations and its role for the improvement of the control and readout of solid-state defect quantum bits. Later on, it has been found that lock-in techniques with usual green-beam excitation can result also in an enhanced PDMR readout contrast [97]. Other parasitic defects producing photoionization bands in the diamond samples could be identified by ab initio calculations [171], which can be rather used as a resource, i.e. PDMR quantum bits, which might not be visible in the ODMR readout but could be effective in the PDMR readout [171].

We note that the carrier capture cross-section rates of NV defect and/or parasitic defects in diamond can be important in understanding the atomistic processes and improving further the PDMR readout. The ab initio framework already exists to calculate these rates [224] but not yet has been applied to the NV center in diamond. The carrier capture cross-section together with the Auger recombination process may play an important role in the electroluminescence of NV defect, which only could successfully employed for the neutral NV defect [225].

3.8 Effect of external perturbation

3.8.2 Electric field and strain

In the presence of electric field and strain, the ground state spin Hamiltonian in Eq. (1) is modified. Again, there are symmetry-conserving interactions that only shifts the levels, whereas a symmetry-breaking interaction may mix the spin states. Accordingly, the z-component of the electric field vector E only shifts the spin levels (transforms as full symmetry A₁ representation; E_z) whereas the x and y components (transform as E representation) can mix the ms=±1 and ms=0, ±1 levels that can be written as

(26)H^el=Hel0+Hel1+Hel2=hdzEzSz2+hd′⊥[Ex(SxSz+SzSx)+Ey(SySz+SzSy)]−hd⊥Ex(Sx2−Sy2)+hd⊥Ey(SxSy+SySx),

where h is the Planck constant. Eq. (26) has been known as a linear Stark shift for C_3vS=1 EPR centers and was derived specifically for the NV center by Doherty et al. [233]. The coefficients d_⊥=17 Hz cm/V and d_z=0.35 Hz cm/V have been inferred in the experiment of Ref. [234]. However, to our knowledge, the coefficient d′⊥ has not been quantified experimentally or theoretically; nevertheless, it is expected [233] to have the same order of magnitude as d_⊥. We note that the effect of d′⊥ is often neglected as in the zero magnetic field it is suppressed by the relatively large D=2.87 GHz gap between the ms=0 and ±1 levels. In contrast, under constant magnetic fields close to the ground-state level anticrossing point, this effect can be dominating, whereas the mixing between ms=±1 (d_⊥ term) will be suppressed.

The first-principles calculations of electric field are not straightforward in the 3D supercell model as the electric field breaks the periodicity of the system in the direction of the electric field. One possible solution is to apply a slab model for modeling the NV center with a surface that is perpendicular to the direction of the electric field. In Ref. [235], the NV center in diamond was modeled by a 990-atom 2-nm-thick diamond (111) slab, where the applied electric field was set to 0 to 0.1 V/Å along the symmetry axis of the defect in diamond slab with –OH termination. In this interval, the response was found to this perturbation to be linear on the calculated D-constant and resulted in d_z=0.76 Hz cm/V (see Ref. [235]), which falls to the order of magnitude deduced from experiments. The disadvantages of this method are that (i) the calculation of the perpendicular component of the coupling parameters requires another slab model with an appropriate surface orientation and termination and (ii) the surface termination may induce “artificial” polarization around the NV center even when is buried deepest in the slab model.

Recently, an alternative approach was used to study the effect of electric fields on the properties of the NV center [236]. In this particular case, the coupling of external electric field to the optical transition was examined, which is also known as Stark shift [237]. The Stark shift is generally an undesirable effect for quantum emitters because it can cause spectral diffusion (e.g. Ref. [20]): during optical excitation of the NV center, nearby defects may be ionized that creates different electric fields around the NV center that results in a shift of ZPL of the NV center in subsequent optical excitation cycles. Another consequence is that zero (magnetic) field ODMR signals of ensemble NV centers will show a splitting due to the various charged defects around the NV centers [238]. For determining the Stark shift of the NV center, the permanent dipole moments should be calculated in the ³E excited state and the ³A₂ ground state.

The alternative approach of calculating this quantity relies on the modern theory of polarization [239], [240], [241], [242] that can be applied to periodic systems. The polarization p can be written as

(27)p=ife8π3∑n=1M∫BZdk〈ukn|∇k|ukn〉,

where M is the number of occupied bands, f is the occupation number, e is the elementary charge of the electron, and the wave functions (u_k ) are cell periodic and periodic in the reciprocal space. Using density functional perturbation theory, ∇_k |u_kn〉can be calculated from the Sternheimer equations with similar self-consistent iterations as in the self-consistent field DFT

(28)(Hk−εkn)∇k|ukn〉=−∂(Hk−εkn)∂k|ukn〉.

By applying this theory to the NV center, the coupling parameter of the electric field to the optical transition is Δp=|〈³E|p|³E〉–〈³A₂|p|³A₂〉|=0.903eÅ≈4.33 Debye, and it directs toward the NV line along the ⟨111⟩ direction. This value corresponds to 2.18 MHz cm/V, which means that a point-like charge (e.g. positively charged N_s) at 50 nm distance from the NV center shifts the ZPL energy by 2.4 GHz. This result agrees well with previous experimental results [243], where |Δp|≈1.5 Debye was deduced for the NV center in diamond.

Basically, the advanced theory of polarization enables to calculate the d_z, d_⊥, and d′⊥ components of Eq. (26) in a 3D supercell model as the electric field can be added to the Hamiltonian at an arbitrary direction. Nevertheless, this has not yet been applied to the NV center in diamond.

The ³E excited state is much more sensitive to the presence of electric field or strain as the degenerate orbitals may split in the presence of symmetry breaking fields, and they are intertwined [24], [25]. Therefore, the effect of strain is first discussed before discussing the effect of electric field in the ³E excited state.

We note that the effect of hydrostatic pressure on the spin levels in the ground state was first studied in experiments [244] and the ab initio theory [200] before determining the spin-stress coupling parameters. In the theoretical study, the lattice constant of the 512-atom supercell was varied to induce the appropriate pressure on the NV center and the D-constant was calculated within pseudo-wave function approximation [200]. The experimentally deduced linear coupling coefficient is 14.58 MHz/GPa, which was measured at room temperature [244]. The results of ab initio simulations are valid at T=0 K and yielded 10.3 MHz/GPa (see Ref. [200]). It was found in the ab initio simulations that the local relaxation of ions and adaption of the wave functions to the new ionic positions upon external pressure plays a crucial role in the final response; thus, ab initio simulations are essential in the calculation of the spin-stress parameters. We further note that the pressure shifts the energy of the ZPL optical transition too by 5.5 meV/GPa [245]. DFT PBE ΔSCF calculations yielded 5.75 meV/GPa shift [246], in good agreement with the experimental data. Furthermore, they predicted that the emission in the phonon sideband will be significantly suppressed and redistributed toward the ZPL emission under giant (100 GPa) pressure [246]. This has not yet been confirmed in experiments to our knowledge. Next, we turn to the spin-strain coupling parameters.

The spin-strain coupling in the ground state of the NV center is more complicated than spin-electric field coupling as strain is a tensor and can result in more coupling coefficients than the electric field does. Indeed, the correct spin-strain coupling coefficients have been only recently derived, which contains six independent real coupling-strength parameters, h₄₁, h₄₃, h₂₅, h₂₆, h₁₅, and h₁₆, and has the following form [237]:

(29a)Hε=Hε0+Hε1+Hε2,

(29b)Hε0/h=[h41(εxx+εyy)+h43εzz]Sz2,

(29c)Hε1/h=12[h26εzx−12h25(εxx−εyy)]{Sx, Sz}+12(h26εyz+h25εxy){Sy, Sz},

(29d)Hε2/h=12[h16εzx−12h15(εxx−εyy)](Sy2−Sx2)+12(h16εyz+h15εxy){Sx, Sy},

where ε_ij=(∂u_i/∂x_j+∂u_j/∂x_i)/2 denotes the strain tensor and u(r) is the displacement field. Similar to Eq. (26), subscripts 0, 1, and 2 here refer to the difference in the electron spin quantum numbers ms connected by the corresponding Hamiltonian. The spin-stress Hamiltonian H_σ is analogous to Eq. (29a), with substitutions ε↦σ and h↦g. The conversion from h to g can be calculated using the stiffness tensor (C) of diamond, which connects the external stress to the internal strain in diamond (e.g. Ref. [237]) with the parameters C₁₁=1076 GPa, C₁₂=125 GPa, and C₄₄=576 GPa in the cubic reference frame.

The six spin-strain coupling-strength parameters of Eq. (29a) were determined using numerical DFT calculations [237]. To model the structure subject to mechanical strain, described by the strain tensor ε , the cubic supercell was deformed to a parallelepiped, whose edge vectors are obtained by transforming the undeformed edge vectors with the matrix 1+ε in the cubic reference frame and allow the atomic positions to relax. For each strain configuration, the elements of the 3×3 ZFS matrix D were calculated using the VASP implementation by Martijn Marsman with the PAW formalism [201].

We illustrate this methodology to obtain the six spin-strain coupling-strength coefficients with the example of h₁₆. To determine h₁₆, the supercell using a strain tensor is deformed whose only nonvanishing element is ε_yz, and the D matrix is obtained from the calculation. Due to Eq. (29a), the chosen strain configuration implies that the Hamiltonian has the form

(30)H=12εyzST⋅(0h160h160h260h260)S.

This, together with the above definition of the D matrix, yields

(31)h16=2∂Dxy∂εyz|ε=0.

The numerical error of our DFT calculations was estimated using a linear fit to the data points that were taken from 11 equidistant values of ε_yz between −0.01 and +0.01 (see Ref. [237]). The final results for all the coupling parameters are summarized in Table 2.

Table 2:

Spin-strain (h) and spin-stress (g) coupling-strength parameters of the ground state as obtained from DFT [237].

Parameter (MHz/strain)		Parameter (MHz/GPa)
h43	2300±200	g43	2.4±0.2
h41	−6420±90	g41	−5.17±0.07
h25	−2600±80	g25	−2.17±0.07
h26	−2830±70	g26	−2.58±0.06
h15	−5700±200	g15	3.6±0.1
h16	−19,660±90	g16	18.98±0.09

1. Results are rounded to significant digits. The negative sign means compression in this convention.

In Table 3, we compare the numerical DFT results of Table 2 to the experimental results of Ref. [247]. In Ref. [247], four of the six independent spin-stress coupling-strength parameters of the spin-stress interaction Hamiltonian were measured. Ref. [247] defined these four spin-stress coupling-strength parameters, denoted as a₁, a₂, b, and c, in a “hybrid” representation, where the spin-stress Hamiltonian is expressed in terms of the NV frame components of the spin vector (S_x, S_y, and S_z) and the cubic frame components of the stress tensor (σ_XX, σ_XY, etc.). To be able to make a comparison between the DFT results and the experimental ones, the notations of Ref. [247] were taken and d, e, Nx, and Ny were introduced to express the spin-stress Hamiltonian H_σ in this hybrid representation:

Table 3:

Spin-stress coupling-strength parameters in the ground state: comparison of DFT [237] and experimental [247] results.

Parameters	Relation	DFT (MHz/GPa)	Experimental results (MHz/GPa)
a1	2g41+g433	−2.66±0.07	−4.4±0.2
a2	−g41+g433	2.51±0.06	3.7±0.2
b	−g15+2g1612	1.94±0.02	2.3±0.3
c	−2g15−2g1612	−2.83±0.03	−3.5±0.3
d	−g25+2g2612	−0.12±0.01	–
e	−2g25−2g2612	0.66±0.01	–

1. Parameters in the hybrid representation (a₁, a₂, etc.) are expressed in terms of the parameters in the NV frame representation (g₄₁, etc.) in the second column. The negative sign means compression in this convention.

(32a)Hσ0/h=ℳzSz2,

(32b)Hσ1/h=Nx{Sx, Sz}+Ny{Sy, Sz},

(32c)Hσ2/h=−ℳx(Sx2−Sy2)+ℳy{Sx, Sy},

where

(33a)ℳz=a1(σXX+σYY+σZZ)+2a2(σYZ+σZX+σXY),

(33b)Nx=d(2σZZ−σXX−σYY)+e(2σXY−σYZ−σZX),

(33c)Ny=3[d(σXX−σYY)+e(σYZ−σZX)],

(33d)ℳx=b(2σZZ−σXX−σYY)+c(2σXY−σYZ−σZX),

(33e)ℳy=3[b(σXX−σYY)+c(σYZ−σZX)].

The relations between the hybrid representation parameters (a₁, a₂, b, c, d, and e) and the NV frame parameters (g₄₁, etc.) are given in the first two columns of Table 3. Importantly, H_σ0 and H_σ2 are identical to the spin-stress Hamiltonian in Eqs. (1) and (2) of Ref. [247]. It was found (Table 3) that the d and e parameters are in the same order of magnitude as the other four parameters; thus, this full six-parameter description is required to deduce accurate parameters from experimental data. From the calculated spin-stress coupling coefficients in the ground state, the ODMR contrast and basic sample dependent data (such as coherence time of the NV center’s spin), the sensitivity of strain quantum sensors can be deduced based on Hahn echo measurements at a given direction of the applied stress [248]. Next, we turn to the spin-strain interaction in the optically allowed excited state.

The ³E excited state can be written in the {A₁, A₂, E_x, E_y, E₁, E₂} manifold, and the strain will affect these as [24], [25]

    (34)

where δE1a=(εxx−εyy)/2,δE2a=(εxy+εyx)/2,δE1b=(εxz+εzx)/2, and δE2b=(εyz+εzy)/2. The degenerate levels will split due to symmetry-breaking strain components and may also mix the corresponding states. The symmetry conserving strain (ε_zz) only shifts all the levels with the same energy. The effect of the electric field is very similar to that of strain in Eq. (32) and may be given as [24], [25]

    (35)

In Ref. [24], DFT PBE simulation on molecular cluster model was applied together with electric fields to estimate the d_z and d_⊥ coefficients. The effect of the electric field was separated to “ionic” effect, i.e. piezo effect, when strain appears due to the presence of electric field, and “electronic” effect, when the cloud of the electrons and the corresponding potential changes due to the applied electric fields (see Ref. [24] and Supplementary Materials in Ref. [235]). In these simulations, the change in the geometry and splitting of the degenerate orbitals as a function of the applied electric field can be used to derive the order of magnitude for d_z and d_⊥. It was found that hd_⊥≈6.6 GHz (MV/m)⁻¹ and hd_⊥≈0.6 GHz (MV/m)⁻¹ (see Appendices C and D in Ref. [24]).

As mentioned above, the electric field is intertwined with strain via the piezo effect; thus, electric field may be applied to suppress the effect of strain [cf. Eqs. (32) and (33)] and tune two individual NV centers’ ZPL energies into the same position [39], [40]. It was also found that if a symmetry-breaking strain is present and a symmetry-breaking electric field is applied to the NV center at the same time, then the optical transition will show quadratic effect as a function of the strength of the electric field with anticrossing feature approaching the zero electric field; however, linear behavior shows up at very low strain fields [24], [38]. This explains the early measurements on the electric field dependence of the ZPL lines where single NV centers showed often linear characteristics upon the applied electric fields, but also quadratic features were observed for other single NV centers [243]. It was also shown that the optical transition from the A₂ state toward the ms=±1 ground state will preserve the circular polarization at low strain fields, which was successfully used in entanglement schemes [34]. Here, we note that, as we have already discussed above, phonons can couple A₁ state into A₂. This coupling is very small at cryogenic temperature [105] but it becomes substantial at T=10 K and above as can be inferred from the deduced ISC rates from the ms=±1 ³E states toward the ¹A₁ state [104]. Thus, the aforementioned conclusion about the strain dependence of the circular polarization in the optical transition between the A₂ state and the ms=±1 ³A₂ state is valid at very low temperatures.