A multiconfigurational study of the negatively charged nitrogen-vacancy center in diamond

In the CASSCF formalism molcas_book , a many-body wave function is described as a linear combination of multiple Slater’s determinants, each of which is made of single-electron molecular orbitals. The coefficients of the Slater’s determinants are referred to as configuration interaction (CI) expansion coefficients. A CASSCF wave function is partitioned into parts from inactive orbitals with double occupancy, virtual orbitals with zero occupancy, and active orbitals with occupancy between zero and two (i.e., 0, 1, or 2). In a CASSCF calculation, for a given spin multiplicity, any possible electron configurations or correlation within the active orbital space are included, while keeping the occupancies of the inactive and virtual orbitals fixed. However, electron excitation or correlations outside the active space are not included. Both the CI coefficients and the molecular orbitals are optimized through self-consistent calculations. Therefore, the choice of the active orbitals is critical for accurate CASSCF calculations. It was shown that the accuracy of CASSCF calculations is greatly improved by including extra molecular orbitals beyond frontier orbitals in the active space molcas_book . CASSCF wave functions are described in terms of spin-free basis states that correspond to all possible configurations generating the maximum $M_{z}$ values, where $M_{z}$ is an eigenvalue of the $S_{z}$ operator (i.e., the $z$ component of the total spin $S$). The state-average is a technique to facilitate convergence of the excited-state CASSCF wave functions molcas_book .

In order to determine the number and character of orbitals to be included in the active space, we start with a qualitative analysis of the electronic structure of an NV^- center from a single-electron point of view. The NV^- center in diamond has four broken dangling bonds, as shown in Fig. 1(a): three dangling bonds of the nearest neighboring carbon atoms to the vacancy ($d_{1}$, $d_{2}$, and $d_{3}$), and the dangling bond of the nitrogen atom to the vacancy ($d_{N}$). They form four single-electron molecular orbitals such as $a_{1}^{C}=(d_{1}+d_{2}+d_{3})/3$, $a_{1}^{N}=d_{N}$, $e_{x}=(2d_{1}-d_{2}-d_{3})$, and $e_{y}=(d_{2}-d_{3})/\sqrt{2}$  Lenef1996 ; MazeNJP011 ; DohertyNJP011 . The first two orbitals transform as a function of the $A_{1}$ irreducible representation (IRRep), and the other two orbitals transform as functions of the $E$ IRRep under the C_3v point group. It is known that the $a_{1}^{N}$ orbital is deeply buried under the valence band of the diamond lattice, whereas the other three orbitals are within the band gap GaliPRB08 ; Lenef1996 ; MazeNJP011 ; DohertyNJP011 . These three states are also referred to as in-gap defect states Bockstedte2018 . Now let us count the total number of electrons in the system. A carbon vacancy within diamond leaves four electrons in four dangling bonds. One of these carbon-atoms is substituted with a nitrogen-atom that has an extra electron (as compared to a carbon atom). The defect further acquires an additional electron and becomes negatively charged, resulting in a total number of six electrons that fill the defect states in accordance with the Hund’s rules. In the spin-triplet ground-state, the nominal occupancy is as follows: the defect state, $a_{1}^{N}$, which lies in the valence band, is doubly occupied, while the remainder of the four electrons are distributed amongst the in-gap states, with $a_{1}^{C}$ being doubly occupied, and the degenerate orbitals, ($e_{x}$ $e_{y}$), being singly occupied.

