# Comprehensive description of color centers by wave function theory: a
CASSCF-NEVPT2 study of the NV defect in diamond

Zsolt Benedek  HUN-REN Wigner Research Centre for Physics, H-1525, Budapest,
Hungary  MTA–ELTE Lendület ”Momentum” NewQubit Research Group, Pázmány Péter
sétány 1/A, 1117 Budapest, Hungary  Ádám Ganyecz  HUN-REN Wigner Research
Centre for Physics, H-1525, Budapest, Hungary  MTA–ELTE Lendület ”Momentum”
NewQubit Research Group, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary
Anton Pershin  HUN-REN Wigner Research Centre for Physics, H-1525, Budapest,
Hungary  Department of Atomic Physics, Institute of Physics, Budapest
University of Technology and Economics, Műegyetem rakpart 3., H-1111 Budapest,
Hungary  Viktor Ivády<EMAIL_ADDRESS>MTA–ELTE Lendület ”Momentum”
NewQubit Research Group, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary
Department of Physics of Complex Systems, Eötvös Loránd University, Egyetem
tér 1-3, H-1053 Budapest, Hungary  Department of Physics, Chemistry and
Biology, Linköping University, SE-581 83 Linköping, Sweden  Gergely Barcza
<EMAIL_ADDRESS>HUN-REN Wigner Research Centre for Physics, H-1525,
Budapest, Hungary  MTA–ELTE Lendület ”Momentum” NewQubit Research Group,
Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

###### Abstract

Paramagnetic point defects in wide-bandgap semiconductors, characterized by
atomic-like in-gap defect states, constitute a unique challenge for ab initio
modeling. In this theoretical study, we aim to devise a wave-function only
computational protocol, exemplified on the prominent nitrogen-vacancy (NV)
center in diamond, which enables the full characterization of future quantum
bit candidates implemented in color centers. We propose the application of the
second order $n$-electron valence state perturbation theory on top of the
complete active space self-consistent field approximation (CASSCF-NEVPT2) to
provide a balanced ab initio level description of the correlation effects
yielded from the defect orbitals and the embedding nanodiamond. By relaxing
the molecular cluster under the compression of the surrounding bulk material,
we manage to model both the vertical and relaxed experimental electronic
spectra within an average error margin of 0.1 eV. Furthermore, the
experimentally observed Jahn-Teller behavior of ${}^{3}E$ and ${}^{1}E$
states, the measured fine structure of the triplet electronic states, as well
as the expected spin-selectivity are quantitatively reproduced by the
presented methodology. Our findings showcase that using conventional
wavefunction-based quantum chemical approaches on carefully crafted cluster
models can be a competing alternative for discussing the energetics of point
defects in solids.

##  I Introduction

Owing to their unique magneto-optical properties, point-like defects in
crystals that act as individual color centers have risen to fame with the
advance of quantum technologies. In the last decade, such solid-state color
centers have been applied as high-resolution nano-sized sensors Degen  _et
al._ (2017); Aslam  _et al._ (2023) thanks to their sensitivity to external
electromagnetic fields, strain, and temperature. Furthermore, a large variety
of single-photon emitters Chunnilall  _et al._ (2014); Aharonovich  _et al._
(2016); Zhang  _et al._ (2020) has been identified in defected solids by now,
which is an integral component in quantum computation Aspuru-Guzik and Walther
(2012) and quantum secure communication Lo  _et al._ (2014). Moreover,
paramagnetic defects that enable spin-selective decay pathways could be used
to create quantum bits Ladd  _et al._ (2010); Awschalom  _et al._ (2013),
controllable through the optically detected magnetic resonance (ODMR)
technique Suter (2020).

From a theoretical point of view, point defects hosted in wide-bandgap
semiconductors behave like atoms featuring localized defect states in a
screening medium of the bulk electrons. Spin-qubit applications relying on
ODMR processes necessitate a complete understanding of the magneto-optical
properties of the defect centers where strongly correlated singlet many-body
states can play a vital role Ádám Gali (2019). Therefore, the proper models of
these color centers require simultaneous high-level treatment of both static
and dynamic correlation effects corresponding to the localized defect and the
embedding solid respectively Ma  _et al._ (2021); Pfäffle  _et al._ (2021);
Muechler  _et al._ (2022).

In the first place, the numerical exploration of the crystalline structures
with point defect Zhang  _et al._ (2020) is typically performed using density
functional theory (DFT) based methods Jones (2015). This approach enables the
computation of many relevant properties of color centers, such as formation
energies, charge transition levels, spin states, hyperfine tensors, zero
phonon lines, and photoluminescence spectra, albeit with varying accuracy, as
summarized in reviews Ref. Freysoldt  _et al._ , 2014; Dreyer  _et al._ ,
2018; Ádám Gali, 2019; Zhang  _et al._ , 2020. However, the widely applied DFT
is an inherently single-determinant method for ground state calculations, and
it has limitations in describing states of strongly multireference nature.
Cohen  _et al._ (2012); Verma and Truhlar (2020) Thus, despite the tremendous
theoretical progress in recent decades in studying correlated electronic
states with DFT Makkar and Ghosh (2021), the quantitative description of solid
state color centers still poses challenges Doherty  _et al._ (2013a); Ádám
Gali (2023).

To further improve the theoretical description of spin-active defects, there
is a strong need for the development of a universally applicable wave function
theory (WFT) based protocol that can accurately handle multiconfigurational
problems. In fact, the defect community has already begun exploring post-DFT
and post-Hartree-Fock methods for this purpose. Without attempting to provide
an exhaustive list, we mention several common methods, including time-
dependent DFT Raghavachari  _et al._ (2002); Gali (2011); Jin  _et al._ (2022,
2023), variational DFT Ivanov  _et al._ (2023), complete active space self-
consistent field (CASSCF)Zyubin  _et al._ (2009); Bhandari  _et al._ (2021),
multireference configuration interaction Zyubin  _et al._ (2009), Monte Carlo
configuration interaction Delaney  _et al._ (2010), configuration interaction
constrained random phase approximation Bockstedte  _et al._ (2018), density
matrix renormalization group approach Barcza  _et al._ (2021), equation of
motion coupled cluster theory Li  _et al._ (2022), GW based approximations Ma
_et al._ (2010); Choi  _et al._ (2012), and various quantum embedding theories
Ma  _et al._ (2020, 2021); Muechler  _et al._ (2022); Sheng  _et al._ (2022);
Chen  _et al._ (2023); Haldar  _et al._ (2023). These methods have mainly been
benchmarked on the negatively charged nitrogen-vacancy (NV-) center in
diamond, which is the most relevant and extensively studied optically active
spin defect to date Ádám Gali (2019). While the distinct models concluded a
largely consistent overall picture of the NV- vertical electronic structure,
there has been less satisfactory quantitative agreement with the experimental
data in most cases. Furthermore, the proper description of magneto-optical
properties, taking into account geometry relaxation effects, that are
important for a detailed understanding of a defect center, has not been fully
addressed within high-level WFT approaches.

In this work, we take a step in this direction by presenting a novel
computational strategy based entirely on conventional quantum chemical
wavefunction approaches. First, we highlight the critical importance of
properly representing the hosting crystalline solid within minimal models of
diamond clusters passivated by hydrogen atoms. By applying the CASSCF approach
on the defect orbitals, we demonstrate that the geometry relaxation only
affects the immediate surroundings of the point defect. The corresponding
CASSCF electronic structure was found to be improved by the NEVPT2
perturbative energy correction, which incorporates the dynamic correlation
effects of the embedding environment. We demonstrate the potential of this
methodology by a comprehensive modeling of the prototypical color center, NV-
in diamond, aligning closely with the most firmly established experimental
observations in the research field.

