Electronic excitations of the charged nitrogen-vacancy center in diamond
obtained using time-independent variational density functional calculations

Aleksei V. Ivanov1$\star$, Yorick L. A. Schmerwitz2, Gianluca Levi2 and Hannes
Jónsson3$\star$

1 Riverlane, St Andrews House, 59 St Andrews Street, Cambridge, CB2 3BZ,
United Kingdom

2 Science Institute of the University of Iceland, VR-III, 107 Reykjavík,
Iceland

3 Faculty of Physical Sciences, University of Iceland, VR-III, 107 Reykjavík,
Iceland

⋆<EMAIL_ADDRESS><EMAIL_ADDRESS>
## Abstract

Elucidation of the mechanism for optical spin initialization of point defects
in solids in the context of quantum applications requires an accurate
description of the excited electronic states involved. While variational
density functional calculations have been successful in describing the ground
state of a great variety of systems, doubts have been expressed in the
literature regarding the ability of such calculations to describe electronic
excitations of point defects. A direct orbital optimization method is used
here to perform time-independent, variational density functional calculations
of a prototypical defect, the negatively charged nitrogen-vacancy center in
diamond. The calculations include up to 511 atoms subject to periodic boundary
conditions and the excited state calculations require similar computational
effort as ground state calculations. Contrary to some previous reports, the
use of local and semilocal density functionals gives the correct ordering of
the low-lying triplet and singlet states, namely
${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$. Furthermore, the more advanced meta
generalized gradient approximation functionals give results that are in
remarkably good agreement with high-level, many-body calculations as well as
available experimental estimates, even for the excited singlet state which is
often referred to as having multireference character. The lowering of the
energy in the triplet excited state as the atom coordinates are optimized in
accordance with analytical forces is also close to the experimental estimate
and the resulting zero-phonon line triplet excitation energy is underestimated
by only 0.15 eV. The approach used here is found to be a promising tool for
studying electronic excitations of point defects in, for example, systems
relevant for quantum technologies.

###### Contents

  1. 1 Introduction
  2. 2 Model and computational method
  3. 3 Results
  4. 4 Discussion and conclusions
  5. A Relation between multideterminant states and single determinants
  6. B Vertical excitation and zero-phonon line
  7. C Smaller cell calculations
  8. D Supplemental information

##  1 Introduction

The negatively charged nitrogen-vacancy center (NV- center) in diamond
exhibits remarkable optical and magnetic properties, making it a promising
candidate for quantum applications, including nanoscale sensing [1, 2, 3, 4],
quantum communication via single photon emission [5, 6, 7] and quantum bits
(qubits) for quantum computing [8, 9, 10, 11, 12]. The applicability of the
NV- center in quantum technologies derives from the possibility of generating
a pure spin state with long coherence time through optical excitation. In
order to optimally control the process and establish a theoretical approach
that can guide the search of other systems for which a pure spin state can be
prepared, accurate modelling of the electronic defect levels and corresponding
excited states is desired.

Early theoretical studies of the NV- center in diamond using density
functional theory (DFT) calculations with local and semilocal Kohn-Sham
functionals[13, 14, 15] have led to contradictory results in that they do not
agree on the ordering of the low-lying triplet and singlet defect states. In
particular, there has been disagreement on which of two singlet states,
${}^{1}E$ or ${}^{1}A_{1}$, is lower in energy and whether or not both
singlets lie between the two lowest triplet states. These conflicting results
have contributed to a long-standing controversy regarding the electronic
states involved in the optical spin-polarization cycle of the NV- center[16,
17, 18]. Only recently has this been resolved thanks to experimental
observations and high-level many-body calculations[19, 20, 21, 22, 23, 24, 25,
26, 27] where the ordering

${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$

of the energy levels has been established. Spin initialization in the NV-
center is realized through optical excitation from the triplet ground state
(${}^{3}A_{2}$) to an excited triplet state (${}^{3}E$), the system then
crossing over to an excited singlet state (${}^{1}A_{1}$), followed by de-
excitation to the lowest singlet state (${}^{1}E$) and finally returning to
the ground triplet state. Optical excitation occurs from both the $m_{s}=0$
and $m_{s}=\pm 1$ sublevels of the triplet ground state, but electronic
relaxation preferentially populates the $m_{s}=0$ state. Therefore, the NV-
center can be initialized in the $m_{s}=0$ state after a few such optical
cycles. The electrons that participate in the electronic transitions of the
optical cycle occupy single-particle states between the valence and conduction
bands and are thus localized at the defect.

As a result of the contradictory calculations mentioned above, there are
several statements in the literature to the effect that the density functional
approach cannot adequately describe the electronic structure of the NV- center
in diamond and that calculations of the singlet electronic states require a
higher level of theory because of multiconfigurational character.

Several different electronic structure methods are typically used for
calculations of band gaps and defect levels in semiconductors for quantum
technologies. A review is given in Ref.[28]. One of the most common methods is
many-body perturbation theory based on the GW approximation [29], which uses
the Green function formalism, together with the Bethe-Salpeter equation (BSE)
[30, 31]. However, it has been shown that the GW+BSE method does not provide a
satisfactory description of the excited singlet ${}^{1}A_{1}$ state of the NV-
center in diamond [16, 26]. This has been attributed to the fact that the wave
function of this state includes contributions of double excitations while
GW+BSE only takes into account singly excited configurations [32]. Moreover,
this approach involves large computational effort, too large for the
supercells needed to accurately describe isolated defect states in
semiconductors.

Time-dependent density functional theory (TDDFT) [33, 34] is another commonly
used method for calculating excited electronic states. Most TDDFT studies of
the excited states of the NV- center in diamond have been carried out within
the adiabatic and linear-response approximations and describe the system with
a molecular cluster model where the surface atoms are saturated with bonds to
hydrogen atoms [35, 13, 27, 36]. Such finite models do not accurately describe
the band structure and defect levels of semiconductors because of quantum
confinement effects and the admixing of artificial surface states [32, 28].
Recently, Galli and co-workers[20] have performed spin-flip TDDFT calculations
of the NV- center in diamond using a supercell approach including periodic
boundary conditions. Using the PBE functional as well as a hybrid functional
that includes exact exchange, the right ordering of the states is obtained,
but PBE gives too low energy for the excited singlet state, and the hybrid
functional gives too high energy for both the singlet and triplet excited
states, which was attributed to the lack of doubly excited configurations in
TDDFT.

Spin defects in semiconductors have also been modelled using quantum embedding
methods [37, 38, 22, 21, 24]. Here, the defect states are included in an
active space described with a many-body effective Hamiltonian, and the
interaction with the environment is taken into account through dielectric
screening evaluated using DFT. This approach is found to estimate accurately
the energy of the vertical excitations. However, quantum embedding
calculations depend on the size of the active space, the method used to avoid
double counting Coulomb interactions, and the procedure for obtaining the bulk
screening, in addition to the choice of the energy functional [38].
Furthermore, since atomic forces are at present not available from quantum
embedding calculations, the approach has so far relied on DFT calculations to
estimate the effect of changes in the atomic coordinates in the excited states
to obtain, for example, the zero-phonon line (ZPL) excitation energy.

While DFT is a ground state theory, various time-independent generalizations
for calculating excited states exist [39, 40, 41, 42, 43]. In practical
calculations, excited states can be found as higher energy stationary points
on the surface describing how the energy varies as a function of the
electronic degrees of freedom. Typically, these stationary solutions
correspond to saddle points, making it necessary to employ an optimization
method that avoids collapse to the ground state [44, 45, 46, 47, 48, 49, 50,
51, 52]. This approach 111In the literature, variational density functional
calculations of excited states are also referred to as delta self-consistent
field ($\Delta$SCF) calculations involves similar computational effort as
ground state calculations and has proven useful for modelling excited states
of both extended [46, 53, 54] and finite systems [55, 45, 56, 47, 57, 51, 58,
59, 60]. Since the excited states are obtained as stationary solutions of a
set of single-particle equations, such as the Kohn-Sham equations [61], the
orbitals are variationally optimized making it possible in principle to
evaluate atomic forces analytically, thereby opening the possibility of
performing atomic structure optimization and simulation of dynamics in the
excited state [55, 57, 62, 58]. Although only a single Slater determinant is
optimized, complex potential energy surfaces for atomic motion have been shown
to be described accurately, including avoided crossings and conical
intersections (where TDDFT calculations are usually problematic), even for
atomic configurations that are typically treated with a multiconfigurational
approach [63, 62, 64]. For example, a conical intersection and avoided
crossing in the ethylene molecule has been shown to be described well when
symmetry breaking of the wave function is allowed for [63]. This is analogous
to calculations of ground states that are inherently multiconfigurational
(sometimes referred to as "strongly correlated") where symmetry breaking of
the wave function gives improved estimate of the energy [65, 66], the
stretched H2 molecule being the classic example.

As mentioned above, DFT calculations of the excited states of the NV- center
in diamond have given contradictory results. On the one hand, early
calculations of Goss et al.[15] using a molecular cluster model and the local
density approximation (LDA) functional [67] give the right ordering of the
energy levels, ${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$. On the other hand,
the calculations of Delaney et al. [13] using larger molecular clusters and
the Becke-Perdew (BP) exchange-correlation functional [68, 69] predict the
${}^{1}A_{1}$ singlet excited state to be higher in energy than the ${}^{3}E$
triplet excited state, and Gali et al. [14] reported calculations based on
periodic boundary conditions and LDA where the singlet ${}^{1}A_{1}$ state is
lower in energy than the singlet ${}^{1}E$ state, approximately isoenergetic
with the ${}^{3}A_{2}$ triplet ground state. This disagreement between
previous DFT calculations of the excited states of the NV- center in diamond
is often taken as an indication of the inability of the method to describe
multiconfigurational ("strongly correlated") states [20, 22, 24, 16]. However,
single-determinant mean-field approaches have in many cases been shown to give
quite good approximations to the energy of multiconfigurational systems, for
example in bond-breaking processes and near avoided crossings, as well as for
challenging molecules such as diradicals and the carbon dimer [65, 63, 66, 70,
71].

Here, the states of the NV- center in diamond are calculated using a recently
developed direct orbital optimization method for a periodic supercell
representation of the system and plane-wave basis set[45]. The functionals
used include LDA, the generalized gradient approximation (GGA) PBE functional
[72], and two meta-GGA functionals, TPSS [73] and r2SCAN [74]. The
calculations using any one of these functionals are found to give the correct
ordering of the states, namely ${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$, with
the meta-GGA functionals providing vertical excitation energy values that are
in remarkably close agreement with results of advanced quantum embedding
calculations [21]. The relaxation energy in the triplet excited state and the
ZPL energy of the optical transition from the triplet ground state are also
found to be in good agreement with experimental estimates. The results
presented here based on variational density functional calculations of the
excited states demonstrate that this approach can indeed give accurate results
even for open shell singlet states that are typically considered to have
multireference character, and thereby provides a useful tool for modeling
excited states of point defects in semiconductors with much smaller
computational effort than the various higher level approaches.

##  2 Model and computational method

The NV- center in diamond consists of a substitutional nitrogen atom and a
nearest-neighbour carbon vacancy and possesses trigonal C${}_{3\text{v}}$
symmetry (see Figure 1). The low-lying triplet and singlet states can be
described with three orbitals: a lower-energy $a_{1}$ orbital and a pair of
higher-energy, degenerate $e$ orbitals, $e_{x}$ and $e_{y}$, that are
localized on the carbon atoms around the vacancy. These orbitals are occupied
by four electrons as shown in Figure 1. The ${}^{3}A_{2}$ triplet ground state
can be represented with the $m_{s}=1$ single Slater determinant
$\ket{e_{y}e_{x}}$, hereafter denoted ${}^{3}\Phi_{1}$, where the $a_{1}$
orbital is doubly occupied and the $e_{x}$ and $e_{y}$ orbitals are singly
occupied with spin up electrons. The ${}^{3}E$ triplet excited state is
obtained by promotion of an electron from the $a_{1}$ orbital to one of the
doubly degenerate $e$ orbitals and can be represented with the single
determinant $\ket{a_{1}e_{x}}$, hereafter denoted ${}^{3}\Phi_{4}$.

The two singlet states, ${}^{1}E$ and ${}^{1}A_{1}$, have multideterminant
character and need to be represented by two or more Slater determinants. The
symmetry-adapted wave functions of these states are [75, 15]

$\displaystyle\Psi\left({}^{1}E\right)=\begin{cases}&\left(\ket{e_{x}\overline{e_{y}}}+\ket{e_{y}\overline{e_{x}}}\right)/\sqrt{2},\\\
&\left(\ket{e_{x}\overline{e_{x}}}-\ket{e_{y}\overline{e_{y}}}\right)/\sqrt{2}.\\\
\end{cases}$ (1)
$\displaystyle\Psi\left({}^{1}A_{1}\right)=\left(\ket{e_{x}\overline{e_{x}}}+\ket{e_{y}\overline{e_{y}}}\right)/\sqrt{2}.$
(2)

Figure 1: Atomic configuration of the NV- center in diamond and representation
of the orbitals corresponding to the defect levels lying within the band gap.
C atoms: grey; N atom: dark grey; vacancy site: orange. The orbitals are
obtained from PBE functional calculations of each one of the Slater
determinants and are rendered at an isovalue of 0.25 Å-3/2. The occupation of
the orbitals in the determinants used to obtain the energy of the low-lying
triplet and singlet states is also indicated. The ground and excited triplet
states, ${}^{3}A_{2}$ and ${}^{3}E$, are calculated as the $m_{s}=1$ single
Slater determinants ${}^{3}\Phi_{1}$ and ${}^{3}\Phi_{4}$, respectively. The
singlet states, ${}^{1}E$ and ${}^{1}A_{1}$, need to be represented by two or
more Slater determinants and their energy is obtained from calculations of the
ground triplet state determinant, ${}^{3}\Phi_{1}$, the spin-symmetry-broken
determinant, ${}^{m}\Phi_{2}$, and the doubly excited spatial-symmetry-broken
determinant, ${}^{1}\Phi_{3}$, according to Eqs. 5 and 6.

By introducing the following linear combinations of the $e_{x}$ and $e_{y}$
orbitals

$\displaystyle e_{-}=\frac{e_{x}-e_{y}}{\sqrt{2}},$ (3) $\displaystyle
e_{+}=\frac{e_{x}+e_{y}}{\sqrt{2}},$ (4)

the energy of the multideterminant singlet states can be calculated from the
energy of single determinants using the formulas (see Appendix A for a
derivation)

$\displaystyle\mathcal{E}[{}^{1}E]=2\mathcal{E}[{}^{m}\Phi_{2}]-\mathcal{E}[{}^{3}\Phi_{1}],$
(5)
$\displaystyle\mathcal{E}[{}^{1}A_{1}]=\mathcal{E}[{}^{3}\Phi_{1}]+2(\mathcal{E}[{}^{1}\Phi_{3}]-\mathcal{E}[{}^{m}\Phi_{2}]),$
(6)

where ${}^{m}\Phi_{2}=\ket{e_{-}\overline{e_{+}}}$ is a determinant with
broken spin symmetry and ${}^{1}\Phi_{3}=\ket{e_{+}\overline{e_{+}}}$ is a
doubly excited determinant with broken spatial symmetry (see Figure 1). Eq.
(5) represents spin-purification [76, 77].

Calculations are carried out using the following density functionals: LDA
[67], PBE [72], TPSS [73], and r2SCAN [74]. The orbitals are represented with
a plane-wave basis set and the projector augumented wave method [78]. The
plane-wave basis corresponds to a 600 eV kinetic energy cutoff. First, the
lattice parameters of a 512-atom supercell of diamond are optimized using the
PBE functional. Then, the NV- defect center is introduced and the resulting
structure optimized in the ground triplet state represented by the
${}^{3}\Phi_{1}$ determinant for each of the chosen functionals until the
largest atomic force is below 0.01 eV/Å. The energy of vertical excitations
for both singlet and triplet states is obtained at the geometry optimized in
the triplet ground state. The excited state determinants, ${}^{1}\Phi_{3}$ and
${}^{3}\Phi_{4}$, correspond to saddle points on the electronic energy
surface. All determinants have been calculated using the direct orbital
optimization method presented in Ref. [45] which makes use of a limited-memory
version of the symmetric rank-one (L-SR1) quasi-Newton algorithm[47] to assist
the convergence on the saddle points. Calculations use integer occupation
numbers and no symmetry constraints are enforced. The ZPL energy for the
optical transition between the two triplet states is obtained by optimizing
the atomic structure in the excited triplet state and then evaluating the
energy difference with respect to the relaxed ground state. All calculations
are performed with the GPAW [79], Libxc [80] and ASE [81] software.
Visualization of orbitals has been carried out with the VMD software [82]. In
order to ensure the 511-atom supercell (with the vacancy) is large enough,
calculations were also carried out with a smaller cell containing 215 atoms,
and the results were found to differ by at most 5 meV (see Tables 2, and 3 in
the Appendix). Given this small difference, it is appropriate to compare the
results obtained here with a 511-atom supercell with the previous periodic
calculations where a 215-atom supercell was used.

##  3 Results

Figure 2 shows results obtained here with the LDA functional (leftmost
column). The right ordering of the energy levels is obtained,
${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$.

Figure 2: Energy of vertical excitations of the NV- center in diamond from the
triplet ground state. The leftmost column shows results obtained here from
variational calculations with the LDA functional. The red horizontal shading
indicates a range of $\pm 0.1$ eV around those values. The ordering of the
states, ${}^{3}A_{2}<{}^{1}E<{}^{1}A_{1}<{}^{3}E$, is in agreement with the
best estimates. For comparison, published results of LDA calculations [14],
also using periodic boundary conditions, are shown in the next column. There,
the ordering of states is different for some unknown reason. The third column
shows results of molecular cluster LDA calculations [15]. The last column
shows results of calculations [13] using the gradient-dependent BP functional,
where the ordering of the excited triplet and excited singlet states is
reversed. aRef.[14], bRef.[15], cRef.[13].

The figure also shows a comparison with two previously reported LDA
calculations (second and third column), one based on periodic boundary
conditions as in the present study [14], and another for a molecular cluster
[15]. The molecular cluster LDA calculations of Goss et al. [15] also get the
right ordering of the energy levels, whereas Gali et al. [14] report a singlet
${}^{1}A_{1}$ energy close to that of the ${}^{3}A_{2}$ triplet ground state,
which is in strong contradiction with high level calculations as well as
experimental measurements. It is unclear what the reason for this discrepancy
is. The figure also shows results obtained with the BP gradient dependent
functional [13] (fourth column) and a cluster model. There, the ordering of
the excited triplet and excited singlet states is reversed,
${}^{3}E<{}^{1}A_{1}$.

Figure 3: Energy of vertical excitations relative to the triplet ground state
of the NV- center in diamond obtained with variational calculations using
different local and semilocal density functional approximations, and
comparison with results of previous calculations based on many-body
approaches: periodic quantum embedding beyond the random phase approximation,
beyond-RPA (bRPA)[21], constrained RPA (cRPA)[24], extended Hubbard model
fitted to GW calculations[16], periodic GW+BSE[26], and molecular cluster
complete active space self-consistent field (CASSCF)[23] calculations. The red
horizontal shadings span $\pm 0.075$ eV around the values obtained with the
beyond-RPA quantum embedding results [21], which are taken to give the best
theoretical estimates. The r2SCAN functional gives results that are remarkably
close, the largest deviation being below 0.06 eV. All four density functionals
used here give the correct ordering of the energy levels of the electronic
states. aRef.[21], bRef.[24], cRef.[16], dRef.[26] eRef.[23].

Figure 3 shows results obtained here with GGA and meta-GGA functionals, in
addition to the LDA results, as well as a comparison with results obtained
using various many-body calculations. The numerical values are listed in Table
2 in the Appendix. The correct ordering of the energy levels is obtained for
all of the functionals. The largest absolute changes in excitation energy with
respect to the density functional are obtained for the excited singlet state,
${}^{1}A_{1}$, but the relative changes of the two singlet states are similar.
The excited triplet state, which is the only state apart from the ground state
that can be described using a single determinant, is affected less by the
choice of functional. The values of the excitation energy increases as the
complexity of the functional increases, in the order LDA $<$ GGA $<$ meta-GGA,
with the recently developed r2SCAN functional also giving slightly larger
values for the energy of the vertical excitations than the TPSS functional.

The most accurate values are believed to be the results of Ma et al. [21]
shown in the fifth column from the left in Figure 3. These are obtained using
quantum embedding calculations beyond the random phase approximation (from now
on referred to as ’beyond-RPA’, but labeled bRPA in the figure) and include
explicit exchange-correlation effects of the environment on the defect energy
levels. Column six shows results of periodic quantum embedding calculations
where the screened Coulomb interactions are evaluated using a constrained
random phase approximation (cRPA), neglecting exchange-correlation effects
[24]. There, the energy of the excited singlet state is significantly lower
than in the beyond-RPA calculation.

The remaining columns in figure 3 show results of an extended Hubbard model
fitted to GW calculations (Hubb.@GW) [16], periodic GW+BSE calculations [26],
and recent molecular cluster complete active space self-consistent field
(CASSCF) calculations with state averaging [23]. In all cases the ordering of
the states is the same, but the two singlet states differ the most, consistent
with the fact that they have multireference character. The GW+BSE calculation
clearly gives a too large value for the excited triplet state.

It is clear from the results shown in figure 3 that the r2SCAN functional
provides the best results of the four functionals used in the present study.
Remarkably close agreement is obtained with the most accurate results coming
from the beyond-RPA calculations. The deviation of the values obtained with
r2SCAN from those of the beyond-RPA values is below 0.06 eV for all the
excited states (see Table 1).

| ${}^{3}A_{2}\leftrightarrow{}^{3}E$ | ${}^{3}A_{2}\leftrightarrow{}^{1}A_{1}$ | ${}^{3}A_{2}\leftrightarrow{}^{1}E$ | ${}^{1}E\leftrightarrow{}^{1}A_{1}$ | ${}^{1}A_{1}\leftrightarrow{}^{3}E$
---|---|---|---|---|---
r2SCAN | 2.057 | 1.804 | 0.621 | 1.183 | 0.253
beyond-RPA (Ref. [21]) | 2.001 | 1.759 | 0.561 | 1.198 | 0.243
Difference | 0.056 | 0.045 | 0.060 | -0.015 | 0.010
Table 1: Calculated energy of vertical excitations (in eV) between states of
the NV- center in diamond obtained with the r2SCAN functional, and taken from
reported theoretical best estimates using a beyond-RPA method [21] as well as
the difference between the two. The r2SCAN results are remarkably close to the
beyond-RPA results.

We note that an earlier version of the SCAN functional[83] gives quite similar
results (see Table 2 in the Appendix). Calculations using the TPSS functional
also provide results in quite good agreement with the beyond-RPA calculations,
with the largest difference being in the excited singlet state, ${}^{1}A_{1}$,
where the vertical excitation energy is underestimated by $\sim$0.25 eV. A
similar deviation is obtained in the Hubbard-model and cRPA many-body
calculations. The largest deviations in the LDA and PBE calculations are also
in the energy of excitation to the ${}^{1}A_{1}$ state, which is
underestimated with respect to the beyond-RPA results by $\sim$0.70 and
$\sim$0.45 eV, respectively. The excitation energy for the ${}^{1}E$ singlet
ground state and ${}^{3}E$ triplet excited state are predicted within $\pm
0.1$ eV from the values of the beyond-RPA reference calculations for all the
density functionals used here. In comparison, the more computationally
intensive GW+BSE calculations of ref. [26] give an accuracy comparable to LDA
for the ${}^{1}A_{1}$ state and have relatively large deviation of 0.3 eV with
respect to the beyond-RPA result for the ${}^{3}E$ triplet excited state. The
CASSCF calculations of ref. [23] also have larger errors for the ${}^{3}E$ and
${}^{1}E$ states than all the density functionals employed here, including
LDA.

Figure 4: Energy of zero-phonon line (ZPL) excitations (solid lines) of the
NV- center in diamond obtained experimentally and from variational
calculations using the four density functionals, as well as the energy of
vertical excitations (dashed lines, same values as in Figure 3) to the singlet
states where energy lowering due to changes in atomic coordinates have not
been included. The green horizontal shading represents the uncertainty in the
experimental value of the ionization energy of the singlet ground state[19].
The results obtained with the r2SCAN underestimate the experimental ZPL
triplet energy by only 0.15 eV, while LDA has the largest error (0.27 eV). The
energy lowering due to relaxation of atomic coordinates in the singlet states
has been estimated recently using spin-flip TDDFT calculations giving 0.06-0.1
eV for ${}^{1}E$ and 0.02 eV for ${}^{1}A_{1}$[20]. Applying these corrections
to the r2SCAN values for the singlets gives ZPL energy of the singlet
transition close to the experimental estimate, 1.19 eV, but underestimates the
difference between the ${}^{1}A_{1}$ and ${}^{3}E$ excited states. The results
using the TPSS functional provide a more accurate value of the
${}^{1}A_{1}-{}^{3}E$ energy difference. aRef.[84, 85, 19].

Figure 4 compares the ZPL energy for the transition between triplet states,
${}^{3}A_{2}\leftrightarrow{}^{3}E$, calculated using the various density
functionals (values reported in Table 2 in the Appendix) with the experimental
ZPL energy, which is 1.945 eV [84]. Again, r2SCAN gives the best results with
a deviation of -0.15 eV with respect to the experimental value. This is only
slightly larger than the error of previous periodic calculations with the
screened hybrid functional Heyd-Scuseria-Ernzerhof (HSE)[86, 87], where it was
found that the ZPL energy is overestimated by 0.1 eV compared to experiment.
Calculations using the TPSS functional underestimate the ZPL energy of the
triplet transition by 0.2 eV, while calculations using the PBE and LDA
functionals give a deviation of 0.27 eV, similar to previous periodic
calculations using local and semilocal functionals[86]. All functionals
reproduce well the experimentally determined energy lowering after excitation
to the ${}^{3}E$ state, 0.235 eV[84], to within 0.035 eV.

Figure 4 also shows the experimentally deduced excitation energy that includes
the lowering of the energy after excitation to the singlet states. The ZPL
energy for the transition between singlets,
${}^{1}E\leftrightarrow{}^{1}A_{1}$, has been determined to be 1.19 eV[85].
The energy of the singlet states with respect to the triplet ground state can
be deduced from the measured ionization energy of the singlet ground state of
1.91-2.25 eV obtained from recent photoluminescence measurements[19].
Theoretical estimates of the lowering of the energy in the singlet states have
recently been presented by Jin et al.[20] who report 0.06-0.1 eV and 0.02 eV
for the ${}^{1}E$ and ${}^{1}A_{1}$ states, respectively, on the basis of
spin-flip TDDFT calculations with the PBE functional. This suggests that the
ZPL singlet energy should differ from the ${}^{1}E\leftrightarrow{}^{1}A_{1}$
vertical excitation energy by at most 0.08 eV. From this, it can be deduced
that r2SCAN predicts the ZPL energy for the singlet transition with an
accuracy of 0.03-0.07 eV with respect to the experimentally deduced value,
while the other functionals give too small ZPL energy because the energy of
the ${}^{1}A_{1}$ state is underestimated. The r2SCAN functional however
underestimates the energy difference of 0.24 eV between the ${}^{1}A_{1}$ and
${}^{3}E$ excited states. There, the TPSS functional appears to give a better
estimate of $\sim$0.2 eV, in good agreement with the experimentally deduced
value.

##  4 Discussion and conclusions

The variational calculations presented here using several density functionals,
ranging from LDA to meta-GGAs, show that the ordering of the four lowest-
energy states of the NV- center in diamond is actually predicted correctly at
this level of theory, contrary to some earlier reports. The energy of vertical
excitations obtained with the meta-GGA functionals are, furthermore, in close
agreement with the results of accurate but much more computationally intensive
many-body beyond-RPA quantum embedding calculations[21]. The best results are
obtained with the r2SCAN functional, with deviations from the beyond-RPA
values of only $\sim$3% and $\sim$1% for the optical transitions between the
triplet and singlet states, ${}^{3}A_{2}\leftrightarrow{}^{3}E$ and
${}^{1}E\leftrightarrow{}^{1}A_{1}$, respectively.

The calculations presented here are based on a recently developed direct
orbital optimization approach and a plane-wave basis set for a large supercell
with up to 511 atoms subject to periodic boundary conditions in order to
ensure convergence with respect to size and without introducing artifacts due
to finite size and truncated surfaces. The single-determinant wave functions
for the singlet states allow for symmetry breaking, as this is known to be an
important feature for obtaining accurate estimates of the energy of
multiconfigurational ("strongly correlated") systems within a mean-field
approximation[65, 63, 66]. We note that the commonly used self-consistent
field optimization of orbitals typically converges on a symmetric solution if
it is started from a symmetric initial wave function.

Since the calculations are variational, the atomic forces in the excited
states can in principle be evaluated analytically. The minimization of the
energy in the ${}^{3}E$ triplet excited state gives an estimate of the ZPL
energy for the triplet transition in good agreement with experimental
estimates, with r2SCAN underestimating the experimental value by $\sim$0.15
eV. This error is similar to the one previously found for the much more
elaborate and computationally intensive HSE hybrid functional [86, 87].
Optimization of atomic coordinates and simulation of the atomic dynamics in
the singlet states ${}^{1}E$ and ${}^{1}A_{1}$ is, however, not
straightforward because it involves determination of the atomic forces for
multiple single-determinant solutions as indicated in Eqs.5 and 6. Moreover,
since the description of the ${}^{1}A_{1}$ state involves a single determinant
that breaks the spatial symmetry of the wave function, optimization of atomic
coordinates in this state may lead to artificial symmetry breaking of the
atomic configuration of the NV- center. This issue can be alleviated by
introducing a basis of complex orbitals expressed as a linear combination of
the real $e_{x}$ and $e_{y}$ orbitals, where the ${}^{1}A_{1}$ state can be
described with a single-determinant solution that does not break the spatial
symmetry[88].

The values of the excitation energy obtained with the LDA and PBE functionals
are found to be underestimates with respect to the beyond-RPA results,
especially for the excited singlet state, ${}^{1}A_{1}$. This underestimation
can be a consequence of the self-interaction error inherent in calculations
with local and semilocal Kohn-Sham functionals, which can be different for the
different states and thereby affect the excitation energy. An explicit self-
interaction correction can be applied using the orbital based approach
proposed by Perdew and Zunger[89]. Such a correction applied to the PBE
functional has, for example, been found to significantly improve the
calculated energy of excitations of the ethylene molecule[63]. This approach
might also give improved results for the NV- center in diamond, but the
calculations are more involved and computationally demanding as they require
an additional optimization of the orbitals due to the lack of unitary
invariance of the corrected functional [45, 90, 91].

Variational density functional calculations where excited states are obtained
as stationary single-determinant solutions have in some articles in the
literature been described as inadequate for describing electronic excitations
of quantum defects. The results of the calculations presented here show
instead that such an approach can provide accurate energetics for the
electronic states involved in the optical spin initialization in the
prototypical NV- center in diamond. This methodology is, therefore, expected
to be a useful tool for characterizing electronic excitations of other point
defects in materials of interest for quantum applications.

## Acknowledgment

We thank Elvar Ö. Jónsson and Sergei Egorov for helpful discussions. This work
was funded by the Icelandic Research Fund (grant nos. 174582, 217734) and by
the Research Fund of the University of Iceland.

##  Appendix A Relation between multideterminant states and single
determinants

Using Eqs. 3 and 4, the single Slater determinants ${}^{m}\Phi_{2}$ and
${}^{1}\Phi_{3}$ can be expanded as

$\displaystyle{}^{m}\Phi_{2}=\ket{e_{-}\overline{e_{+}}}=(\ket{e_{x}\overline{e_{x}}}-\ket{e_{y}\overline{e_{y}}})/2+(\ket{e_{x}\overline{e_{y}}}-\ket{e_{y}\overline{e_{x}}})/2=$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\Psi\left({}^{1}E\right)+\Psi\left({}^{3}A_{2}\right)\right).$
(7)
$\displaystyle{}^{1}\Phi_{3}=\ket{e_{+}\overline{e_{+}}}=(\ket{e_{x}\overline{e_{x}}}+\ket{e_{y}\overline{e_{y}}})/2+(\ket{e_{x}\overline{e_{y}}}+\ket{e_{y}\overline{e_{x}}})/2=$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\Psi\left({}^{1}A_{1}\right)+\Psi\left({}^{1}E\right)\right).$
(8)

where it has been used that the many-body singlet wave functions are given by
the expressions in Eqs. 1 and 2 and that
$(\ket{e_{x}\overline{e_{y}}}-\ket{e_{y}\overline{e_{x}}})/2$ is the $m_{s}=0$
triplet ground state wave function. Eq. A shows that the single determinant
${}^{m}\Phi_{2}$ is a mixture of the multideterminant ground singlet and
triplet wave functions, $\Psi\left({}^{1}E\right)$ and
$\Psi\left({}^{3}A_{2}\right)$, which leads to spin-symmetry breaking. Eq. A
shows that the single determinant ${}^{1}\Phi_{3}$ is a mixture of the
multideterminant ground and excited singlet wave functions,
$\Psi\left({}^{1}A_{1}\right)$ and $\Psi\left({}^{1}E\right)$, which leads to
spatial-symmetry breaking.

Since the Hamiltonian matrix elements between wave functions of different spin
or spatial symmetry are zero, the energy of the single determinants
${}^{m}\Phi_{2}$ and ${}^{1}\Phi_{3}$ can be expressed as

$\displaystyle\mathcal{E}\left[{}^{m}\Phi_{2}\right]=\left(\mathcal{E}\left[{}^{1}E\right]+\mathcal{E}\left[{}^{3}A_{2}\right]\right)/2,$
(9)
$\displaystyle\mathcal{E}\left[{}^{1}\Phi_{3}\right]=\left(\mathcal{E}\left[{}^{1}E\right]+\mathcal{E}\left[{}^{1}A_{1}\right]\right)/2.$
(10)

Therefore, in the absence of orbital relaxation, the energy of the
multideterminant ground and excited singlet states is given in terms of the
energy of the single determinants ${}^{m}\Phi_{2}$, ${}^{3}\Phi_{1}$ and
${}^{1}\Phi_{3}$ by

$\displaystyle\mathcal{E}\left[{}^{1}E\right]$
$\displaystyle=2\mathcal{E}\left[{}^{m}\Phi_{2}\right]-\mathcal{E}\left[{}^{3}A_{2}\right]=2\mathcal{E}\left[{}^{m}\Phi_{2}\right]-\mathcal{E}\left[{}^{3}\Phi_{1}\right],$
(11) $\displaystyle\mathcal{E}\left[{}^{1}A_{1}\right]$
$\displaystyle=2\mathcal{E}\left[{}^{1}\Phi_{3}\right]-\mathcal{E}\left[{}^{1}E\right]=\mathcal{E}[{}^{3}\Phi_{1}]+2(\mathcal{E}[{}^{1}\Phi_{3}]-\mathcal{E}[{}^{m}\Phi_{2}]),$
(12)

where the energy of the ground triplet state is evaluated from the $m_{s}=1$
single determinant ${}^{3}\Phi_{1}$ instead of the $m_{s}=0$ state. This is a
good approximation since the splitting between the $m_{s}=0$ and $m_{s}=\pm 1$
levels in the triplet ground state is $\sim$2.88 GHz ($\sim$10-5 eV) as
determined by electron paramagnetic resonance measurements[10, 92].

##  Appendix B Vertical excitation and zero-phonon line

|  ${}^{3}A_{2}\leftrightarrow{}^{3}E$ (ZPL) | ${}^{3}A_{2}\leftrightarrow{}^{1}A_{1}$ | ${}^{3}A_{2}\leftrightarrow{}^{1}E$ | ${}^{1}E\leftrightarrow{}^{1}A_{1}$ | ${}^{1}A_{1}\leftrightarrow{}^{3}E$
---|---|---|---|---|---
LDA | 1.889 (1.675) | 1.080 | 0.387 | 0.693 | 0.809
PBE | 1.917 (1.691) | 1.300 | 0.457 | 0.843 | 0.617
TPSS | 1.981 (1.732) | 1.500 | 0.532 | 0.968 | 0.481
SCAN | 2.069 (1.832) | 1.980 | 0.654 | 1.326 | 0.089
r2SCAN | 2.057 (1.789) | 1.804 | 0.621 | 1.183 | 0.253
beyond-RPAa | 2.001 ( - ) | 1.759 | 0.561 | 1.198 | 0.243
Table 2: Vertical excitation energy (in eV) of the NV- center in diamond
obtained from variational calculations using various local and semilocal
density functionals. For the singlet states, the energy has been calculated
according to Eqs. 5 and 6, taking into account the multideterminant character
of these states. The values in parentheses correspond to the calculated zero-
phonon line (ZPL) energy of the optical transition between the triplet states.
aRef.[21].

##  Appendix C Smaller cell calculations

| ${}^{3}A_{2}\leftrightarrow{}^{3}E$ | ${}^{3}A_{2}\leftrightarrow{}^{1}A_{1}$ | ${}^{3}A_{2}\leftrightarrow{}^{1}E$
---|---|---|---
PBE | 1.921 | 1.295 | 0.459

Table 3:  Vertical excitation energy (in eV) obtained from variational
calculations with the PBE functional on a 215-atom supercell of the NV- center
in diamond. The atomic structure of this supercell corresponds to the
structure optimized with PBE in the ground triplet state in Ref. [21]. A
kinetic energy cutoff of 600 eV is used. The values of excitation energy
differ by at most 5 meV from the values obtained using PBE with the larger
511-atom supercell shown in Table 2

.

##  Appendix D Supplemental information

Data related to the results presented here is available at Zenodo [93].

## References

  * [1] R. Schirhagl, K. Chang, M. Loretz and C. L. Degen,  _Nitrogen-Vacancy Centers in Diamond: Nanoscale Sensors for Physics and Biology_ ,  Annual Review of Physical Chemistry 65(1), 83 (2014),  10.1146/annurev-physchem-040513-103659.
  * [2] J. R. Maze, P. L. Stanwix, J. S. Hodges, S. Hong, J. M. Taylor, P. Cappellaro, L. Jiang, M. V. G. Dutt, E. Togan, A. S. Zibrov, A. Yacoby, R. L. Walsworth _et al._ ,  _Nanoscale magnetic sensing with an individual electronic spin in diamond_ ,  Nature 455(7213), 644 (2008),  10.1038/nature07279.
  * [3] G. Balasubramanian, I. Y. Chan, R. Kolesov, M. Al-Hmoud, J. Tisler, C. Shin, C. Kim, A. Wojcik, P. R. Hemmer, A. Krueger, T. Hanke, A. Leitenstorfer _et al._ ,  _Nanoscale imaging magnetometry with diamond spins under ambient conditions_ ,  Nature 455(7213), 648 (2008),  10.1038/nature07278.
  * [4] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth and M. D. Lukin,  _High-sensitivity diamond magnetometer with nanoscale resolution_ ,  Nature Physics 4(10), 810 (2008),  10.1038/nphys1075.
  * [5] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri _et al._ ,  _Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres_ ,  Nature 526(7575), 682 (2015),  10.1038/nature15759.
  * [6] I. Aharonovich, S. Castelletto, D. A. Simpson, C.-H. Su, A. D. Greentree and S. Prawer,  _Diamond-based single-photon emitters_ ,  Reports on Progress in Physics 74(7), 076501 (2011),  10.1088/0034-4885/74/7/076501.
  * [7] A. Beveratos, R. Brouri, T. Gacoin, A. Villing, J.-P. Poizat and P. Grangier,  _Single Photon Quantum Cryptography_ ,  Physical Review Letters 89(18), 187901 (2002),  10.1103/PhysRevLett.89.187901.
  * [8] G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann and J. Wrachtrup,  _Quantum error correction in a solid-state hybrid spin register_ ,  Nature 506(7487), 204 (2014),  10.1038/nature12919.
  * [9] T. H. Taminiau, J. Cramer, T. van der Sar, V. V. Dobrovitski and R. Hanson,  _Universal control and error correction in multi-qubit spin registers in diamond_ ,  Nature Nanotechnology 9(3), 171 (2014),  10.1038/nnano.2014.2.
  * [10] M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup and L. C. Hollenberg,  _The nitrogen-vacancy colour centre in diamond_ ,  Physics Reports 528(1), 1 (2013),  10.1016/j.physrep.2013.02.001.
  * [11] J. R. Weber, W. F. Koehl, J. B. Varley, A. Janotti, B. B. Buckley, C. G. Van de Walle and D. D. Awschalom,  _Quantum computing with defects_ ,  Proceedings of the National Academy of Sciences 107(19), 8513 (2010),  10.1073/pnas.1003052107.
  * [12] P. Neumann, N. Mizuochi, F. Rempp, P. Hemmer, H. Watanabe, S. Yamasaki, V. Jacques, T. Gaebel, F. Jelezko and J. Wrachtrup,  _Multipartite Entanglement Among Single Spins in Diamond_ ,  Science 320(5881), 1326 (2008),  10.1126/science.1157233.
  * [13] P. Delaney, J. C. Greer and J. A. Larsson,  _Spin-Polarization Mechanisms of the Nitrogen-Vacancy Center in Diamond_ ,  Nano Letters 10(2), 610 (2010),  10.1021/nl903646p.
  * [14] A. Gali, M. Fyta and E. Kaxiras,  _Ab initio supercell calculations on nitrogen-vacancy center in diamond: Electronic structure and hyperfine tensors_ ,  Physical Review B 77(15), 155206 (2008),  10.1103/PhysRevB.77.155206.
  * [15] J. P. Goss, R. Jones, S. J. Breuer, P. R. Briddon and S. Öberg,  _The Twelve-Line 1.682 eV Luminescence Center in Diamond and the Vacancy-Silicon Complex_ ,  Physical Review Letters 77(14), 3041 (1996),  10.1103/PhysRevLett.77.3041.
  * [16] S. Choi, M. Jain and S. G. Louie,  _Mechanism for optical initialization of spin in NV - center in diamond_ ,  Physical Review B 86(4), 041202 (2012),  10.1103/PhysRevB.86.041202.
  * [17] V. M. Acosta, A. Jarmola, E. Bauch and D. Budker,  _Optical properties of the nitrogen-vacancy singlet levels in diamond_ ,  Physical Review B - Condensed Matter and Materials Physics 82(20), 2 (2010),  10.1103/PhysRevB.82.201202.
  * [18] L. J. Rogers, S. Armstrong, M. J. Sellars and N. B. Manson,  _Infrared emission of the NV centre in diamond: Zeeman and uniaxial stress studies_ ,  New Journal of Physics 10 (2008),  10.1088/1367-2630/10/10/103024.
  * [19] S. Wolf, I. Meirzada, G. Haim and N. Bar-Gill,  _Nitrogen-vacancy singlet manifold ionization energy_ ,  Physical Review Applied 19, 034076 (2023),  10.1103/PhysRevApplied.19.034076.
  * [20] Y. Jin, M. Govoni and G. Galli,  _Vibrationally resolved optical excitations of the nitrogen-vacancy center in diamond_ ,  npj Computational Materials 8(1), 1 (2022),  10.1038/s41524-022-00928-y.
  * [21] H. Ma, M. Govoni and G. Galli,  _Quantum simulations of materials on near-term quantum computers_ ,  npj Computational Materials 6(1), 85 (2020),  10.1038/s41524-020-00353-z.
  * [22] H. Ma, N. Sheng, M. Govoni and G. Galli,  _Quantum Embedding Theory for Strongly Correlated States in Materials_ ,  Journal of Chemical Theory and Computation 17(4), 2116 (2021),  10.1021/acs.jctc.0c01258.
  * [23] C. Bhandari, A. L. Wysocki, S. E. Economou, P. Dev and K. Park,  _Multiconfigurational study of the negatively charged nitrogen-vacancy center in diamond_ ,  Physical Review B 103(1), 014115 (2021),  10.1103/PhysRevB.103.014115.
  * [24] M. Bockstedte, F. Schütz, T. Garratt, V. Ivády and A. Gali,  _Ab initio description of highly correlated states in defects for realizing quantum bits_ ,  npj Quantum Materials 3(1), 31 (2018),  10.1038/s41535-018-0103-6.
  * [25] A. Ranjbar, M. Babamoradi, M. Heidari Saani, M. A. Vesaghi, K. Esfarjani and Y. Kawazoe,  _Many-electron states of nitrogen-vacancy centers in diamond and spin density calculations_ ,  Physical Review B 84(16), 165212 (2011),  10.1103/PhysRevB.84.165212.
  * [26] Y. Ma, M. Rohlfing and A. Gali,  _Excited states of the negatively charged nitrogen-vacancy color center in diamond_ ,  Physical Review B 81(4), 041204 (2010),  10.1103/PhysRevB.81.041204.
  * [27] A. S. Zyubin, A. M. Mebel, M. Hayashi, H. C. Chang and S. H. Lin,  _Quantum chemical modeling of photoadsorption properties of the nitrogen-vacancy point defect in diamond_ ,  Journal of Computational Chemistry 30(1), 119 (2009),  10.1002/jcc.21042.
  * [28] C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van de Walle,  _First-principles calculations for point defects in solids_ ,  Reviews of Modern Physics 86(1), 253 (2014),  10.1103/RevModPhys.86.253.
  * [29] L. Hedin,  _New method for calculating the one-particle green’s function with application to the electron-gas problem_ ,  Phys. Rev. 139, A796 (1965),  10.1103/PhysRev.139.A796.
  * [30] S. Albrecht, L. Reining, R. Del Sole and G. Onida,  _Ab initio calculation of excitonic effects in the optical spectra of semiconductors_ ,  Phys. Rev. Lett. 80, 4510 (1998),  10.1103/PhysRevLett.80.4510.
  * [31] M. Rohlfing and S. G. Louie,  _Electron-hole excitations in semiconductors and insulators_ ,  Phys. Rev. Lett. 81, 2312 (1998),  10.1103/PhysRevLett.81.2312.
  * [32] A. Gali,  _Ab initio theory of the nitrogen-vacancy center in diamond_ ,  Nanophotonics 8(11), 1907 (2019),  10.1515/nanoph-2019-0154.
  * [33] E. Runge and E. K. U. Gross,  _Density-Functional Theory for Time-Dependent Systems_ ,  Physical Review Letters 52(12), 997 (1984),  10.1103/PhysRevLett.52.997.
  * [34] P. Hohenberg and W. Kohn,  _Inhomogeneous Electron Gas_ ,  Physical Review 136(3B), B864 (1964),  10.1103/PhysRev.136.B864.
  * [35] A. Gali,  _Time-dependent density functional study on the excitation spectrum of point defects in semiconductors_ ,  Physica Status Solidi (b) 248(6), 1337 (2011),  10.1002/pssb.201046254.
  * [36] C.-K. Lin, Y.-H. Wang, H.-C. Chang, M. Hayashi and S. H. Lin,  _One- and two-photon absorption properties of diamond nitrogen-vacancy defect centers: A theoretical study_ ,  The Journal of Chemical Physics 129(12), 124714 (2008),  10.1063/1.2987717.
  * [37] Y. Chen, T. Jiang, H. Chen, E. Han, A. Alavi, K. Yu, E.-G. Wang and J. Chen,  _Multi-configurational nature of electron correlation within nitrogen vacancy centers in diamond_ ,  10.48550/ARXIV.2302.13730 (2023).
  * [38] L. Muechler, D. I. Badrtdinov, A. Hampel, J. Cano, M. Rösner and C. E. Dreyer,  _Quantum embedding methods for correlated excited states of point defects: Case studies and challenges_ ,  Physical Review B 105(23), 235104 (2022),  10.1103/PhysRevB.105.235104.
  * [39] P. W. Ayers, M. Levy and Á. Nagy,  _Time-independent density functional theory for degenerate excited states of Coulomb systems_ ,  Theor. Chem. Acc. 137(11), 152 (2018),  10.1007/s00214-018-2352-7.
  * [40] P. W. Ayers, M. Levy and Á. Nagy,  _Communication: Kohn-Sham theory for excited states of Coulomb systems_ ,  J. Chem. Phys. 143(19), 191101 (2015),  10.1063/1.4934963.
  * [41] P. W. Ayers, M. Levy and A. Nagy,  _Time-independent density-functional theory for excited states of Coulomb systems_ ,  Phys. Rev. A 85(4), 042518 (2012),  10.1103/PhysRevA.85.042518.
  * [42] M. Levy and Á. Nagy,  _Variational Density-Functional Theory for an Individual Excited State_ ,  Phys. Rev. Lett. 83(21), 4361 (1999),  10.1103/PhysRevLett.83.4361.
  * [43] A. Görling,  _Density-functional theory beyond the Hohenberg-Kohn theorem_ ,  Phys. Rev. A 59(5), 3359 (1999),  10.1103/PhysRevA.59.3359.
  * [44] Y. L. A. Schmerwitz, G. Levi and H. Jónsson,  _Calculations of excited electronic states by converging on saddle points using generalized mode following_ ,  10.48550/ARXIV.2302.05912 (2023).
  * [45] A. V. Ivanov, G. Levi, E. Ö. Jónsson and H. Jónsson,  _Method for Calculating Excited Electronic States Using Density Functionals and Direct Orbital Optimization with Real Space Grid or Plane-Wave Basis Set_ ,  Journal of Chemical Theory and Computation 17(8), 5034 (2021),  10.1021/acs.jctc.1c00157.
  * [46] L. E. Daga and L. Maschio,  _Electronic Excitations in Crystalline Solids through the Maximum Overlap Method_ ,  Journal of Chemical Theory and Computation 17(10), 6073 (2021),  10.1021/ACS.JCTC.1C00427.
  * [47] G. Levi, A. V. Ivanov and H. Jónsson,  _Variational Density Functional Calculations of Excited States via Direct Optimization_ ,  Journal of Chemical Theory and Computation 16(11), 6968 (2020),  10.1021/acs.jctc.0c00597.
  * [48] G. Levi, A. V. Ivanov and H. Jónsson,  _Variational calculations of excited states via direct optimization of the orbitals in DFT_ ,  Faraday Discussions 224, 448 (2020),  10.1039/D0FD00064G.
  * [49] D. Hait and M. Head-Gordon,  _Excited state orbital optimization via minimizing the square of the gradient: General approach and application to singly and doubly excited states via density functional theory_ ,  J. Chem. Theory Comput. 16, 1699 (2020),  10.1021/acs.jctc.9b01127,  1911.04709.
  * [50] K. Carter-Fenk and J. M. Herbert,  _State-Targeted Energy Projection: A Simple and Robust Approach to Orbital Relaxation of Non-Aufbau Self-Consistent Field Solutions_ ,  J. Chem. Theory Comput. 16(8), 5067 (2020),  10.1021/acs.jctc.0c00502.
  * [51] G. M. Barca, A. T. Gilbert and P. M. Gill,  _Simple models for difficult electronic excitations_ ,  J. Chem. Theory Comput. 14, 1501 (2018),  10.1021/acs.jctc.7b00994.
  * [52] A. T. Gilbert, N. A. Besley and P. M. Gill,  _Self-consistent field calculations of excited states using the maximum overlap method (MOM)_ ,  J. Phys. Chem. A 112(50), 13164 (2008),  10.1021/jp801738f.
  * [53] J. Gavnholt, T. Olsen, M. Engelund and J. Schiøtz,  _$\Delta$ self-consistent field method to obtain potential energy surfaces of excited molecules on surfaces_,  Phys. Rev. B 78(7), 75441 (2008),  10.1103/PhysRevB.78.075441.
  * [54] A. Hellman, B. Razaznejad and B. I. Lundqvist,  _Potential-energy surfaces for excited states in extended systems_ ,  J. Chem. Phys. 120(10), 4593 (2004),  10.1063/1.1645787.
  * [55] E. Vandaele, M. Mališ and S. Luber,  _The $\Delta$SCF method for non-adiabatic dynamics of systems in the liquid phase_,  J. Chem. Phys 156, 130901 (2022).
  * [56] D. Hait and M. Head-Gordon,  _Orbital optimized density functional theory for electronic excited states_ ,  J. Phys. Chem. Lett. 12, 4517 (2021),  10.1021/acs.jpclett.1c00744.
  * [57] G. Levi, E. Biasin, A. O. Dohn and H. Jónsson,  _On the interplay of solvent and conformational effects in simulated excited-state dynamics of a copper phenanthroline photosensitizer_ ,  Physical Chemistry Chemical Physics 22(2), 748 (2020),  10.1039/c9cp06086c.
  * [58] G. Levi, M. Pápai, N. E. Henriksen, A. O. Dohn and K. B. Møller,  _Solution structure and ultrafast vibrational relaxation of the ptpop complex revealed by $\delta$ scf-qm/mm direct dynamics simulations_,  J. Phys. Chem. C 122, 7100 (2018),  10.1021/acs.jpcc.8b00301.
  * [59] T. Kowalczyk, S. R. Yost and T. V. Voorhis,  _Assessment of the $\Delta$SCF density functional theory approach for electronic excitations in organic dyes_,  J. Chem. Phys. 134(5) (2011),  10.1063/1.3530801.
  * [60] C.-L. Cheng, Q. Wu and T. Van Voorhis,  _Rydberg energies using excited state density functional theory_ ,  J. Chem. Phys. 129(12), 124112 (2008),  10.1063/1.2977989.
  * [61] W. Kohn and L. J. Sham,  _Self-Consistent Equations Including Exchange and Correlation Effects_ ,  Physical Review 140(4A), A1133 (1965),  10.1103/PhysRev.140.A1133.
  * [62] E. Pradhan, K. Sato and A. V. Akimov,  _Non-adiabatic molecular dynamics with $\delta$scf excited states_,  J. Phys. Condens. Matter 30, 484001 (2018),  10.1088/1361-648X/aae864.
  * [63] Y. L. A. Schmerwitz, A. V. Ivanov, E. O. Jónsson, H. Jónsson and G. Levi,  _Variational Density Functional Calculations of Excited States: Conical Intersection and Avoided Crossing in Ethylene Bond Twisting_ ,  The Journal of Physical Chemistry Letters 13(18), 3990 (2022),  10.1021/acs.jpclett.2c00741.
  * [64] R. J. Maurer and K. Reuter,  _Assessing computationally efficient isomerization dynamics: $\Delta$SCF density-functional theory study of azobenzene molecular switching_,  Journal of Chemical Physics 135(22), 1 (2011),  10.1063/1.3664305,  1202.4671.
  * [65] J. P. Perdew, S. Tanvir, R. Chowdhury, C. Shahi, A. D. Kaplan, D. Song and E. J. Bylaska,  _Symmetry Breaking with the SCAN Density Functional Describes Strong Correlation in the Singlet Carbon Dimer_ ,  J. Phys. Chem. A 127, 389 (2023),  10.1021/acs.jpca.2c07590.
  * [66] J. P. Perdew, A. Ruzsinszky, J. Sun, N. K. Nepal and A. D. Kaplan,  _Interpretations of ground-state symmetry breaking and strong correlation in wavefunction and density functional theories_ ,  Proceedings of the National Academy of Sciences of the United States of America 118(4) (2021),  10.1073/pnas.2017850118.
  * [67] J. P. Perdew and Y. Wang,  _Accurate and simple analytic representation of the electron-gas correlation energy_ ,  Physical Review B 45(23), 13244 (1992),  10.1103/PhysRevB.45.13244.
  * [68] A. D. Becke,  _Density-functional exchange-energy approximation with correct asymptotic behavior_ ,  Phys. Rev. A 38, 3098 (1988),  10.1103/PhysRevA.38.3098.
  * [69] J. P. Perdew,  _Density-functional approximation for the correlation energy of the inhomogeneous electron gas_ ,  Phys. Rev. B 33, 8822 (1986),  10.1103/PhysRevB.33.8822.
  * [70] A. J. Cohen, W. Y. P. Mori-Sánchez, P. Mori-Sánchez and W. Yang,  _Insights into Current Limitations of Density Functional Theory_ ,  Science 321(5890), 792 (2008),  10.1126/science.1158722.
  * [71] D. Cremer,  _Density functional theory: Coverage of dynamic and non-dynamic electron correlation effects_ ,  Molecular Physics 99(23), 1899 (2001),  10.1080/00268970110083564.
  * [72] J. P. Perdew, K. Burke and M. Ernzerhof,  _Generalized Gradient Approximation Made Simple_ ,  Physical Review Letters 77(18), 3865 (1996),  10.1103/PhysRevLett.77.3865.
  * [73] J. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria,  _Climbing the Density Functional Ladder: Nonempirical Meta–Generalized Gradient Approximation Designed for Molecules and Solids_ ,  Physical Review Letters 91(14), 146401 (2003),  10.1103/PhysRevLett.91.146401.
  * [74] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew and J. Sun,  _Accurate and numerically efficient r2scan meta-generalized gradient approximation_ ,  The Journal of Physical Chemistry Letters 11(19), 8208 (2020),  10.1021/acs.jpclett.0c02405,  PMID: 32876454.
  * [75] A. Lenef and S. C. Rand,  _Electronic structure of the N- V center in diamond: Theory_ ,  Physical Review B 53(20), 13441 (1996),  10.1103/PhysRevB.53.13441.
  * [76] T. Ziegler, A. Rauk and E. J. Baerends,  _On the calculation of multiplet energies by the hartree-fock-slater method_ ,  Theoretica Chimica Acta 43(3), 261 (1977),  10.1007/BF00551551.
  * [77] U. von Barth,  _Local-density theory of multiplet structure_ ,  Physical Review A 20(4), 1693 (1979),  10.1103/PhysRevA.20.1693.
  * [78] P. E. Blöchl,  _Projector augmented-wave method_ ,  Physical Review B 50(24), 17953 (1994),  10.1103/PhysRevB.50.17953.
  * [79] J. Enkovaara, C. Rostgaard, J. J. Mortensen, J. Chen, M. Dułak, L. Ferrighi, J. Gavnholt, C. Glinsvad, V. Haikola, H. A. Hansen, H. H. Kristoffersen, M. Kuisma _et al._ ,  _Electronic structure calculations with GPAW: a real-space implementation of the projector augmented-wave method_ ,  Journal of Physics: Condensed Matter 22(25), 253202 (2010),  10.1088/0953-8984/22/25/253202.
  * [80] S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques,  _Recent developments in libxc — A comprehensive library of functionals for density functional theory_ ,  SoftwareX 7, 1 (2018),  10.1016/j.softx.2017.11.002.
  * [81] A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings _et al._ ,  _The atomic simulation environment—a Python library for working with atoms_ ,  Journal of Physics: Condensed Matter 29(27), 273002 (2017),  10.1088/1361-648X/aa680e.
  * [82] W. Humphrey, A. Dalke and K. Schulten,  _VMD: Visual molecular dynamics_ ,  Journal of Molecular Graphics 14(1), 33 (1996),  10.1016/0263-7855(96)00018-5.
  * [83] J. Sun, A. Ruzsinszky and J. P. Perdew,  _Strongly Constrained and Appropriately Normed Semilocal Density Functional_ ,  Physical Review Letters 115(3), 036402 (2015),  10.1103/PhysRevLett.115.036402.
  * [84] G. Davies and M. F. Hamer,  _Optical studies of the 1.945 eV vibronic band in diamond_ ,  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 348(1653), 285 (1976),  10.1098/rspa.1976.0039.
  * [85] L. J. Rogers, S. Armstrong, M. J. Sellars and N. B. Manson,  _Infrared emission of the NV centre in diamond: Zeeman and uniaxial stress studies_ ,  New Journal of Physics 10(10), 103024 (2008),  10.1088/1367-2630/10/10/103024.
  * [86] A. Gali, E. Janzén, P. Deák, G. Kresse and E. Kaxiras,  _Theory of Spin-Conserving Excitation of the NV-Center in Diamond_ ,  Physical Review Letters 103(18), 186404 (2009),  10.1103/PhysRevLett.103.186404.
  * [87] A. Alkauskas, B. B. Buckley, D. D. Awschalom and C. G. Van de Walle,  _First-principles theory of the luminescence lineshape for the triplet transition in diamond NV centres_ ,  New Journal of Physics 16(7), 073026 (2014),  10.1088/1367-2630/16/7/073026.
  * [88] G. Thiering and A. Gali,  _Theory of the optical spin-polarization loop of the nitrogen-vacancy center in diamond_ ,  Physical Review B 98(8), 1 (2018),  10.1103/PhysRevB.98.085207.
  * [89] J. P. Perdew and A. Zunger,  _Self-interaction correction to density-functional approximations for many-electron systems_ ,  Physical Review B 23(10), 5048 (1981),  10.1103/PhysRevB.23.5048.
  * [90] S. Lehtola and H. Jónsson,  _Unitary optimization of localized molecular orbitals_ ,  J. Chem. Theory Comput. 9, 5365 (2013).
  * [91] S. Lehtola, M. Head-Gordon and H. Jónsson,  _Complex orbitals, multiple local minima and symmetry breaking in Perdew-Zunger self-interaction corrected density-functional theory calculations_ ,  J. Chem. Theory Comput. 12, 3195 (2016).
  * [92] J. Loubser and J. van Wyk,  _Optical spin-polarisation in a triplet state in irradiated and annealed type 1b diamonds_ ,  Diamond Research pp. 11–15 (1977).
  * [93] A. V. Ivanov, Y. L. A. Schmerwitz, G. Levi and H. Jónsson,  _Data for the manuscript "Electronic excitations of the charged nitrogen-vacancy center in diamond obtained using time-independent variational density functional calculations"_ ,  10.5281/zenodo.7699825 (2023).