Inspired by the single-electron picture, we initially perform CASSCF calculations using the active space consisting of six electrons and the four dangling bond orbitals ($a_{1}^{N}$, $a_{1}^{C}$, $e_{x}$, and $e_{y}$) for the 70-atom and 162-atom diamond clusters with $C_{3v}$ symmetry, shown in Fig. 2(a) and (b). These calculations are referred to as CASSCF(6,4) following the number of electrons and orbitals used in the active space. The excitation energies obtained via CASSCF(6,4) calculations are highly overestimated as compared to experiment and the excited-state wave functions are found to be inconsistent with physical and chemical intuitions. As a result, we expand the active space by including extra unoccupied defect-localized states. The most common practice is to identify these extra states in the virtual space of the converged CASSCF(6,4) result. However, no such defect orbitals are found in the virtual space. Therefore, we introduce a series of CASSCF calculations discussed in the Appendix (Fig. 6) in order to identify and include extra defect orbitals in the active space. With this systematic CASSCF procedure, we find two unoccupied defect orbitals with $E$ IRRep. In order to distinguish them from the dangling bond orbitals, $e_{x}$ and $e_{y}$, they are, henceforth, referred to as $e_{x}^{\prime}$ and $e_{y}^{\prime}$ (Fig. 2(a) and (b)). The $e_{x}^{\prime}$ and $e_{y}^{\prime}$ orbitals are expected to lie in the conduction band region (i.e., defect orbitals are resonant with the conduction band). With these two extra unoccupied orbitals, as well as, the four dangling bond orbitals, we form an active space consisting of six electrons and six orbitals, and carry out CASSCF(6,6) calculations for both the total spin $S=1$ and $S=0$. Furthermore, in order to achieve high accuracy and exact numerical degeneracy (up to $\sim 10$ neV) in states with $E$ symmetry, we carefully maintain the IRRep symmetry of all of the molecular orbitals and remove the surface-dominant orbitals in the self-consistent calculations.

For the ground ${}^{3}A_{2}$ state, the first-order SOC effect on the zero-field splitting vanishes and higher-order terms are negligibly small due to weak SOC. However, for the first-excited ${}^{3}E$ state, the first-order SOC effect becomes important within the subspace of degenerate states and the SOC-induced splitting turns out to be non-negligible. Therefore, for the most accurate calculation of SOC-induced splitting, we need to describe degenerate states the most accurately. In order to achieve this, state average is carried out only over the first-excited triplet pair (${}^{3}E$) of the CASSCF(6,6) wave functions. Then SOC is included in the converged CASSCF(6,6) spin-triplet wave functions within the atomic mean-field approximation Hess1996 , using the restricted active space state interaction (RASSI) method rassi implemented in OpenMolcas. For the CASSCF(6,6) energy eigenvalues and the SOC-induced zero-field splitting, OpenMolcas is used because it provides more accurate results due to purely symmetric orbitals and removal of surface-dominated orbitals (see the Appendix).

The zero-field splitting by the SSC is expected for all spin-triplet states. This feature is computed for the CASSCF(6,6) wave functions using ORCA because it is not available in OpenMolcas. The SSC is calculated as the two-electron direct SSC over the CASSCF(6,6) wave functions using first-order perturbation theory Neese2006 , as implemented in ORCA. The CASSCF(6,6) wave functions using ORCA are obtained by following the CASSCF procedure sketched in the Appendix without orbital symmetrization, SUPERSYMMETRY keyword, and removal of surface orbitals, because they are not available in ORCA. We confirm that the zero-field splitting induced by SSC is not sensitive to technical details of the calculations (i.e., the cluster size, the size of the active space and the number of roots included in the state average).

Figure 3 shows schematic level diagrams of our calculated spin-triplet and spin-singlet states for the two cluster sizes using OpenMolcas (quantum-chemistry methods, CASSCF(6,6)). Note that we use the ground-state geometry without phonon modes and that we do not consider structural relaxation of the electronic excited states. An experimental absorption spectrum of an NV^- center in diamond consists of a sharp zero-phonon line with a broad spectrum of phonon side bands with several peaks DaviesExpt76 ; Kehayias2013 . With significant electron-phonon coupling, a zero-phonon absorption energy can noticeably differ from a vertical excitation energy. The latter energy is always higher than the former energy. The latter energy is commonly experimentally obtained from the maximum-intensity peak of the broad phonon side-band spectrum. The broadness of the phonon side bands provides some uncertainty in the maximum-intensity peak energy, which renders uncertainty in the experimental vertical excitation energy. For comparison to experiment, we provide both experimental zero-phonon absorption energies and experimental vertical excitation energies in Fig. 3.