The rest of the our paper is organized as follows. In Sec. II., we briefly
present the theoretical background of the applied ab initio methods for
interested readers, while the actual details of our performed computations are
given in Sec. III. In Sec. IV, the construction of the cluster models for the
NV- defect center in cubic diamond crystal is presented. Sec. V discusses our
numerical results comparing them to the experimental predictions. Finally, we
conclude on the practical applicability of the presented WFT calculations in
Sec. VI.

##  II Ab initio methods

Our computational approach integrates the CASSCF and NEVPT2 methods to account
for both static and dynamic correlation effects on the electronic structure.
In the following, we briefly summarize the features of the applied approaches,
highlighting the potential advantages of the CASSCF-NEVPT2 method for studying
color centers with highly correlated electronic states arising from the
atomic-like defect orbitals.

###  II.1 Static correlation (CASSCF)

The complete active space self-consistent field (CASSCF) method Siegbahn  _et
al._ (1980); Roos  _et al._ (1980); Siegbahn  _et al._ (1981) captures the
full range of correlation effects within a specific set (”active space”) of
molecular orbitals (MOs). During the CASSCF calculations, a two-step cycle of
the full configuration interaction (FCI) solution and orbital optimization
(i.e. mixing active and external orbitals to minimize the energy) is repeated
until reaching convergence.

Thus this method is capable of providing a practically exact solution in case
all orbitals with non-negligible correlation effects are included in the
active space. Unfortunately, owing to the exponential scaling of the approach
with respect to the number of active orbitals, only a handful of orbitals (up
to 20) can be set to active without further approximations. The remaining
orbitals are kept frozen on the Hartree–Fock level during the solution of the
configuration interaction problem, and only influence the correlation energy
through some orbital mixing during the orbital optimization procedure.

Note that the CASSCF approach is suitable for studying not only the ground
state, but also the many-body excitations, where the orbitals are optimized
either for an ensemble of target states in a state-averaged manner (SA-
CASSCF), or for one selected state in a state-specific manner (SS-CASSCF). The
latter approach is also suitable for geometry relaxation of excited electronic
states.

Overall, the method provides a highly accurate description of the static
correlation (i.e. the mixing of electronic states of different orbital
occupation patterns) but it fails to assess dynamic correlation effect of the
frozen inactive (doubly occupied) and virtual (unoccupied) orbitals. The
resulting CASSCF wave function is qualitatively correct, and it can be used
for geometry optimization and vibrational analysis. However, neglecting the
dynamic correlation generally results in dramatic errors in energy differences
(gaps) between different electronic states.

In computational chemistry Tóth and Pulay (2020), the proper choice of the
active space is not always straightforward. Nevertheless, point defects in
crystals typically yield a few of so-called defect orbitals, which define a
chemically intuitive CAS. Note that these defect-localized MOs lie inside a
large ($>2$ eV) band gap shaping the low-energy excitations. In the case of
the NV center in diamond Doherty  _et al._ (2013b), four relevant defect
orbitals can be identified, specifically those that originate from the
dangling bonds of the three carbons and a nitrogen adjacent to the vacancy.

###  II.2 Dynamic correlation (NEVPT2)

While pristine diamond is known for its lattice structure formed by single
covalent bonds, it is essential to take into account dynamic correlation to
reproduce experimental data Stoll (1992); Booth  _et al._ (2013); McClain  _et
al._ (2017). Accordingly, handling the correlation effects of the crystalline
hosting the point-defect, which are completely omitted by the CASSCF approach,
is critical.

In the literature, various concepts have been developed to provide an _a
posteriori_ correction to the CASSCF solution. Building on our previous work
on point defects in hexagonal boron-nitride layers Babar  _et al._ (2021);
Benedek  _et al._ (2023), we employ the second-order $n$–electron valence
state perturbation theory (NEVPT2) Angeli  _et al._ (2001, 2002, 2007). This
method, a type of multireference perturbative approach, can be thought of as
an extension of the second-order Møller-Plesset perturbation theory Møller and
Plesset (1934) to multireference systems. Owing to the construction of the
zero-order Dyall Hamiltonian Dyall (1995) and to the choice of perturbers, the
NEVPT2 can provide a size-consistent theory on top of the CASSCF reference,
that is free of intruder states and spin-contamination issues.

To speed up the NEVPT2 calculations by using only a limited selection of
perturbers, various contraction schemes were formulated that carefully balance
accuracy with practicality. Angeli  _et al._ (2001) Recent technical
developments, such as the application of the cumulant approximation Guo  _et
al._ (2021), prescreening Kollmar  _et al._ (2021), resolution of identity and
domain-based local pair natural orbitals Guo  _et al._ (2016), have enabled
the routine treatment of molecular systems with thousands of orbitals.

By now, in case of sizeable active spaces, the computational bottleneck of the
CASSCF-NEVPT2 procedure may often no longer be attributed to the NEVPT2 energy
correction itself, but rather to the preceding CASSCF calculation which scales
exponentially with respect to the number of active orbitals. Guo  _et al._
(2023)

###  II.3 Fine structure and spin properties

The splitting of the triplet electronic states induced by secondary magnetic
effects, which is neglected by the non-relativistic CASSCF solution, can also
be estimated a posteriori by quasi-degenerate perturbation theory (QDPT)
Roemelt  _et al._ (2013). Here, the NEVPT2 correction is used to adjust the
Born-Oppenheimer energy contributions of the diagonal QDPT matrix elements.

In numerical implementations, a mean-field approximation of the Breit-Pauli
operator is suitable to describe the second-order spin-orbit coupling (SOC)
operator in case of lighter elements. Neese (2005) With the spin-sublevel
dependent SOC matrix elements at hand, the QDPT matrix is diagonalized,
yielding spin-coupled states as eigenvectors and the respective energy levels
as eigenvalues. Then, the SOC contribution to the zero-field splitting (ZFS)
parameters is extracted as the energy differences. Additionally, the spin-spin
coupling (SSC) generated by the dipolar spin-spin interaction is added to the
mean-field SSC operator to account for its contribution to ZFS. Note that the
inclusion of SSC does not introduce new states, but rather modifies the energy
levels of the obtained QDPT eigenstates.

Figure 1: Construction of the studied molecular models of the NV center.

##  III Technical details of the calculations

In this study, we applied the quantum chemical program package ORCA Neese
(2022) (version 5.0.3) for all computations. The samples of ORCA input files
can be found in the Supporting Information.

For all MODELs, the geometry of all six relevant electronic states were
optimized at state-specific CASSCF(6e,4o)/cc-pVDZ level, by following the
gradients of the root of interest. At the resulting equilibrium geometries,
numerical vibrational analyses (CASSCF(6e,4o)/cc-pVDZ) and single-point energy
calculations (CASSCF(6e,4o)-NEVPT2/cc-pVDZ) were performed.

Calculations at the CASSCF level were performed assuming $C_{1h}$ point group
as it is applicable for both Jahn-Teller stable and Jahn-Teller distorted
structures. Additionally, $E_{x}$ ad $E_{y}$ states can be conveniently
distinguished according to their irreducible representations ($A^{\prime}$ and
$A^{\prime\prime}$, respectively), which is mandatory for state-specific
geometry optimizations. In case when state-specific features such as geometry
relaxation were investigated, the weight of the root of interest was set to
100%, while other roots were kept at 0% (meaning their wavefunctions were
computed but they did not contribute to orbital optimization). On the other
hand, when properties involving multiple states (e.g. transition matrix
elements) were calculated, the six roots were weighed evenly. Vibrational
analysis performed at the CASSCF level of theory was carried out using
numerical gradients without symmetry considerations as available in ORCA.

The CASSCF-NEVPT2 framework was tested on smaller cluster models using both
the strongly contracted and the fully internally contracted implementation of
the ORCA suit, yielding a spectrum with tolerable discrepancy. The results
discussed in the main text were obtained using the fully internally contracted
approach.

All presented results were obtained in cc-pVDZ basis Dunning (1989). Note
that, according to our preliminary analysis performed using cc-pVDZ, cc-pVTZ
and cc-pVQZ basis sets Dunning (1989), only the spectrum of MODEL-1, see Fig.
1 was found to be slightly susceptible to the choice of basis, while larger
clusters were insensitive up to 0.05 eV in average.

To speed up calculations, we took advantage of the resolution of identity
approximation (RI). The construction of Coulomb- and exchange integrals was
carried out in the RIJCOSX framework with def2/J auxiliary basis set. In the
RI-NEVPT2 calculations, cc-pVDZ/C auxiliary basis set was applied. Note that
RI-NEVPT2 implementation of ORCA does not benefit from symmetry. Nevertheless,
the symmetry-adapted CASSCF orbitals were used as input.

The main text focuses solely on the results obtained with MODEL-3. The
complete numerical analysis was also performed for MODEL-1 and 2, confirming
the rapid convergence of our cluster models with respect to the size of
nanodiamond.

##  IV Optimization of cluster models

The NV defect in diamond comprises a substitutional nitrogen atom paired with
an adjacent vacancy. In its negative charge state, it manifests a triplet
ground state characterized by C3v symmetry. Ádám Gali (2019) In our molecular
investigation, we employ quantum chemical models to simulate NV center
embedded within nanodiamonds terminated with hydrogen at the surface. By
progressively scaling up the cluster size, our objective is to accurately
replicate the essential characteristics of the defected bulk crystal.

Large-scale bulk calculations, e.g., Ref. Takács and Ivády, 2024, indicate
clearly that the defect center perturbs the perfect diamond crystal structure
only in its close vicinity. Thus, to reflect the observed stiffness of the
surrounding solid in our cluster models, we optimized atomic positions only
near the vacancy while enforcing the perfect diamond structure in the outer
shells of the cluster.

In the following, we discuss this scheme, depicted in Fig. 1, in more details.
Initially, a sizeable pristine diamond of cubic crystal structure was formed
with C-C bond distances, corresponding to the experimentally measured value of
1.54 Å. After removing a carbon atom close to the center of this structure,
the remaining carbons were divided into ”shells” according to their position
relative to the vacancy: atoms at $n$ chemical bonds distance were assigned to
the $n$th shell. Note that the first and second shells contain 4 and 12
carbons, respectively. In our most compact defected cluster model, denoted as
MODEL-1, we consider a structure where one C atom of the first shell is
substituted by nitrogen and the surface were capped by hydrogens in the
following manner. Each carbon atom of the adjacent third shell were replaced
by either 1, 2, or 3 hydrogens to saturate all dangling bonds of the second
shell of the model. Motivated by hybrid theories Dapprich  _et al._ (1999) and
by preserving the diamond crystal structure, we opt to adjust the length of
these single C-H bonds to the conventional value of 1.09 Å without altering
the positions of the hydrogens any further. During the geometry optimization
process, the position of the terminating hydrogens and the carbons in the
outermost shell are kept fixed, while allowing relaxation only for the atoms
in the inner shells.

Note that this approach requires at least 2 shells around the vacancy; hence,
MODEL-1 is the minimal conceivable size. Larger clusters, MODEL-2 and MODEL-3,
were constructed analogously considering models of carbons up to the third and
fourth shell, respectively.

Figure 2: Left: shape of the frontier Hartree–Fock molecular orbitals which
form the active space for CASSCF(6e,4o) calculations. Right: Expected dominant
electron configurations of NV- electronic states involved in the spin
polarization loop.

By constraining the position of the outer atoms, our model not only considers
the compression of the bulk material, but also the relatively limited number
of atoms to be relaxed, giving the chance to perform geometry optimization on
a high-level theory for all states of interest. While the accurate modeling of
transition energies requires the comprehensive treatment of correlation
effects, the geometry optimization, which involves bond stretching processes,
can be readily performed by capturing merely the static correlations.
Therefore, we applied the CASSCF approach for the constrained relaxation of
the cluster models. These calculations were feasible in practice since in
MODEL-1, 2, and 3 one only has to optimize the position of 4, 16, and 40
atoms, respectively.

##  V Results and discussion

###  V.1 Defect orbitals

We begin our investigation by examining the Hartree–Fock orbitals of the
defective nanodiamond, computed for the ${}^{3}\\!A_{2}$ ground state. The
four defect orbitals of $sp^{3}$ hybrid character situated within the band gap
of the nanodiamond are illustrated in Fig. 2 (left). These orbitals correspond
to the dangling bonds of the three carbons and the nitrogen adjacent to the
vacancy. Among these, the $a_{1}$ is a bonding orbital shared between the
three carbon and nitrogen atoms, while $a_{1}^{*}$ is its anti-bonding pair.
The degenerate orbitals $e_{x}$ and $e_{y}$, which are responsible for the
spin-density distribution, are localised exclusively on the carbon atoms.

Collectively, these four in-gap orbitals are occupied by six valence
electrons, the distribution of which is expected to characterize the low-
energy spectrum Maze  _et al._ (2011). More specifically, the spin-
polarization loop of NV- progresses through the six lowest-lying singlet and
triplet electronic states, i.e., ${{}^{3}}\\!A_{2}$, ${{}^{1}}\\!A_{1}$, and
the doubly degenerate ${{}^{3}}\\!E$, ${{}^{1}}\\!E$ states. Thiering and Gali
(2018) The dominant electron configurations of these states, determined
through group theoretical considerations Maze  _et al._ (2011), are presented
on the right part of Fig. 2. While single-reference methods can approximate
these many-body states to a certain extent, state mixing in the singlets can
become crucial, a complexity addressed by the CASSCF method in our
computational model.

###  V.2 Vertical electronic spectrum

Figure 3: Composition of electronic states at ${}^{3}A_{2}$ geometry,
expressed in the basis of SA-CASSCF natural orbitals. Vertical excitation
energies (CASSCF(6e,4o)-NEVPT2/cc-pVDZ). Accurate experimental data is
available for the vertical singlet and triplet optical transitions 1.26
Kehayias  _et al._ (2013) and 2.18 Davies and Hamer (1976) eV respectively.
The energy shift of the singlet spectrum relative to ground state can be
estimated Bhandari  _et al._ (2021) using the experimentally predicted
${}^{3}E$–${}^{1}A_{1}$ adiabatic gap which was found in the range of
0.32-0.41 eVGoldman  _et al._ (2015a, b, 2017). Note that the
${}^{3}E$-${}^{1}A_{1}$ transition is not directly observable in experiments,
which can lead to an uncertainty in the respective experimental value.

Furthermore, we used the natural orbitals from the state-averaged
CASSCF(6e,4o)/cc-pVDZ solution obtained at the ground state geometry to study
the composition of the electronic states. As shown in Fig. 3, the triplet
states are clearly single-reference, meaning they can be effectively described
by a single electron configuration. On the other hand, the singlet states show
evidence of state-mixing in line with our expectations. Specifically, the
${}^{1}A_{1}$ CASSCF eigenstate contains not only the dominant configurations
proposed by group theory (76%) but also the configuration with an empty
$a_{1}^{*}$ orbital (24%). The ${}^{1}E$ CASSCF eigenstates admixes the spin-
singlet counterpart of the ${}^{3}E$ states (18%) to the group-theoretic
states (82%). Note that the observed state mixing, especially for
${}^{1}A_{1}$, cannot be treated by conventional DFT based methods.

Now, we turn the focus to the discussion of the excitation energies, see Fig.
3. Analyzing the numerical data, we observed that the application of the
NEVPT2 correction on top of the CASSCF theory had a large impact by reducing
the raw CASSCF gaps by 0.2-0.7 eV. In fact, the resulted CASSCF-NEVPT2 optical
transitions are in great agreement with the experimental data deviating by
only $\sim$0.1 eV. Specifically, we predict 2.28 and 1.19 eV for the vertical
triplet and singlet optical gaps while experiments provided 2.18 and 1.26 eV
energies, respectively. From our vertical calculation, we deduce 0.61 eV for
the ${}^{3}E$–${}^{1}A_{1}$ energy gap which is comparable to the indirectly
obtained adiabatic experimental findings of 0.32-0.41 eV.