Our calculations show that the first-excited spin-triplet ${}^{3}E$ state is separated from the ground state (${}^{3}A_{2}$) by 1.93 and 2.14 eV for the 70-atom and 162-atom clusters, respectively. This energy separation does not depend much on the cluster size and it is close to the experimental energies of zero-phonon absorption, 1.945 eV, and of vertical excitation, 2.18 eV DaviesExpt76 . We find that the lowest-energy singlet state has character of ${}^{1}E$ and that the first-excited singlet ${}^{1}A_{1}$ state is located at 1.07 eV and 1.35 eV above the ${}^{1}E$ state for the 70-atom and 162-atom clusters, respectively. The ordering and the character of the singlet states agree with experiment, considering the experimental energies of zero-phonon absorption, 1.190 eV RogersNJP08 , and of vertical excitation, 1.26 eV Kehayias2013 . Our results also reveal the energy differences between the triplet and singlet states. The ${}^{1}E$ state lies at 0.34 eV and 0.25 eV above the ${}^{3}A_{2}$ state for the 70-atom and 162-atom clusters, respectively. As a result, the energy gap between the ${}^{3}E$ and ${}^{1}A_{1}$ states becomes 0.52 and 0.54 eV for the 70-atom and 162-atom clusters, respectively. Although the energy gap between the ${}^{3}A_{2}$ state and the ${}^{1}E$ state has not been directly experimentally measured, the separation between the ${}^{3}E$ state and the ${}^{1}A_{1}$ state was measured to be 0.321-0.414 eV Goldman2015 ; Goldman2017 , which is in good agreement with our results. The second-excited (third-excited) triplet state has characteristics of ${}^{3}A_{1}$ (${}^{3}E$). The second-excited singlet ${}^{1}A_{2}$ state appears even above the third-excited triplet ${}^{3}E$ state. There are no experimental reports on the higher-energy levels or separations.

Our calculated results show that for the four lowest states (${}^{3}A_{2}$, ${}^{3}E$, ${}^{1}E$, and ${}^{1}A_{1}$) the energy eigenvalues do not depend much on the cluster size. However, we find that the cluster-size dependence becomes more apparent for higher-energy states, especially for the second- and third-excited triplet states (${}^{3}A_{1}$ and ${}^{3}E$). Depending on the cluster size, the energy separations change but the ordering of the states does not change. A similar trend of the cluster-size dependence was reported in the complete-active space approach, using DFT Kohn-Sham orbitals and density-matrix renormalization group Barcza2020 . This trend can be understood by the fact that higher-energy levels have stronger electron correlations which requires inclusion of more empty orbitals in the active space. Since experimental data are available for mainly up to the first-excited triplet ${}^{3}E$ state, we do not further study an effect of cluster size on the electronic structure.

Let us now compare our calculated energies of the spin-triplet and spin-singlet states (${}^{3}A_{2}$, ${}^{3}E$, ${}^{1}E$, and ${}^{1}A_{1}$) to the previous ab-initio theoretical studies. See Table 1 and Fig. 4. In our analysis, we focus on the four lowest states because only the level separations among them were experimentally measured and because higher-energy states are more sensitive to the cluster size and the size of the active space. (For example, the higher-energy ${}^{1}E^{\prime}$ state that many-body theory studies predicted Bockstedte2018 ; Ma2010 ; ChoiPRB012 has not been observed Kehayias2013 .) We first discuss comparison to other DFT calculations and then to other quantum-chemistry studies as well as many-body theory studies, separately.