###  V.3 Relaxed electronic spectrum

Figure 4: Relaxed geometries (CASSCF(6e,4o)/cc-pVDZ), relaxed excitation
energies (CASSCF(6e,4o)-NEVPT2/cc-pVDZ) with ZPE correction at
CASSCF(6e,4o)/cc-pVDZ level. The complementing experimental data for the
adiabatic singlet, triplet optical and ${}^{3}E$–${}^{1}A_{1}$ transitions is
1.19 Rogers  _et al._ (2008), 1.95 Davies and Hamer (1976) and 0.32-0.41
eVGoldman  _et al._ (2015a, b, 2017) eV respectively. The overall energy
spectrum was constructed as given in Ref. 29.

To determine adiabatic energy differences between the states that correspond
to the experimental zero-phonon lines we optimized the geometry of each state
using the state-specific CASSCF(6e,4o)/cc-pVDZ method. The accurate electronic
energies were computed at CASSCF(6e,4o)-NEVPT2/cc-pVDZ level of theory while
the zero point energy (ZPE) corrections were obtained using CASSCF(6e,4o)/cc-
pVDZ.

The relaxed geometries and the excitation energies are summarized in Fig. 4.
We find that ${}^{1}A_{1}$ conserves the $C_{3v}$ symmetry of the
${}^{3}A_{2}$ ground state while the C-C distances are stretched to 2.70 Å.
However, for the $E$ states JT distortion is observed lowering the symmetry to
$C_{1h}$ and forming A’ and A” states. In the former case, one inner C-C
distance elongates relative to the other two, but the system remains close to
$C_{3v}$ symmetry - only 0.02-0.04 Ådeviations were observed. In the latter
case, elongation occurs, and the two longer C-C distances exceed the short C-C
side by 0.09-0.14 Å. The more apparent symmetry breaking of A”, as well as the
orbital picture (Fig. 2), suggest that the $E_{x}\rightarrow A^{\prime\prime}$
configuration is the more stable JT state; namely, the $e_{x}$ orbital, which
has bonding character is more populated in these structures than the
antibonding $e_{y}$.

Nevertheless, the energy difference between $A^{\prime}$ and
$A^{\prime\prime}$ states, which is in the range of ten of _milli_ eVs, cannot
be reliably predicted by CASSCF-NEVPT2 - partly because this small difference
lies within the error margin of NEVPT2, and partly due to the fact that the
optimization was carried out at CASSCF (rather than CASSCF-NEVPT2) level.
Thus, the order of JT states in Fig. 4 reflects mere chemical intuition at
this point, and an average energy level is provided for the $E$ states. A more
sophisticated discussion of JT energy levels can be found in the next section.

Studying the adiabatic spectrum of the color center depicted in Fig. 4, we
found that the CASSCF-NEVPT2 excitation energies are close to the available
experimental values with a deviation of less than 0.1 eV. In particular, for
the transitions ${}^{3}E$–${}^{3}A_{1}$, ${}^{1}E$–${}^{1}A_{1}$ and
${}^{3}E$–${}^{1}A_{1}$ our model neglecting the JT effects predicts 1.90,
1.10 and 0.30 eV which compares excellently to the adiabatic experimental data
1.95, 1.19 and 0.32-0.41 eV respectively.

###  V.4 Jahn-Teller behavior of $E$ states

Figure 5: Schematic representation of potential energy hypersurfaces in the
vicinity of the minimum-energy crossing points (MECPs) between ${}^{3}E$
states (up) and ${}^{1}E$ states (down). The relative energy levels of
$A^{\prime}$, $A^{\prime\prime}$ and MECP structures were calculated at the
level of geometry optimization, i.e. state-specific CASSCF(6e,4o). Jahn-Teller
parameters ($E_{\rm JT},\delta_{\rm JT}$) were obtained from SS-CASSCF energy
differences, as shown, while the energy of the Jahn-Teller phonon ($h\nu_{e}$)
was derived from the vibrational analysis of the MECP geometry, setting 50-50%
weight for the two degenerate $E$ states.

When studying the Jahn-Teller distorted ${}^{3}E$ and ${}^{1}E$ states, it is
of fundamental interest whether the system is trapped in a single potential
energy valley corresponding to an $A^{\prime\prime}$ structure (static Jahn-
Teller effect), or continuously oscillates among three spatially degenerate
$A^{\prime\prime}$ minima (dynamic Jahn-Teller effect). Namely, in the latter
case, the system appears to be of high ($C_{3v}$) symmetry in experiments,
even if the energetically most favorable geometry is distorted.

Static and dynamic Jahn-Teller systems can be distinguished by comparing the
Jahn-Teller barrier ($\delta_{\rm JT}$) to the zero-point energy level of the
two distortion-driving $e$ vibrational modes ($h\nu_{e}$). Here, $\delta_{\rm
JT}$ refers to the energy difference between the two $C_{1h}$ symmetrical
configurations, i.e. $A^{\prime\prime}$ an $A^{\prime}$ (the latter acts as a
transition state in Jahn-Teller oscillation). The criterion for a static JT
effect is $\delta_{\rm JT}>h\nu_{e}$; in the opposite case, the energy of the
vibration is sufficiently large to form a dynamic JT system. Bersuker (2006)