Earlier DFT studies of an NV^- center in diamond clusters and periodic supercells GaliPRB08 ; Goss1996 ; DelaneyNL10 ; LinJCP08 ; LarsonPRB08 ; Galireview2019 showed that the calculated excitation energy of the ${}^{3}E$ state more or less agrees with our result and experiment except for Ref. DelaneyNL10 . However, DFT-calculated energies of the singlet states are scattered in a wide range and the ordering of the triplet and singlet states is inconsistent with recent experiment. This trend is understandable considering that DFT poorly describes the singlet states due to the well-known spin contamination effect.

In the previous quantum-chemistry studies of an NV^- center in diamond clusters DelaneyNL10 ; ZyubinJPC07 ; LinJCP08 , either the excitation energies are significantly different from experiment, or the ordering of the singlet and triplet states is reversed. More specifically, CASSCF(6,8) calculations discussed in Ref. LinJCP08 showed that the excitation energy of the ³E state is 0.5-0.6 eV (0.3-0.4 eV) higher than our result (experiment). The singlet states were not investigated in that work. In the CASSCF(8,11) calculations presented in Ref. ZyubinJPC07 , the excitation energy of the ³E state is about 1.0 eV lower than our result or experiment, and the singlet ${}^{1}A_{1}$ state is slightly above the triplet ${}^{3}E$ state, which does not agree with our result or recent experiment. Their multireference configuration interaction (MRCI) calculations ZyubinJPC07 somewhat increase the energies of the triplet and singlet states with the correct ordering of the excited triplet and singlet states. However, the energy of the ${}^{3}E$ state remains lower than our value by about 0.6 eV. Monte Carlo configuration interaction (MCCI) studies DelaneyNL10 showed the energy of the ${}^{3}E$ state in agreement with our result and experiment. However, the ordering of the ${}^{3}E$ and ${}^{1}A_{1}$ states is reversed. See Table 1 and Fig. 4. The discrepancies between our results and all of the earlier quantum-chemistry calculations arise from the choice of orbitals in the active space. One of the most common ways to choose active orbitals is to use single-electron molecular orbitals in the vicinity of the band gap such as orbitals near the highest occupied molecular orbitals (HOMO) and lowest unoccupied molecular orbitals (LUMO). For a hydrogen-passivated diamond cluster with an NV^- center, either this common practice within CASSCF or MRCI, or an automatic choice of the active space in MCCI may result in non-physical surface-dominated orbitals in the CI basis set. As shown in Fig. 2(a), in our case, all six orbitals in the active are localized near the vacancy defect.

An earlier many-body perturbation study Ma2010 based on the $GW$ approximation with Bethe-Salpeter equation (BSE) provided the singlet-singlet (${}^{1}A_{1}$ - ${}^{1}E$) energy difference about 0.6 eV lower than our result and recent experiment Kehayias2013 ; Goldman2015 , although the energy of the ${}^{1}E$ state relative to the ground state, as well as the ordering of the two singlet states are in agreement with the recent experimental data. On the other hand, recent many-body studies Bockstedte2018 ; Ma2020 ; ChoiPRB012 showed more promising results by using effective many-body model Hamiltonians with parameters obtained from (or fitted to) ab-initio calculations in order to properly include many-body character in the wave functions. For example, additional unoccupied defect states (resonant to the conduction band) and doubly occupied defect states (in the valence band) were included in the configuration interaction constrained random phase approximation (CI-CRPA) method Bockstedte2018 . This is analogous to our inclusion of unoccupied level 5 and 6 and doubly occupied level 1 [Fig. 2(b)] in the active space for proper treatment of electron correlation. Their results are closest to our result among the previous studies that we have discussed (see Fig. 4 and Table 1). Yet, there are some differences. In the fitting of $GW$-calculated bands to model Hamiltonian ChoiPRB012 (in the CI-CRPA method Bockstedte2018 ), the singlet-singlet energy difference is about 0.2-0.3 eV (0.3-0.4 eV) lower than our result and experiment. In the beyond-RPA implemented in the quantum embedding theory Ma2020 , the energy difference between the ${}^{3}E$ and ${}^{1}A_{1}$ states is somewhat smaller than our result and experiment. This discrepancy may arise from missing orbital configurations in the ${}^{1}E$ and ${}^{1}A_{1}$ states in Refs. Bockstedte2018 ; Ma2020 ; ChoiPRB012 that are discussed in Sec.V.3. Here we stress that it does not seem to be straightforward to include effects of SOC and SSC within the formalisms used in Refs. Bockstedte2018 ; Ma2020 ; ChoiPRB012 in contrast to the quantum chemistry methods where such effects can be added to the many-body wave functions without an introduction of new fitting parameters (see Sec.V.4).

We now discuss characteristics of our calculated triplet and singlet energy eigenstates (Table 2). Here we use configuration basis states which are all possible states generating the maximum $M_{z}$ value from the six active orbitals for a given total spin $S$, where $M_{z}$ is an eigenvalue of the $S_{z}$ operator. The total wave functions in terms of true $S_{z}$ eigenstates are obtained when SOC is applied to the many-body (CASSCF) wave functions within the RASSI method rassi using the Wigner-Eckart theorem. The SOC effect is discussed later in Sec. V.4.

For the ground and first-excited triplet states, the configurations of our calculated eigenstates are similar to those identified from the phenomenological molecular models based on group theory MazeNJP011 ; DohertyNJP011 ; Tamarat2008 , as long as we focus on the configurations with weights greater than 10%. However, for the singlet states, we find that the following additional configurations significantly contribute: ${\tiny\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{$~{}~{}\uparrow{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{$~{}~{}\downarrow{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}}$ and ${\tiny\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{$~{}~{}\uparrow{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{$~{}~{}\downarrow{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}}$ with 12% each for the ${}^{1}E$ state and ${\tiny\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{$\uparrow\downarrow$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}\framebox[9.03374pt]{${\color[rgb]{1,1,1}\uparrow}{\color[rgb]{1,1,1}\uparrow}$}}$ with 20% for the ${}^{1}A_{1}$ state. Refer to Table 2 for the notations. The former states indicate single-excitations from the doubly occupied $a_{1}^{C}$ level, while the latter state indicates a double-excitation from the $a_{1}^{C}$ level. These configurations have not been considered before in the literature. Their inclusion in our work may have given rise to the discrepancy between our result and those obtained in Refs. Bockstedte2018 ; Ma2020 ; ChoiPRB012 and it may also affect the intersystem crossing.

Furthermore, above the first-excited triplet ${}^{3}E$ state and the first-excited ${}^{1}A_{1}$ state, we find the triplet ${}^{3}A_{1}$ and ${}^{3}E$ states and the singlet ${}^{1}A_{2}$ state. Due to the lack of experimental data beyond the four lowest states, we only briefly mention these higher-energy states. Our higher-energy states differ from those in the literature Bockstedte2018 ; ZyubinJPC07 ; Ma2010 ; ChoiPRB012 ; MazeNJP011 ; DohertyNJP011 . As shown in Table 2, the main contributions to these states originate from single excitations from the $a_{1}^{C}$, $e_{x}$, or $e_{y}$ orbital to beyond the dangling bond orbitals ($e_{x}^{\prime}$ and $e_{y}^{\prime}$) (Fig. 2(b)). On the other hand, the previous many-body and molecular-model studies Ma2010 ; ChoiPRB012 ; MazeNJP011 ; DohertyNJP011 were mostly obtained considering only three or four dangling bond orbitals ($a_{1}^{N}$, $a_{1}^{C}$, $e_{x}$, and $e_{y}$). As discussed earlier, the higher-energy states are more sensitive to the size of active space and cluster size than the four lowest states due to stronger electron correlation. Note that the ${}^{1}E^{\prime}$ state predicted in the literature has not been experimentally observed Kehayias2013 .

All of the spin-triplet states that we discussed earlier are split due to SOC and/or SSC. Note that SOC plays an important role in the zero-field splitting only for the degenerate levels in this system because of the weak SOC. Since experimental data do not exist for higher-energy states, we present calculated zero-field splitting values of the ground ${}^{3}A_{2}$ state and the first-excited triplet ${}^{3}E$ state only. Figure 5 and Table 3 show our calculated level splitting by SOC alone and by SOC in combination with SSC (SOC+SSC) for the ${}^{3}A_{2}$ and ${}^{3}E$ states, separately, compared to experimental data BatalovPRL90 . Table 4 lists the corresponding eigenvectors $\Psi_{1,...,9}$ obtained from the quantum chemistry calculations including SOC and SSC. The SOC-induced level splitting is obtained for the 70-atom and 162-atom clusters, while the SSC-induced splitting is obtained for the 70-atom cluster. Regarding the SOC+SSC induced splitting for the 162-atom cluster, we use the SSC-induced splitting for the 70-atom cluster since the SSC-induced splitting does not depend much on cluster size. Let us now discuss the ${}^{3}A_{2}$ and ${}^{3}E$ states separately.

The SOC does not split the ${}^{3}A_{2}$ state to the first order and its splitting by higher-order SOC is negligible. However, we find that the SSC splits the ${}^{3}A_{2}$ state into one lower non-degenerate level with $M_{z}=0$ and one higher doubly degenerate level with $M_{z}=\pm 1$ by $-$1.9 GHz and 0.8 GHz, respectively. (See the eigenvectors $\Psi_{1,2,3}$ in Table 4.) Therefore, the energy separation between them is about 2.7 GHz, which is in excellent agreement with the experimental value of 2.88 GHz BatalovPRL90 as well as a previous DFT calculation Viktor014 .

On the other hand, the SOC splits the ${}^{3}E$ state into three (degenerate) groups, each of which has eigenvalues of the $z$ component of orbital angular momentum, $L_{z}$, of $\pm$0.46 ($\pm$0.53) for the 70-atom (162-atom) cluster. The separation of the levels is about 6.5 GHz for the 70-atom cluster (Fig. 5) and about 8.1 GHz for the 162-atom cluster. Our calculated level splitting values show a weak cluster-size dependence and they are somewhat larger than the experimental value of 5.3 GHz BatalovPRL90 . A possible reason for this is the dynamic Jahn-Teller effect Fu2009 ; Abtew2011 and the resulting quenching of SOC (i.e., Ham reduction factor Ham1968 ; Thierling2017 ; Goldman2017 ) Note that our calculations are done for zero strain without electron-phonon coupling. Quantum-chemistry calculations of electron-phonon coupling and the dynamic Jahn-Teller effect are worth investigating in the future. In addition to the SOC-induced splitting, the SSC further shifts the lowest degenerate level upward by 0.8 GHz ($\Psi_{4}$, $\Psi_{5}$ in Table 4) and moves the second degenerate level downward by 1.9 GHz ($\Psi_{6}$, $\Psi_{7}$). In this case, the degeneracy still holds. Interestingly, the amount of the downward level shift is almost twice that of the upward shift. The trend of the level-shift direction as well as the ratio between the downward and upward shift amount, are in good agreement with experiment BatalovPRL90 , although our shifted values are off by a factor of 2 compared to experiment. We also find that the SSC splits the third doubly degenerate level into two separate levels ($\Psi_{8}$, $\Psi_{9}$): one level shifts downward by 3.1 GHz and the other moves upward by 5.3 GHz. Again, the trend of the level shift agrees with experiment BatalovPRL90 , although the calculated shift amount is greater than experiment by a factor of 2 or 3. This overestimated SSC contribution may partially arise from our first-order perturbation treatment of SSC.