As visualized in Figure 5, $\delta_{\rm JT}$ can be calculated at as an energy
difference $A^{\prime\prime}$ and $A^{\prime}$ geometries. For this
calculation, owing to the expected extremely small gap between the two
equilibrium geometries, we used the level of geometry optimization. Namely,
the SS-CASSCF potential energy surface is expected to be accurate in the
vicinity of the equilibrium geometry, given that a single electronic state is
studied. (A NEVPT2 correction would be disadvantageous for this calculation as
the geometries were not optimized at NEVPT2 level. This discrepancy is
negligible if energy differences of $\approx$ 1 eV size are studied, but
causes large errors in the order of magnitude of meVs.

The Jahn-Teller barrier $\delta_{\rm JT}$ was determined indirectly, as
follows. In the fist step, we searched for the minimum-energy crossing point
(MECP) between $E$ states by setting the weight of both $E_{x}$ and $E_{y}$ to
50% in the CASSCF optimization run. The resulting MECP (geometry for ${}^{3}E$
and ${}^{1}E$ is depicted in Fig. 5, see black frames) is the most stable
geometry where the $E_{x}$ and $E_{y}$ are degenerate, which is only possible
at $C_{3v}$ symmetry. Importantly, MECP ca be considered as a part of both
potential energy surfaces. As indicated in Fig. 5, the energy difference
between MECP and the bottom of the A” valley corresponds to the Jahn-Teller
stabilization energy ($E_{JT}$), the role of which in the modeling is
discussed below.

To gain accurate energy differences, we re-calculated the energy of MECP in
SS-CASSCF calculations, setting the weight of $E_{x}$ or $E_{y}$ to 100%.
These SS-CASSCF calculations target the same energy level in theory, but the
result depends on the selected electronic state due to the slight alteration
of active-space orbital shapes. Then, these energies can be directly compared
to the SS-CASSCF energies for $A^{\prime\prime}$ and $A^{\prime}$, which were
obtained at the end of the geometry optimization run. $E_{JT}$ can be
calculated as the energy difference between MECP and $A^{\prime\prime}$ (using
state-specific orbitals optimized for $E_{x}$); the difference between MECP
and $A^{\prime}$ gives $E_{JT}-\delta_{JT}$ (using state-specific orbitals
optimized for $E_{y}$); finally, $\delta_{JT}$ arises as the difference
between the latter two amounts.

In the end, we obtained $\delta_{\rm JT}$ = 26 and 41 meV for ${}^{3}E$ and
${}^{1}E$, respectively. To determine the nature of the JT effect, we first
assumed static JT where $A^{\prime\prime}$ vibrates in its own potential
valley. The vibrational analysis on this geometry resulted in $h\nu_{e}$
values of 84-88 meV for ${}^{3}A^{\prime\prime}$ and 75-80 meV for
${}^{1}A^{\prime\prime}$ (the degeneracy of $e$ phonons is broken due to the
decreased symmetry in $A^{\prime\prime}$). This, however, clearly contradicts
the initial assumption - the system has enough vibrational energy, even at 0
K, to cross the barrier. Therefore, both ${}^{3}E$ and ${}^{1}E$ will be
handled as dynamic JT systems in the following.

In dynamic JT, the strength of the vibronic coupling, i.e. the proportion of
the Jahn-Teller stabilization energy ($E_{\rm JT}$) and the Jahn-Teller phonon
energy ($h\nu_{e}$) characterizes the behavior of the system; it determines
how the observed properties, such as the fine structure (_vide infra_) is
altered compared to the Born-Oppenheimer picture.

$E_{\rm JT}$ refers to the energy which is released upon symmetry breaking
($C_{3v}\rightarrow C_{1h}$). As it has been mentioned above, it equals to the
energy difference between MECP and $A^{\prime\prime}$. At SS-CASSCF level, we
found $E_{\rm JT}$ = 67 and 143 meV for ${}^{3}E$ and ${}^{1}E$, respectively.

$h\nu_{e}$ for the dynamic system was obtained from a vibrational analysis at
MECP. (We recall that 50-50% weighing was used in CASSCF, under the condition
of which MECP represents a minimum on PES. The actual JT surface is not
differentiable at MECP.) The calculated zero-point JT vibrational energies are
67 and 71 meV for for ${}^{3}E$ and ${}^{1}E$, respectively, as shown in Fig.
5 (black wavy lines).

Putting these data together, both systems turn out to be strongly coupled, as
$h\nu_{e}$ and $E_{\rm JT}$ are commensurable. In the case of ${}^{3}E$, where
$E_{\rm JT}$ is below $h\nu_{e}$, the $C_{3v}$ symmetrical MECP is expected to
serve as a reasonable molecular model for calculating properties. For
${}^{1}E$, on the other hand, stronger vibronic coupling is expected and
computational results from MECP are less relevant.

###  V.5 Fine structure of triplet states

Figure 6: Visualization of the definition of zero-field splitting parameters $D$, $\lambda_{z}$, $D_{es}$ and $\Delta$. The corresponding energy parameters are found in Tab. 1. Some explanation: ${}^{3}E_{x}$, ${}^{3}E_{y}$: two degenerate ${}^{3}E$ states (excitation to $e_{x}$ or $e_{y}$ orbital). Blue numbers - $m_{S}$ quantum numbers. In the case of degenerate states, any linear combination can be an eigenfunction of the respective Hamiltonian. Table 1: Zero-field splitting parameters of triplet electronic states measured in GHz. State | Parameter | Theory | Experiment
---|---|---|---
${}^{3}A_{2}$ | D | 3.10 | 2.88 Loubser and van Wyk (1978)
${}^{3}E$ | $\lambda_{z}$ | 24.7 | -
$p\lambda_{z}$ | 4.20 | 5.3Batalov  _et al._ (2009)
$D_{es}$ | 2.17 | 1.42 Batalov  _et al._ (2009)
$\Delta$ | 2.92 | -
$q\Delta$ | 1.71 | 1.55 Batalov  _et al._ (2009)

By applying quasi-degenerate perturbation theory on the six obtained
CASSCF(6,4)/cc-pVDZ eigenstates, we computed the fine-structure splitting of
the spin-triplet states. The data discussed below were calculated at the
ground state (${}^{3}A_{2}$) equilibrium geometry, and the MECP geometry
between ${}^{3}E_{x}$ and ${}^{3}E_{y}$ in the case of ${}^{3}E$.

The D tensor of the ${}^{3}A_{2}$ state was obtained as the energy difference
between the ground QDPT state (corresponding to ${}^{3}A_{2}$, $m_{S}=0$) and
the two lowest excited states (${}^{3}A_{2}$, $m_{S}\pm 1$). We note that
spin-spin coupling (SSC) effects dominate in ground-state zero-field splitting
- QDPT with only SOC results in three degenerate sublevels for ${}^{3}A_{2}$.
As shown in the top row of Table 1, a near-exact reproduction of experimental
data (D = 3.10 GHz vs the measured D = 2.88 GHz) was achieved.

Next, we investigated the fine structure of the ${}^{3}E$ electronic states.
In the $C_{3v}$ symmetry of the MECP geometry, two orbitally degenerate states
are split into six spin sublevels by SOC and SSC, the relative energy levels
of which can be characterized by three parameters: $\lambda_{z}$, $D_{e}s$ and
$\Delta$, see Fig. 6 for visual explanation. The modeling of this splitting
is, however, more complicated than simply extracting the raw QDPT energies -
some manual processing was required, as described below.

Firstly, as $E_{x}$ and $E_{y}$ electronic states are not entirely degenerate
at the level of computation due to numerical noises, $\approx 0.002eV$ energy
difference is observed at the $C_{3v}$ symmetrical MECP geometry. Therefore,
the electronic energy of the two ${}^{3}E$ states of interest were manually
adjusted to their average value before the ORCA QDPT treatment.

More importantly, the dynamic Jahn-Teller instability of the ${}^{3}E$ states
attenuates the bare $\lambda_{z}$ and $\Delta$ splitting parameters obtained
at MECP by the so-called Ham reduction factors Ham (1968). In our
calculations, the corresponding multiplicators were found to be $p$=0.17 and
$q$=0.59.

Having the proper values at hand, the damped SOC elements $p\lambda_{z}$ and
$q\Delta$ can be readily compared to the available experimental data as
summarized in Table 1. While raw $\lambda_{z}$ is more than 5 times larger
than the experimental value, the Ham-reduced result ($p\lambda_{z}$ = 4.20
GHz) is in reasonable agreement with the measurements (5.3 GHz). Similarly,
after taking into account the $q$ reduction factor, the experimental $\Delta$
splitting of 1.55 GHz is reproduced by our calculations (1.71 GHz) with high
accuracy. We note that the ab initio reproduction of these parameters is
especially challenging due to the exponentialized $\frac{E_{\rm
JT}}{h\nu_{e}}$ term in the reduction factors, which easily introduces
considerable errors. As for the $D_{es}$ parameter, we obtained 2.17 GHz
slightly overestimating the experimental 1.42 GHz.

In total, within the applied computational framework, experimental zero-field-
splitting parameters were recovered with absolute errors of 0.1-1.1 GHz. The
relative error are below 20%, with the sole exception of $D_{es}$, which is
off by a factor of 1.5.

###  V.6 Spin-orbit coupling between singlet and triplet states

Table 2: Spin-orbit coupling matrix elements [GHz], calculated at the equilibrium geometry of the specified initial state. $\Delta E$ refers to the energy difference between states at initial-state geometry, which indicates the distance from the minimum-energy crossing point. Transition | Initial state | Final state | $\Delta E[eV]$ | Channel | SOC matrix element | Experiment Robledo  _et al._ (2011)
---|---|---|---|---|---|---
${}^{3}E~{}\rightarrow~{}^{1}A_{1}$ | ${}^{3}E_{x}~{}\rightarrow~{}^{3}A^{\prime\prime}$ | ${}^{1}A_{1}$ | 0.18 | $m_{s}=\pm 1$ | 6.3 GHz |
$m_{S}=0$ | 5.1 GHz | high selectivity
${}^{3}E_{y}~{}\rightarrow~{}^{3}A^{\prime}$ | ${}^{1}A_{1}$ | 0.18 | $m_{s}=\pm 1$ | 6.0 GHz | (dominant: $m_{s}=\pm 1$ )
$m_{s}=0$ | 0.0 GHz |
${}^{1}E~{}\rightarrow~{}^{3}A_{2}$ | ${}^{1}E_{x}~{}\rightarrow~{}^{1}A^{\prime\prime}$ | ${}^{3}A_{2}$ | 0.40 | $m_{S}=\pm 1$ | 12.3 GHz |
$m_{S}=0$ | 0.0 GHz | low selectivity
${}^{1}E_{y}~{}\rightarrow~{}^{1}A^{\prime}$ | ${}^{3}A_{2}$ | 0.43 | $m_{S}=\pm 1$ | 14.1 GHz | (favored: $m_{S}=0$)
$m_{S}=0$ | 7.5 GHz |
Figure 7: CASSCF wavefuncions of the ${}^{1}E_{y}$ electronic state, on the
basis of SA-CASSCF natural orbitals. Left: At ground-state geometry, where no
spin-orbit coupling is observed with ${}^{3}A_{2}$ at $m_{S}=0$ channel.
Right: At the Jahn-Teller distorted optimized geometry, where the two
determinants of pure ${}^{1}E_{y}$ are no longer equal in weight. This can be
understood as a state mixing between ${}^{1}E_{y}$ and ${}^{1}A_{1}$, which
induces significant spin-orbit coupling.

The probability of transition between singlet and triplet electronic states is
determined in leading order by the corresponding spin-orbit coupling matrix
element (SOCME). We recall that SOCME is a three-dimensional vector, the
square of which is directly proportional to the rate of transition according
to Fermi’s golden rule for intersystem crossing (ISC). Furthermore, the spin
selectivity required for polarization and readout is also introduced to the
theoretical modeling by SOCME, as ISC to/from the $m_{s}=0$ sublevel of the
triplet depends on its $z$-axis (”parallel”) projection, while the other
channel ($m_{s}=\pm 1$) is related to the $xy$-plane (”perpendicular”)
component. While the ab initio calculation of absolute transition rates by
excited state dynamics is a complicated process and is out of the scope of
this paper, the investigation of SOCMEs presented herein already provides some
information about spin-polarization behavior of the system.

When calculating the SOCMEs, an important aspect is the choice of geometry. As
the transition between electronic states occurs most likely in the vicinity of
the minimum-energy crossing point (MECP) between the singlet and triplet
potential energy surfaces, the MECP geometry would be the most rational
choice. Nevertheless, MECP optimization is not as straightforward as the
geometry optimizations for specific electronic states - for example, dynamic
correlation is mandatory, in order to grasp the distance between the two
crossing PESs correctly. As analytical gradients are usually not implemented
for post-CASSCF methods in quantum chemical codes, MECP structures can only be
obtained at large computational cost, which either derives from numerical
gradient or a thorough potential energy surface scan.

Herein, we present an initial investigation, where we assume that the
equilibrium geometry of the initial, higher-energy state is reasonably close
to MECP. In the ISC processes of the polarization cycle of NV-, where the
energy gap between the electronic states involved in ISC is already small
(${}^{3}E\rightarrow^{1}A_{1}:\approx 0.3eV$;
(${}^{1}E\rightarrow^{3}A_{2}:\approx 0.5eV$), this is a reasonable starting
point. Thus, SOCMEs at initial-state geometries will be discussed in the
following.

Firstly, we examine the upper branch of the polarization cycle, i.e. the
${}^{3}E\rightarrow^{1}A_{1}$ transition (Tab. 2, top row). As ${}^{3}E$ is a
DJT system, the coupling from both ${}^{3}E_{x}$ and ${}^{3}E_{y}$ (each at
the respective Jahn-Teller distorted equilibrium geometry) is considered at
state-average CASSCF(6e,4o) level. In both cases, we found significant (6̃
GHz) coupling to ${}^{1}A_{1}$ at the $m_{s}=\pm 1$ channel, while transition
through $m_{s}=0$ is only conceivable from ${}^{3}A^{\prime\prime}$ formed
from ${}^{3}E_{y}$, and at a slightly lower strength of 5.4 GHz. These numbers
are in agreement with experiments, which indicate that $m_{s}=0$ is the
brighter state of NV-, and the intersystem crossing proceeds selectively
through $m_{s}=\pm 1$. The energy difference between the initial and final
states of ISC (${}^{3}E$ and ${}^{1}A_{1}$) is 0.18 eV at initial-state
geometry, which is commensurable to the error margin of NEVPT2. Therefore,
these results closely represent the behavior at MECP, i. e. the actual
geometry of transition where the energy difference ($\Delta E$) is zero.

Next, we moved on to the lower branch, where the ${}^{1}E\rightarrow^{3}A_{2}$
transition closes the polarization circle (Tab. 2, top row). Conspicuously,
the selectivity of $A^{\prime}$ ad $A^{\prime\prime}$ states is reversed
compared to the case of ${}^{3}E$:
${}^{1}E_{x}\rightarrow^{1}$$A^{\prime\prime}$ is a highly spin-selective
process, clearly favoring $m_{s}=\pm 1$ (SOCME = 12.3 GHz), while
${}^{1}E_{y}$ couples to ${}^{1}A^{\prime\prime}$ at both channels. The
underlying reason is that the ground state has $A^{\prime\prime}$ irreducibe
representation in $C_{1h}$ symmetry (in contrast to ${}^{1}A_{1}$, which turns
to $A^{\prime}$), as a result of which the selection rule are reversed. Even
though the $m_{s}=0$ channel, which is required for closing the polarization
circle, is not favored even from the JT-distorted ${}^{1}E_{y}$ (SOCME = 7.5
GHz vs 14.1 GHz for 0 and $\pm 1$ sublevels, respectively), this does not
necessarily mean that our calculations are inaccurate. The relatively high gap
of 0̃.4 eV between ${}^{1}E_{y}$ and ${}^{3}A_{2}$ at the geometry of the
former strongly suggests that we are far from the more relevant MECP
structure. Thus, MECP optimization, which is a subject of our ongoing
research, seems crucial to characterize the lower branch of the nonradiative
pathway.

It is interesting to analyze why the $m_{s}=0$ channel is only accessible
through one of the $E$ states. We illustrate the electronic structure related
reasons on the example of ${}^{1}E_{y}$, the wavefunction of which is provided
in Fig. 7. (For the sake of clarity and consistency, we use SA-CASSCF
wavefuctions and natural orbitals in this representation.) At $C_{3v}$
symmetrical geometry (left side of the Figure), the determinants of the pure
${}^{1}E_{y}$ state (black color) do not provide any coupling to $3A_{2}$, and
even the determinants from state mixing (blue) interact only through
$m_{s}=\pm 1$. On the other hand, after symmetry breaking by Jahn-Teller
distortion (right side), the weight of closed-shell determinants is no longer
equal in the wave function, which can be grasped an additional state mixing
between ${}^{1}E_{y}$ and ${}^{1}A_{1}$ states (see red dashed frames), both
of which possess $A^{\prime}$ irreducible representation at $C_{1h}$ symmetry.
It is this particular mixing (also termed as pseudo-Jahn-Teller effect in the
literature Bersuker (2021) that makes the transition through $m_{s}=0$
feasible. Importantly, $E_{x}$ has a different representation
($A^{\prime\prime}$) even at reduced symmetry, the coupling between $E_{x}$
and ${}^{1}A_{1}$ remains zero, regardless of the geometry.

Altogether, as our wave-function based methodology takes all SOC-inducing
state mixing effects by construction, the calculated spin selectivity of ISC
processes matches the expectations from experiments. Future work will focus on
the ab initio computation of absolute ISC (and photoluminescence) rates.

##  VI Conclusion

In this theoretical study, we aimed to devise a wave-function only
computational protocol which enables the full characterization of future
quantum bit candidates implemented by point defects in semiconductors. The
proposed methodology was tested on the notorious NV- center in diamond, by
reproducing the available experimental data with tolerable error margins.

Our molecular-cluster based framework accounts for static electron correlation
at CASSCF level, placing the chemically relevant defect orbitals to the active
space. Dynamic correlation of the environment was computed at NEVPT2 level
which provides a perturbative correction. The cluster geometries were
optimized on CASSCF level of theory considering the stiffness of the hosting
crystal.

We expect that the presented modeling approach can be routinely applied with
similar accuracy for color centers in covalent crystals featuring up to 6-10
active defect orbitals.

## Acknowledgments

We acknowledge fruitful discussions with Rohit Babar, Gergő Thiering, Ádám
Gali and Michael Roemelt. Helpful computational tweaks learned from the
community of the ORCA Forum are also greatly appreciated.

This research was supported by the National Research, Development, and
Innovation Office of Hungary within the Quantum Information National
Laboratory of Hungary (Grant No. 2022-2.1.1-NL-2022-00004) and within grants
FK 135496 and FK 145395. V.I. also acknowledges the support from the Knut and
Alice Wallenberg Foundation through WBSQD2 project (Grant No. 2018.0071). A.P.
acknowledges the financial support of János Bolyai Research Fellowship of the
Hungarian Academy of Sciences.

The computations were enabled by resources provided by the National Academic
Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National
Infrastructure for Computing (SNIC) at NSC, partially funded by the Swedish
Research Council through grant agreements no. 2022-06725 and no. 2018-05973.
We also acknowledge KIFÜ for awarding us computational resources at the
Komondor supercomputer in Hungary.

## References

## References

  * Degen  _et al._ (2017) C. L. Degen, F. Reinhard, and P. Cappellaro, Rev. Mod. Phys. 89, 035002 (2017).
  * Aslam  _et al._ (2023) N. Aslam, H. Zhou, E. K. Urbach, M. J. Turner, R. L. Walsworth, M. D. Lukin, and H. Park, Nature Reviews Physics 5, 157 (2023).
  * Chunnilall  _et al._ (2014) C. J. Chunnilall, I. P. Degiovanni, S. Kück, I. Müller, and A. G. Sinclair, Optical Engineering 53, 081910 (2014).
  * Aharonovich  _et al._ (2016) I. Aharonovich, D. Englund, and M. Toth, Nature Photonics 10, 631 (2016).
  * Zhang  _et al._ (2020) G. Zhang, Y. Cheng, J.-P. Chou, and A. Gali, Applied Physics Reviews 7, 031308 (2020), https://pubs.aip.org/aip/apr/article-pdf/doi/10.1063/5.0006075/14576561/031308_1_online.pdf .
  * Aspuru-Guzik and Walther (2012) A. Aspuru-Guzik and P. Walther, Nature Physics 8, 285 (2012).
  * Lo  _et al._ (2014) H.-K. Lo, M. Curty, and K. Tamaki, Nature Photonics 8, 595 (2014).
  * Ladd  _et al._ (2010) T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L. O’Brien, Nature 464, 45 (2010).
  * Awschalom  _et al._ (2013) D. D. Awschalom, L. C. Bassett, A. S. Dzurak, E. L. Hu, and J. R. Petta, Science 339, 1174 (2013).
  * Suter (2020) D. Suter, “Optical detection of magnetic resonance,”  (2020).
  * Ádám Gali (2019) Ádám Gali, Nanophotonics 8, 1907 (2019).
  * Ma  _et al._ (2021) H. Ma, N. Sheng, M. Govoni, and G. Galli, Journal of Chemical Theory and Computation 17, 2116 (2021), https://doi.org/10.1021/acs.jctc.0c01258 .
  * Pfäffle  _et al._ (2021) W. Pfäffle, D. Antonov, J. Wrachtrup, and G. Bester, Phys. Rev. B 104, 104105 (2021).
  * Muechler  _et al._ (2022) L. Muechler, D. I. Badrtdinov, A. Hampel, J. Cano, M. Rösner, and C. E. Dreyer, Phys. Rev. B 105, 235104 (2022).
  * Jones (2015) R. O. Jones, Rev. Mod. Phys. 87, 897 (2015).
  * Freysoldt  _et al._ (2014) C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, Rev. Mod. Phys. 86, 253 (2014).
  * Dreyer  _et al._ (2018) C. E. Dreyer, A. Alkauskas, J. L. Lyons, A. Janotti, and C. G. Van de Walle, Annual Review of Materials Research 48, 1 (2018).
  * Cohen  _et al._ (2012) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Chemical Reviews 112, 289 (2012), https://doi.org/10.1021/cr200107z .
  * Verma and Truhlar (2020) P. Verma and D. G. Truhlar, Trends in Chemistry 2, 302 (2020), special Issue - Laying Groundwork for the Future.
  * Makkar and Ghosh (2021) P. Makkar and N. N. Ghosh, RSC Adv. 11, 27897 (2021).
  * Doherty  _et al._ (2013a) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. Hollenberg, Physics Reports 528, 1 (2013a), the nitrogen-vacancy colour centre in diamond.
  * Ádám Gali (2023) Ádám Gali, Nanophotonics 12, 359 (2023).
  * Raghavachari  _et al._ (2002) K. Raghavachari, D. Ricci, and G. Pacchioni, The Journal of Chemical Physics 116, 825 (2002), https://pubs.aip.org/aip/jcp/article-pdf/116/2/825/19324682/825_1_online.pdf .
  * Gali (2011) A. Gali, physica status solidi (b) 248, 1337 (2011), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.201046254 .
  * Jin  _et al._ (2022) Y. Jin, M. Govoni, and G. Galli, npj Computational Materials 8, 238 (2022).
  * Jin  _et al._ (2023) Y. Jin, V. W.-z. Yu, M. Govoni, A. C. Xu, and G. Galli, Journal of Chemical Theory and Computation 19, 8689 (2023), https://doi.org/10.1021/acs.jctc.3c00986 .
  * Ivanov  _et al._ (2023) A. V. Ivanov, Y. L. A. Schmerwitz, G. Levi, and H. Jónsson, SciPost Phys. 15, 009 (2023).
  * Zyubin  _et al._ (2009) A. S. Zyubin, A. M. Mebel, M. Hayashi, H. C. Chang, and S. H. Lin, Journal of Computational Chemistry 30, 119 (2009), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.21042 .
  * Bhandari  _et al._ (2021) C. Bhandari, A. L. Wysocki, S. E. Economou, P. Dev, and K. Park, Phys. Rev. B 103, 014115 (2021).
  * Delaney  _et al._ (2010) P. Delaney, J. C. Greer, and J. A. Larsson, Nano Letters 10, 610 (2010), https://doi.org/10.1021/nl903646p .
  * Bockstedte  _et al._ (2018) M. Bockstedte, F. Schütz, T. Garratt, V. Ivády, and A. Gali, npj Quantum Materials 3, 31 (2018).
  * Barcza  _et al._ (2021) G. Barcza, V. Ivády, T. Szilvási, M. Vörös, L. Veis, A. Gali, and O. Legeza, Journal of Chemical Theory and Computation 17, 1143 (2021).
  * Li  _et al._ (2022) M. Li, R. Kobayashi, R. D. Amos, M. J. Ford, and J. R. Reimers, Chemical Science 13, 1492 (2022).
  * Ma  _et al._ (2010) Y. Ma, M. Rohlfing, and A. Gali, Phys. Rev. B 81, 041204 (2010).
  * Choi  _et al._ (2012) S. Choi, M. Jain, and S. G. Louie, Phys. Rev. B 86, 041202 (2012).
  * Ma  _et al._ (2020) H. Ma, M. Govoni, and G. Galli, npj Computational Materials 6, 85 (2020).
  * Sheng  _et al._ (2022) N. Sheng, C. Vorwerk, M. Govoni, and G. Galli, Journal of Chemical Theory and Computation 18, 3512 (2022), pMID: 35648660, https://doi.org/10.1021/acs.jctc.2c00240 .
  * Chen  _et al._ (2023) Y. Chen, T. Jiang, H. Chen, E. Han, A. Alavi, K. Yu, E. Wang, and J. Chen, Phys. Rev. B 108, 045111 (2023).
  * Haldar  _et al._ (2023) S. Haldar, A. Mitra, M. R. Hermes, and L. Gagliardi, The Journal of Physical Chemistry Letters 14, 4273 (2023), pMID: 37126760, https://doi.org/10.1021/acs.jpclett.3c00551 .
  * Siegbahn  _et al._ (1980) P. Siegbahn, A. Heiberg, B. Roos, and B. Levy, Physica Scripta 21, 323 (1980).
  * Roos  _et al._ (1980) B. O. Roos, P. R. Taylor, and P. E. Sigbahn, Chemical Physics 48, 157 (1980).
  * Siegbahn  _et al._ (1981) P. E. M. Siegbahn, J. Almlöf, A. Heiberg, and B. O. Roos, The Journal of Chemical Physics 74, 2384 (1981), https://pubs.aip.org/aip/jcp/article-pdf/74/4/2384/18926917/2384_1_online.pdf .
  * Tóth and Pulay (2020) Z. Tóth and P. Pulay, Journal of Chemical Theory and Computation 16, 7328 (2020), https://doi.org/10.1021/acs.jctc.0c00123 .
  * Doherty  _et al._ (2013b) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. Hollenberg, Physics Reports 528, 1 (2013b).
  * Stoll (1992) H. Stoll, Phys. Rev. B 46, 6700 (1992).
  * Booth  _et al._ (2013) G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Nature 493, 365 (2013).
  * McClain  _et al._ (2017) J. McClain, Q. Sun, G. K.-L. Chan, and T. C. Berkelbach, Journal of Chemical Theory and Computation 13, 1209 (2017), https://doi.org/10.1021/acs.jctc.7b00049 .
  * Babar  _et al._ (2021) R. Babar, G. Barcza, A. Pershin, H. Park, O. B. Lindvall, G. Thiering, O. Legeza, J. H. Warner, I. A. Abrikosov, A. Gali, and V. Ivády, arXiv:2111.09589 [cond-mat] (2021), 2111.09589 .
  * Benedek  _et al._ (2023) Z. Benedek, R. Babar, Á. Ganyecz, T. Szilvási, Ö. Legeza, G. Barcza, and V. Ivády, npj Computational Materials 9, 187 (2023).
  * Angeli  _et al._ (2001) C. Angeli, R. Cimiraglia, S. Evangelisti, T. Leininger, and J.-P. Malrieu, The Journal of Chemical Physics 114, 10252 (2001).
  * Angeli  _et al._ (2002) C. Angeli, R. Cimiraglia, and J.-P. Malrieu, The Journal of Chemical Physics 117, 9138 (2002), https://doi.org/10.1063/1.1515317 .
  * Angeli  _et al._ (2007) C. Angeli, M. Pastore, and R. Cimiraglia, Theoretical Chemistry Accounts 117, 743 (2007).
  * Møller and Plesset (1934) C. Møller and M. S. Plesset, Phys. Rev. 46, 618 (1934).
  * Dyall (1995) K. G. Dyall, The Journal of Chemical Physics 102, 4909 (1995), https://pubs.aip.org/aip/jcp/article-pdf/102/12/4909/19185508/4909_1_online.pdf .
  * Guo  _et al._ (2021) Y. Guo, K. Sivalingam, and F. Neese, The Journal of Chemical Physics 154, 214111 (2021), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0051211/16673662/214111_1_online.pdf .
  * Kollmar  _et al._ (2021) C. Kollmar, K. Sivalingam, Y. Guo, and F. Neese, The Journal of Chemical Physics 155, 234104 (2021), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0072129/16704053/234104_1_online.pdf .
  * Guo  _et al._ (2016) Y. Guo, K. Sivalingam, E. F. Valeev, and F. Neese, J. Chem. Phys. 144, 094111 (2016).
  * Guo  _et al._ (2023) Y. Guo, F. Pavosevic, K. Sivalingam, U. Becker, E. Valeev, and F. Neese, The Journal of Chemical Physics 158 (2023), 10.1063/5.0144260.
  * Roemelt  _et al._ (2013) M. Roemelt, D. Maganas, S. DeBeer, and F. Neese, The Journal of Chemical Physics 138, 204101 (2013), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4804607/13280319/204101_1_online.pdf .
  * Neese (2005) F. Neese, The Journal of Chemical Physics 122, 034107 (2005), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.1829047/10872245/034107_1_online.pdf .
  * Neese (2022) F. Neese, Wiley Interdisciplinary Reviews: Computational Molecular Science 12, e1606 (2022).
  * Dunning (1989) T. H. Dunning, The Journal of Chemical Physics 90, 1007 (1989), https://doi.org/10.1063/1.456153 .
  * Takács and Ivády (2024) I. Takács and V. Ivády, “Accurate hyperfine tensors for solid state quantum applications: Case of the nv center in diamond,” (2024), arXiv:2309.03983 [quant-ph] .
  * Dapprich  _et al._ (1999) S. Dapprich, I. Komáromi, K. Byun, K. Morokuma, and M. J. Frisch, Journal of Molecular Structure: THEOCHEM 461-462, 1 (1999).
  * Maze  _et al._ (2011) J. R. Maze, A. Gali, E. Togan, Y. Chu, A. Trifonov, E. Kaxiras, and M. D. Lukin, New Journal of Physics 13, 025025 (2011).
  * Thiering and Gali (2018) G. m. H. Thiering and A. Gali, Phys. Rev. B 98, 085207 (2018).
  * Kehayias  _et al._ (2013) P. Kehayias, M. W. Doherty, D. English, R. Fischer, A. Jarmola, K. Jensen, N. Leefer, P. Hemmer, N. B. Manson, and D. Budker, Phys. Rev. B 88, 165202 (2013).
  * Davies and Hamer (1976) G. Davies and M. F. Hamer, Proceedings of the Royal Society of London. Series A 348, 285 (1976).
  * Goldman  _et al._ (2015a) M. L. Goldman, A. Sipahigil, M. W. Doherty, N. Y. Yao, S. D. Bennett, M. Markham, D. J. Twitchen, N. B. Manson, A. Kubanek, and M. D. Lukin, Phys. Rev. Lett. 114, 145502 (2015a).
  * Goldman  _et al._ (2015b) M. L. Goldman, M. W. Doherty, A. Sipahigil, N. Y. Yao, S. D. Bennett, N. B. Manson, A. Kubanek, and M. D. Lukin, Phys. Rev. B 91, 165201 (2015b).
  * Goldman  _et al._ (2017) M. L. Goldman, M. W. Doherty, A. Sipahigil, N. Y. Yao, S. D. Bennett, N. B. Manson, A. Kubanek, and M. D. Lukin, Phys. Rev. B 96, 039905 (2017).
  * Rogers  _et al._ (2008) L. J. Rogers, S. Armstrong, M. J. Sellars, and N. B. Manson, New Journal of Physics 10, 103024 (2008).
  * Bersuker (2006) I. Bersuker,  _The Jahn-Teller Effect_ (University of Texas, Austin, 2006).
  * Loubser and van Wyk (1978) J. H. N. Loubser and J. A. van Wyk, Reports on Progress in Physics 41, 1201 (1978).
  * Batalov  _et al._ (2009) A. Batalov, V. Jacques, F. Kaiser, P. Siyushev, P. Neumann, L. J. Rogers, R. L. McMurtrie, N. B. Manson, F. Jelezko, and J. Wrachtrup, Phys. Rev. Lett. 102, 195506 (2009).
  * Ham (1968) F. S. Ham, Phys. Rev. 166, 307 (1968).
  * Robledo  _et al._ (2011) L. Robledo, H. Bernien, T. van der Sar, and R. Hanson, New Journal of Physics 13, 025013 (2011).
  * Bersuker (2021) I. B. Bersuker, Chemical Reviews 121, 1463 (2021), https://doi.org/10.1021/acs.chemrev.0c00718 .