# Fluorescence spectrum and charge state control of divacancy qubits via
illumination at elevated temperatures in 4H silicon carbide

A. Csóré  Department of Atomic Physics, Budapest University of Technology and
Economics, Budafoki út 8., H-1111, Budapest, Hungary  I. G. Ivanov  Department
of Physics, Chemistry and Biology, Linköping University, SE-58183 Linköping,
Sweden  N. T. Son  Department of Physics, Chemistry and Biology, Linköping
University, SE-58183 Linköping, Sweden  A. Gali  Department of Atomic Physics,
Budapest University of Technology and Economics, Budafoki út 8., H-1111,
Budapest, Hungary  Wigner Research Centre for Physics, PO. Box 49, Budapest
H-1525, Hungary

(August 27, 2024)

###### Abstract

Divacancy in its neutral charge state (V${}_{\text{C}}$V${}_{\text{Si}}^{0}$)
in 4H silicon carbide (SiC) is a leading quantum bit (qubit) contender. Owing
to the lattice structure of 4H SiC, four different
V${}_{\text{C}}$V${}_{\text{Si}}$ configurations can be formed. The ground and
the optically accessible excited states of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ configurations exhibit a high-spin
state, and the corresponding optical transition energies are around $\approx
1.1$ eV falling in the near-infrared wavelength region. Recently,
photoluminescence (PL) quenching has been experimentally observed for all
V${}_{\text{C}}$V${}_{\text{Si}}$ configurations in 4H SiC at cryogenic
temperatures. It has been shown that V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ is
converted to V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ and it remains in this
shelving dark state at cryogenic temperatures until photoexcitation with the
threshold energies or above is applied to convert
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ back to
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$. In this study, we demonstrate both in
experiments and theory that the threshold energy for reionization is
temperature dependent. We carry out density functional theory (DFT)
calculations in order to investigate the temperature dependent reionization
spectrum, i.e., the spectrum of the
V${}_{\text{C}}$V${}_{\text{Si}}^{-}\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ process. We find that simultaneous
optical reionization and qubit manipulation can be carried out at room
temperature with photoexcitation at the typical excitation wavelength used for
readout of the divacancy qubits in 4H SiC, in agreement with our experimental
data. We also provide the analysis of the PL spectrum of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$, characteristic for each
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ configuration in 4H SiC, using the
Huang-Rhys theory, and find that one configuration in 4H SiC stands out in
terms of the strength of coherent emission among the four configurations.

###### pacs:

Valid PACS appear here

##  I Introduction

Spin active point defects in wide band gap semiconductors have proven to be
promising candidates for numerous applications in the rapidly emerging field
of quantum technology. Leading contender of this family is the negatively
charged nitrogen-vacancy defect, i.e., the NV center in diamond with
exhibiting many desirable magneto-optical properties Doherty  _et al._ (2013);
Gali (2019). Nevertheless, further alternative point defect color centers have
been searched for in the past decade embedded in technologically more mature
materials than diamond.

A suitable host candidate to this end is silicon carbide (SiC) introducing
different crystal modifications called polytypes. The technologically most
relevant polytype is the hexagonal 4H SiC which hosts different types of
quantum bits (qubits) Koehl  _et al._ (2011); Castelletto  _et al._ (2013);
Christle  _et al._ (2015); Widmann  _et al._ (2015); Wolfowicz  _et al._
(2020) being pivotal in quantum information processing Lohrmann  _et al._
(2015); Radulaski  _et al._ (2017); Spindlberger  _et al._ (2019); Nagy  _et
al._ (2019); Babin  _et al._ (2022), ultrasensitive nanosensors in
magnetometry Lee  _et al._ (2015); Kraus  _et al._ (2014); Simin  _et al._
(2016, 2015); Niethammer  _et al._ (2016); Cochrane  _et al._ (2016) or
thermometry Kraus  _et al._ (2014); Anisimov  _et al._ (2016).

Divacancy is one of the most successful qubits in 4H SiC Gali (2011); Koehl
_et al._ (2011) which consists of neighbor carbon and silicon vacancies,
V${}_{\text{C}}$V${}_{\text{Si}}$. 4H SiC consists of Si-C bilayers exhibiting
hexagonal ($h$) or quasicubic ($k$) sites. Consequently,
V${}_{\text{C}}$V${}_{\text{Si}}$ possesses four different defect
configurations in 4H SiC because of the two possible defect sites for both
V${}_{\text{C}}$ and V${}_{\text{Si}}$. Henceforward, we denote the different
configurations with the lattice sites’ labels of the C and Si vacancies,
respectively, such as V${}_{\text{C}}$V${}_{\text{Si}}$($hh$),
V${}_{\text{C}}$V${}_{\text{Si}}$($kk$),
V${}_{\text{C}}$V${}_{\text{Si}}$($hk$) and
V${}_{\text{C}}$V${}_{\text{Si}}$($kh$). The $hh$ and $kk$ configurations are
called axial or on-axis configurations exhibiting C${}_{\text{3v}}$ symmetry,
while the $hk$ and $kh$ configurations are basal or off-axis configurations
with C${}_{\text{1h}}$ symmetry. Lattice structure of 4H SiC and the
V${}_{\text{C}}$V${}_{\text{Si}}$ defect configurations are shown in Fig. 1.

Figure 1: (a) Lattice structure of 4H SiC. The [0001] crystal axis ($c$-axis)
and the $h$/$k$ bilayers are depicted. (b) Microscopic structure of the
divacancy defect. (c) The four possible V${}_{\text{C}}$V${}_{\text{Si}}$
configurations embedded into 4H SiC. Color code of the atoms and the vacant
sites is indicated.

The divacancy in 4H SiC exhibits several desirable magneto-optical properties
making it a promising qubit. In particular, electronic structure of the
neutral V${}_{\text{C}}$V${}_{\text{Si}}$ center
(V${}_{\text{C}}$V${}_{\text{Si}}^{0}$) introduces high-spin ($S=1$) ground
and optically accessible excited states Magnusson and Janzén (2005); Baranov
_et al._ (2005); Son  _et al._ (2006); Koehl  _et al._ (2011). The
corresponding no-phonon optical transitions, i.e., the zero-phonon lines
(ZPLs) are around 1100 nm falling into the near-infrared (NIR) region where
optical fibers and organic tissues have relatively low absorption window. This
property makes the divacancy a promising candidate as a single photon emitter
with quantum memory for quantum communication and as an ultrasensitive
biomarker Somogyi and Gali (2014); Somogyi  _et al._ (2012); Beke  _et al._
(2020). Indeed, high-fidelity infrared spin-photon interface Christle  _et
al._ (2017) and entaglement and control of single nuclear spins Bourassa  _et
al._ (2020) have been demonstrated with divacancies in 4H SiC. The ZPLs are
now labeled Koehl  _et al._ (2011) and identified with the corresponding
configurations Falk  _et al._ (2014) as PL1 -
V${}_{\text{C}}$V${}_{\text{Si}}$($hh$), PL2 -
V${}_{\text{C}}$V${}_{\text{Si}}$($kk$), PL3 -
V${}_{\text{C}}$V${}_{\text{Si}}$($hk$) and PL4 -
V${}_{\text{C}}$V${}_{\text{Si}}$($kh$) that were previously labeled as UD2
centers in 4H SiC in the literature Magnusson and Janzén (2005).

The charge state control of defect qubits is indispensable for robust quantum
information processing Wolfowicz  _et al._ (2021). Photoexcitation of target
point defect qubits may lead to their unwanted charge switching with losing
the qubit state. According to previous ab initio studies Gali  _et al._
(2006); Gali (2012); Gordon  _et al._ (2015); Magnusson  _et al._ (2018), all
the divacancy configurations could appear in the positive, neutral, negative
and double negative charge states depending on the position of the Fermi-level
in 4H SiC which may be actively controlled by doping and electrical switching
in a 4H SiC diode Anderson  _et al._ (2019). Nevertheless, photoexcitation of
the neutral divacancies may lead to the ionization of the defect. Indeed,
photoluminescence (PL) quenching of divacancy defects in 4H SiC has been
recently observed, i.e., the PL signals of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defects do not appear upon certain
excitation energies Zwier  _et al._ (2015); Wolfowicz  _et al._ (2017); Golter
and Lai (2017); Magnusson  _et al._ (2018). On the other hand, it has been
demonstrated that applying a second laser (repump laser) with a higher photon
energy completely recovers the PL intensity of the
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ centers Zwier  _et al._ (2015);
Wolfowicz  _et al._ (2017); Golter and Lai (2017); Magnusson  _et al._ (2018).
Here, we note that the PL quenching phenomenon is not unique for
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ in 4H SiC but was observed for several
other defects (e.g., see Ref. Gali, 2019 and references therein for NV center
in diamond). The physical mechanism behind the quenching phenomenon may be the
charge state switching of the corresponding defect upon illumination as
already demonstrated by previous experimental and theoretical works Zwier  _et
al._ (2015); Wolfowicz  _et al._ (2017); Golter and Lai (2017); Siyushev  _et
al._ (2013); Magnusson  _et al._ (2018). In this case, the role of the repump
laser is to reionize the defect, i.e., to maintain the neutral charge state.
In particular, the ’dark’ charge state has been attributed to either the
singly positive Golter and Lai (2017) or the singly negative charge state
Wolfowicz  _et al._ (2017) of divacancies by earlier experimental works. In
order to unravel this issue, we carried out a combined experimental and
density functional theory (DFT) study Magnusson  _et al._ (2018) and found
that the singly negative charge state [V${}_{\text{C}}$V${}_{\text{Si}}^{-}$]
is the dark state because the calculated and observed photoexcitation
threshold energies at cryogenic temperatures agreed well for this charge state
for all the divacancy configurations in 4H SiC.

In this work, we investigate the temperature dependent charge state control of
the V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defect configurations in 4H SiC. To
this end, we calculate the V${}_{\text{C}}$V${}_{\text{Si}}^{0}$-related PL
lineshape and the V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ reionization spectra based on the Huang-
Rhys (HR) theory Huang  _et al._ (1950); Alkauskas  _et al._ (2014); Gali  _et
al._ (2016), and monitor the apperance of the divacancy PL spectra as a
function of the applied repump laser’s wavelength and temperature in high-
purity semi-insulating 4H SiC samples. We find an excitation scheme with
simultaneous charge state control and qubit manipulation of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ by a single color laser beam at elevated
temperatures. Our paper is organized as follows: we describe the computational
methodology in Sec. II, where we briefly describe the HR theory in Sec. II.1,
the temperature dependence of the charge transition levels in Sec. II.2 and
the methodology of the ab initio DFT calculations in Sec. II.3. We report and
discuss our results in Sec. III. In particular, PL spectra at cryogenic
temperature of the V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defect configurations
are discussed in Sec. III.1, the simulated
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ temperature dependent reionization
spectra are reported within various levels of approximation in Secs.
III.2-III.3. The experimental data are presented and interpreted by our DFT
simulations in Sec. III.4. We conclude our work in Sec. IV.

##  II Methodology

###  II.1  Huang-Rhys theory

The PL spectrum [$I^{\text{EM}}(\hbar\omega)$] within HR theory can be
expressed as Alkauskas  _et al._ (2014)

$I^{\text{EM}}(\hbar\omega)=C^{\text{EM}}\omega^{3}\mu_{\text{EG}}^{2}\sum_{n,m}c_{m}(T)\Big{|}\Big{\langle}\Theta_{\text{ES}}^{m}(Q)\Big{|}\Theta_{\text{GS}}^{n}(Q)\Big{\rangle}\Big{|}^{2}\delta(E_{\text{ES}}^{m}-E_{\text{GS}}^{n}-\hbar\omega)\text{,}$
(1)

where the transitions from the corresponding vibronic substates of the excited
state (ES), indexed by $m$, to those of the ground state (GS) with the index
$n$ are also considered establishing the phonon sideband (PSB) besides the
ZPL, i.e., $m=n=0$ in $I^{\text{EM}}(\hbar\omega)$ yields the intensity of
ZPL. Generalized coordinate $Q$ appears in both the ES and GS vibrational
function $\Theta$ according to the HR theory. The temperature dependence
appears by the factor $c_{m}(T)$ with temperature $T$ which will be defined
later. In Eq. 1 $\boldsymbol{\mu}_{\text{EG}}$ is the corresponding matrix
element (which is now a three component vector with the length of
$\mu_{\text{EG}}$) of the transition dipole moment operator,
$\boldsymbol{\hat{\mu}}_{\text{EG}}=q\sum_{j}\mathbf{\hat{r}}_{j}$ defined as

$\boldsymbol{\mu}_{\text{EG}}=\Big{\langle}\Psi_{\text{ES}}(\mathbf{r}_{i})\Big{|}\boldsymbol{\hat{\mu}}_{\text{EG}}\Big{|}\Psi_{\text{GS}}(\mathbf{r}_{i})\Big{\rangle}$
(2)

and $\hbar\omega$ is the photon energy. The prefactor $\omega^{3}$ in Eq. 1
consists of the photon density of states (DOS) causing spontaneous emission
($\sim\omega^{2}$) and the perturbing field of those photons ($\sim\omega$).
In contrast, the absorption spectrum is a linear function of $\omega$ since no
spontaneous emission is involved. Thus the absorption spectrum can be
expressed as

$I^{\text{ABS}}(\hbar\omega)=C^{\text{ABS}}\mu_{\text{EG}}^{2}\omega\mathop{\raisebox{-20.00024pt}{\scalebox{2.0}{$\displaystyle\sum$}}}_{n,m}c_{n}(T)\Big{|}\Big{\langle}\Theta^{n}_{\text{GS}}(Q)\Big{|}\Theta^{m}_{\text{ES}}(Q)\Big{\rangle}\Big{|}^{2}\delta(E^{m}_{\text{ES}}-E^{n}_{\text{GS}}-\hbar\omega)\text{,}$
(3)

where $C^{\text{EM}}$ and $C^{\text{ABS}}$ are the corresponding material and
measurement dependent constants of the emission and absorption in Eqs. 1 and
3, respectively. Since $C^{\text{EM/ABS}}$ and
$\boldsymbol{\hat{\mu}}_{\text{EG}}$ are constants across the spectrum, we
used Eqs. 4 and 5 as the luminescence and absorption lineshapes
[$L^{\text{EM/ABS}}(\hbar\omega)$], respectively, defined as

$L^{\text{EM}}(\hbar\omega)=\omega^{3}\mathop{\raisebox{-20.00024pt}{\scalebox{2.0}{$\displaystyle\sum$}}}_{n,m}c_{m}(T)\Big{|}\Big{\langle}\Theta^{n=0}_{\text{GS}}(Q)\Big{|}\Theta^{m}_{\text{ES}}(Q)\Big{\rangle}\Big{|}^{2}\delta(E^{m}_{\text{ES}}-E^{n}_{\text{GS}}-\hbar\omega)$
(4)

and

$L^{\text{ABS}}(\hbar\omega)=\omega\mathop{\raisebox{-20.00024pt}{\scalebox{2.0}{$\displaystyle\sum$}}}_{n,m}c_{n}(T)\Big{|}\Big{\langle}\Theta^{n}_{\text{GS}}(Q)\Big{|}\Theta^{m}_{\text{ES}}(Q)\Big{\rangle}\Big{|}^{2}\delta(E^{m}_{\text{ES}}-E^{n}_{\text{GS}}-\hbar\omega)\text{.}$
(5)

The vibrational overlap integral can be calculated by applying the HR theory
Huang  _et al._ (1950) as implemented previously Gali  _et al._ (2016). In the
HR framework three basic assumptions are used: ($i$) the normal modes and
($ii$) the vibrational frequencies are identical in the initial and final
electronic states, and ($iii$) the equilibrium configuration is shifted by
$\Delta\mathbf{Y}$ in the final state with respect to the initial nuclear
configuration. Generalized nuclear configuration weighted by the nuclear
masses ($\mathbf{Y}$) can be defined as
$\mathbf{Y}=\mathbf{M}^{\frac{1}{2}}\mathbf{X}$, where
$\mathbf{X}=(\mathbf{R}_{1},...,\mathbf{R}_{N})=(X_{1},Y_{1},Z_{1},...,X_{N},Y_{N},Z_{N})$
is a vector constructed by the nuclear coordinates and

$\mathbf{M}=\begin{bmatrix}M_{1}&&&&&\\\ &M_{1}&&&&\\\ &&M_{1}&&&\\\
&&&M_{2}&&\\\ &&&&\ddots&\\\ &&&&&M_{N}\end{bmatrix}$ (6)

is the mass tensor. The partial HR factor ($S_{k}$) can be calculated as the
scalar product of the normal vectors of mode $k$ ($\mathbf{Y}_{0}^{k}$) and
the displacement vector of $\Delta\mathbf{Y}$ as

$S_{k}=(\Delta\mathbf{Y}^{\text{T}}\mathbf{Y}_{0}^{k})^{2}\text{.}$ (7)

The (total) HR factor ($S$) can be defined as the sum of the partial ones,

$S=\sum_{k}^{3N-6}S_{k}\text{.}$ (8)

While the partial HR factor indicates the weight of phonon mode $k$ in the
corresponding electronic transition, total HR factor gives a measure for the
strength of the overall electron-phonon coupling. By using the HR factors,
overlapping integral of the initial and final vibrational wavefunctions
describing the lineshape of the spectrum can be written as

$\Big{\langle}\Theta^{m}_{\text{f}}(Q)\Big{|}\Theta^{n}_{\text{i}}(Q)\Big{\rangle}=\prod_{k}\langle
m_{k}|n_{k}\rangle,$ (9)

where $\Theta^{m}_{\text{i}}(Q)$ and $\Theta^{n}_{\text{f}}(Q)$ stand for the
vibrational wavefunctions of the initial and final states, respectively
$m_{k}$ and $n_{k}$ stand for the phonon occupation of mode $k$ in the final
and initial states with the corresponding $n$ and $m$ quantum numbers,
respectively, and

$\Big{|}\langle
m_{k}|n_{k}\rangle\Big{|}=e^{-S_{k}}\frac{S_{k}^{m-n}}{(m-n)!}\text{.}$ (10)

At $T=0$ K, the lowest energy phonon state is occupied in the electronic
excited state in the PL process, i.e., $m=0$. As the temperature raised higher
energy $m$ phonon states are occupied. Boltzmann distribution function may be
applied at a given temperature for calculating the occupation of $\omega_{k}$
phonon state which is given by $c_{m}(T)$ in Eqs. 1 and 4. By defining
$\varepsilon$ as the vibronic excited state’s ($m$) energy in the electronic
excited state,

$\varepsilon_{m}=E_{\text{ES}}^{m}-E_{\text{ES}}^{0}\text{,}$ (11)

$c_{m}(T)$ can be given as

$c_{m}=\frac{\exp\Bigl{(}\frac{-\varepsilon_{m}}{k_{B}T}\Bigr{)}}{\sum_{m}\exp\Bigl{(}\frac{-\varepsilon_{m}}{k_{B}T}\Bigr{)}}\text{.}$
(12)

By this way, the temperature dependent PL lineshape may be obtained. In the
case of absorption spectrum, $c_{n}(T)$ can be defined analogously to
$c_{m}(T)$ but it rather represents the occupation of the vibronic excites
state in the electronic ground state. We note that all HR spectra presented in
this work are normalized.

Experimental determination of the $S$ factor can be carried out via the so-
called Debye-Waller (DW) factor ($w$) Debye (1913); Waller (1923) which can be
directly read out from the corresponding spectra. DW factor represents the
ratio of the ZPL and total intensity in the corresponding spectrum and can be
expressed using $S$ as $w=\text{exp}(-S)$. In this work, we present normalized
HR spectra, i.e., probability density functions (PDFs) of the corresponding
electronic transition. By this definition, the area under the curve is one in
all cases, and thus the DW factor is the area under the ZPL peak. We note that
the temperature dependent width of ZPL emission is not incorporated in this
theory and the long wavelength acoustic phonons are not sampled in our
calculations due to the finite size of the supercell.

###  II.2  Charge transition levels

Formation energy as defined in Ref. Csóré  _et al._ , 2021 can be expressed
for the V${}_{\text{C}}$V${}_{\text{Si}}$ defect as

$E_{\text{form}}^{q}=E_{\text{tot}}^{q}-\frac{n_{\text{Si}}+n_{\text{C}}}{2}\mu_{\text{SiC}}+qE_{\text{F}}+\Delta
V(q),$ (13)

where $E_{\text{tot}}^{q}$ is the total energy of the defective system in the
$q$ charge state, $n_{\text{Si/C}}$ is the number of the Si/C atoms in the
supercell, $\mu_{\text{SiC}}$ is the chemical potential of Si-C unit in the
perfect 4H SiC crystal, $E_{\text{F}}$ represents the Fermi-level and $\Delta
V(q)$ stands for the charge correction term. To determine $\Delta V(q)$, we
use the Freysoldt charge correction scheme Freysoldt  _et al._ (2009). For
calculating the $(-/0)$ charge transition level referenced to the conduction
band minimum (CBM) one may use the following formula

$E^{\text{CBM}}_{q+1/q}=E_{\text{CBM}}-E_{q+1/q}\text{,}$ (14)

where $E_{\text{CBM}}$ is the CBM energy, and $E_{q+1/q}$ is the charge
transition level calculated as

$E_{q+1/q}=E_{\text{tot}}^{q}-E_{\text{tot}}^{q+1}+\Delta V(q)-\Delta
V(q+1)\text{.}$ (15)

By increasing the temperature, the charge transition levels may shift due to
the thermal evolution of the band edges and the thermal shift of the
vibrational free energy.

Effect of the thermal shift of the CBM manifests in Eq. 14, where
$E_{\text{CBM}}$ is in fact temperature dependent. Since $E_{\text{CBM}}$
decreases upon increasing the temperature Cannuccia and Gali (2020), the
$(-/0)$ charge transition level referenced to the CBM and hence the
reionization threshold energy will also decrease upon introducing this
correction as $E^{\text{CBM}^{*}}_{-/0}=E^{\text{CBM}}_{-/0}-\Delta
E^{\text{CBM}}(T)$.

On the other hand, thermal shift of the charge transition levels originating
from the vibrational entropy of the defect is governed by the temperature
dependent vibrational free energy in the corresponding charge state ($F^{q}$)
formulated as Wickramaratne  _et al._ (2018)

$F^{q}(T)=\sum_{i}\Big{\\{}\frac{1}{2}\hbar\omega_{i}+k_{\text{B}}T\text{ln}\Big{[}1-\text{exp}\Big{(}-\frac{\hbar\omega_{i}}{k_{\text{B}}T}\Big{)}\Big{]}\Big{\\}},$
(16)

where $\hbar$ and $k_{\text{B}}$ are the reduced Planck constant and Boltzmann
constant, respectively, and $\omega_{i}$ is the frequency of the $i$th phonon
mode. Since free energy contributes to the total energy of the system, the
$E_{-/0}$ charge transition level shifts by the free energy difference of the
negative and neutral charge states [$\Delta F(T)$]. The free energy difference
is defined as $\Delta F(T)=F^{-}(T)-F^{0}(T)$, so the free energy-shifted
charge transition level can be written as
$E^{\text{CBM}^{**}}_{-/0}=E^{\text{CBM}}_{-/0}+\Delta F(T)$. In order to
calculate $\Delta F(T)$, we considered all the quasilocal phonon modes
associated with the motion of the atoms within the second neighbor shells
around the V${}_{\text{C}}$V${}_{\text{Si}}$ defect both in the negative and
neutral charge states.

Combining the two correction terms, the reionization energy threshold energy
is defined as

$E^{\text{CBM}}_{-/0}(\text{corr})=E^{\text{CBM}}_{-/0}-\Delta
E^{\text{CBM}}(T)+\Delta F(T)\text{.}$ (17)

###  II.3  Computational methodology

In order to obtain the neutral PL spectra for all defect configurations, we
calculated the neutral ground and optically accessible excited states by DFT.
Excited state calculations were carried out by applying the $\Delta$SCF method
Gali  _et al._ (2009). In addition, we calculated the ground state of
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ for the
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ reionization spectra. Nuclear
coordinates of the global energy minimum in the initial and final electronic
states, i.e., the fully relaxed geometries were obtained by minimizing the
quantum mechanical forces between the ions falling below the threshold of 0.01
eV/Å. The corresponding electronic structure for all
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defect configurations in the neutral
ground/excited state and that for V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ in the
ground state is shown in Fig. 2. Here we note that both excited state of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ and ground state of
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ introduce a degenerate $e$ KS state to
the band gap occupied by three electrons giving rise to the so-called
dynamical Jahn-Teller (JT) effect Bersuker (2006). Accordingly, the symmetry
will be reduced from C${}_{\text{3v}}$ to C${}_{\text{1h}}$ (even for axial
configurations) by coupling to phonons in order to split the corresponding $e$
level removing the degeneracy in total energy. In calculation of the
corresponding spectra we use static JT distorted geometries exhibiting the
lowest total energy.

Figure 2: Electronic structure of the [(a), (d)] ground and [(b), (e)] excited
state of V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defects and [(c), (f)] ground
state of V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ defects. Figures (a)-(c)
corresponds to the axial defects exhibitng C${}_{\text{3v}}$ symmetry and
those of (d)-(f) show the electronic structures of the basal defect
configurations possessing C${}_{\text{1h}}$ symmetry. Black arrows represent
electrons while gray arrows stand for the reminiscent holes after optical
excitation. Valence band (VB) and conduction band (CB) of 4H SiC are also
indicated along with the characters of the corresponding Kohn-Sham states.

All investigated V${}_{\text{C}}$V${}_{\text{Si}}$ defect configurations were
modelled in a 576-atom 4H supercell. In order to reach sufficient accuracy in
the total energies, we employed $\Gamma$ – point centered 2 $\times$ 2
$\times$ 2 $k$-point mesh for all ground and excited state calculations. Here
we note that only $\Gamma$-point sampling of the Brillouin-zone fails to
achieve sufficient accuracy manifesting in the incorrect order in the total
energies of the different V${}_{\text{C}}$V${}_{\text{Si}}$ configurations as
pointed out in Ref. Magnusson  _et al._ , 2018. For the electronic structure
calculations we employed the HSE06 range-separated hybrid functional Heyd  _et
al._ (2003). Kohn-Sham (KS) wavefunctions were expanded in plane wave basis
set with the cutoff energy of 420 eV. Only the valence electrons were treated
explicitly, the core-electrons with ionic potentials were considered in the
framework of projector augmented wave (PAW) method Blöchl (1994) as
implemented in the Vienna Ab-Initio Simulation Package (VASP) Kresse and
Furthmüller (1996).

For the vibrational modes, we calculated the corresponding dynamical matrix
containing the second order derivatives of the total energy by means of the
Perdew-Burke-Ernzerhof (PBE) Perdew  _et al._ (1996) functional. We also
determined the corresponding eigenvectors. Here we note that the calculated
dynamical matrix and hence the HR and DW factors depend on the supercell size
where the former may be a bit underestimated because of neglect of long
wavelength phonons.

Localization of the phonon modes can be addressed by the inverse participation
ratio (IPR) Alkauskas  _et al._ (2014); Bell  _et al._ (1970) defined as

$\text{IPR}_{k}=\frac{N\sum_{i}\mathbf{u}_{i}^{4}}{\big{(}\sum\mathbf{u}_{i}^{2}\big{)}^{2}},$
(18)

where $\mathbf{u}_{i}$ is the vector displacement amplitude of the $i$th atom
in the $k$th phonon mode and $N$ is the number of atoms in the defective
supercell that is $N$=574 in our case. Consequently, the IPR falls in the
region of $[1,N]$ and it is equal to the number of atoms vibrating in the
certain phonon mode.

##  III Results and Discussion

In Secs. III.1 \- III.3 we provide our numerical results for the
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ PL spectra and the temperature dependent
reionization spectra for the V${}_{\text{C}}$V${}_{\text{Si}}^{-}\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ transition for all
V${}_{\text{C}}$V${}_{\text{Si}}$ defect configurations. We also discuss the
charge state control and simultaneous qubit manipulation of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ in 4H SiC with comparing our novel
experimental and simulation results in Sec. III.4.

###  III.1 PL spectra of V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ configurations

In order to verify the methodology for calculating the PL spectra Gali  _et
al._ (2009, 2016), we first compare the total PL spectrum (i.e., including all
contributions from the individual V${}_{\text{C}}$V${}_{\text{Si}}$
configurations) obtained from theory with the non-quenched
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ spectrum reported in Sec. III.4. The
non-quenched PL spectra were recorded at $T=3.5$ K with using two different
excitation wavelengths, 976 nm and 1053 nm (c.f., Fig. 10) along with a 532-nm
repump laser in order to maintain the neutral charge state (i.e., to avoid
quenching). We note that our simulations correspond to the $T=0$ K condition
which is close to the experimental one. Excitation energies can be taken into
account in the calculations by multiplying the obtained numerical spectra with
the optical absorption cross-sections ($\sigma$) of the individual
V${}_{\text{C}}$V${}_{\text{Si}}$ configurations, respectively, as
$L^{\text{EM}}(E_{\text{ex}},\hbar\omega)=\sigma(E_{\text{ex}})L^{\text{EM}}(\hbar\omega)$,
where $L^{\text{EM}}(E_{\text{ex}},\hbar\omega)$ denotes the PL lineshape upon
the excitation energy of $E_{\text{ex}}$. In our calculations, we identified
$\sigma(E_{\text{ex}})$ as the transition probability at $E_{\text{ex}}$ in
the absorption PDF, i.e.
$\sigma(E_{\text{ex}})=\text{PDF}_{\text{abs}}(E_{\text{ex}})$. We report the
corresponding values in Table 1.

Table 1: Absorption cross-sections ($\sigma_{\text{th}}$) and ratios of the ZPL intensities ($R_{\text{th}}$) in the numerical PL spectrum for the individual V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ configurations. We also report the ZPL ratios for the experimental non-quenched spectra recorded at $T=3.5$ K ($R_{\text{exp}}$). We show values for both 976-nm and 1053-nm excitations.  Config. | $\sigma_{\text{th}}$ | $R_{\text{th}}$ | $R_{\text{exp}}$
---|---|---|---
976 nm | 1053 nm | 976 nm | 1053 nm | 976 nm | 1053 nm
$hh$ | 0.0015 | 0.0013 | 1.00 | 1.00 | 1.00 | 1.00
$kk$ | 0.0018 | 0.0017 | 1.26 | 1.33 | 1.16 | 1.33
$hk$ | 0.0021 | 0.0019 | 1.42 | 1.42 | 2.62 | 2.55
$kh$ | 0.0022 | 0.0015 | 1.53 | 1.18 | 3.90 | 1.75

The total numerical PL emission spectra (cf. Fig. 3) upon 976-nm and 1053-nm
excitation are then obtained by adding up the PL spectra obtained for the
individual V${}_{\text{C}}$V${}_{\text{Si}}$ configurations weighted by the
relative strength of absorption probability ($\sigma_{\text{th}}$) at those
excitation wavelengths. Numerical and experimental emission spectra can be
well-compared via the ZPL intensities. To this end, we listed the ratios of
the numerical and experimental ZPL intesities in Table 1. We find that the
trend in the ZPL intensity ratios agrees well; however, there are differences
between the absolute values. This discrepancy may arise from the fact that,
while polarization properties manifest in the experimental spectrum (i.e., the
experimental spectrum is the sum of the individual polarization contributions
parallel and perpendicular to the $c$-axis), all the emitted photons are
considered in the theoretical spectrum. This implies that optical alignment of
the PL measurement plays a key role in the recorded lineshape. In order to
allow detection of both polarizations simultaneously, the experimental PL
spectrum was registered through the edge of the sample as given in Sec. III.4.
Thus, the origin of the discrepancy may be caused by other processes that are
not considered in our model.

Figure 3: Calculated PL spectra of divacancies in 4H SiC upon (a) 976-nm and (b) 1053-nm excitation wavelengths. Figure 4: Calculated Huang-Rhys emission spectrum for the optical transition of V${}_{\text{Si}}$V${}_{\text{C}}^{0}$ for all defect configurations (see textbox) assuming resonant excitation. All the spectra are normalized, i.e., the area under the curves is one, hence, the $y$-axis represents the electronic transition probability from the optically accessible excited state to the ground state. Table 2: Calculated and experimental (Ref. Falk  _et al._ , 2013) ZPLs for all V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defect configurations in 4H SiC. We apply the notations for the centers observed in experiments from Ref. Koehl  _et al._ , 2011. The calculated Huang-Rhys ($S$) and Debye-Waller ($w$) factors are also listed. Config. (Center) |  $E_{\text{ZPL}}^{\text{exp}}$ (eV) |  $E_{\text{ZPL}}^{\text{calc}}$ (eV) | $S$ | $w$
---|---|---|---|---
$hh$ (PL1) | 1.095 | 1.142 | 2.856 | 0.057
$kk$ (PL2) | 1.096 | 1.151 | 3.380 | 0.034
$hk$ (PL3) | 1.119 | 1.167 | 3.428 | 0.032
$kh$ (PL4) | 1.150 | 1.187 | 3.540 | 0.029

The calculated HR emission PDFs of V${}_{\text{C}}$V${}_{\text{Si}}^{0}$
configurations are shown in Fig. 4. ZPLs in all spectra are shifted to the
experimental ZPL positions as reported in Table 2. We note that the
fluorescence spectra of V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ $hh$ and $hk$
configurations have recently been computed in Ref. Hashemi  _et al._ , 2021
reporting similar lineshapes. Here we provide the PL lineshape for each
configuration with high resolution in the PSB for further analysis. IPR
analyis of the phonons shows that phonon modes at $\approx 17-35$ meV are
mostly localized on the three Si atoms neighboring the V${}_{\text{C}}$ with
the IPR of $\approx 1-5$ while phonon modes at $\approx 35-115$ meV region are
mostly localized on the three C atoms around the V${}_{\text{Si}}$ with a
similar IPR range [Fig. 5(b)]. The sum of the individual PL spectra is shown
in Fig. 5(a).

The corresponding HR and DW factors along with the experimental and calculated
ZPLs are listed in Table 2. The lowest ZPL energy corresponds to
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ along with the smallest HR factor
and thus the largest DW factor. Here we note that the exceptionally large DW
factor of the V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ was also found in a
previous study Hashemi  _et al._ (2021). Interestingly, ZPL energy of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(kk)$ is extremely close to that of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$, i.e., the difference between the
experimental values is only 1 meV ($\approx 0.1$ %). On the other hand,
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(kk)$ exhibits by $\approx 19$ % larger HR
factor and $\approx 40$ % lower DW factor than those of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$. We note that the DW factor of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ is close to that of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ in the cubic 3C SiC (0.07) Csóré  _et
al._ (2022). We already showed a counterintuitive issue on the example of
vanadium defect in 4H SiC Spindlberger  _et al._ (2019) that the hexagonal
site in 4H SiC has a cubic-like environment in the second and third neighbor
shells unlike the quasicubic site. Therefore, $hh$ configuration of divacancy
in 4H SiC indeed shows similar properties to those of the divacancy in 3C SiC.
The largest ZPL energy and HR factor along with the lowest DW factor occur for
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(kh)$ in 4H SiC.

Figure 5: (a) Probability density function (PDF) of the
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ defects in 4H SiC related to (a) total
and (b) individual PL spectra. The total spectrum is generated by adding up
the individual PDFs of V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ configurations.
ZPL energies are shifted to $0$ meV, in order to read the phonon energy.
Localization of phonon modes is indicated by shaded area in (b) caused by the
three Si atoms around the V${}_{\text{C}}$ (yellow shaded area) and by the
three C atoms around V${}_{\text{Si}}$ (blue shaded area).

###  III.2 Temperature dependent reionization spectra from phonon-assisted
transitions

We conclude that the HR theory V${}_{\text{C}}$V${}_{\text{Si}}^{0}$
reproduces well the experimental spectrum, therefore we can apply this theory
for understanding the temperature dependent reionization process. Ionization
energy of point defects can be approximated by their charge transition level
which is the position of the Fermi-level in the band gap for which the
formation energies of the defect in the two charge states are equal. We showed
in our previous study Magnusson  _et al._ (2018) that the energy separation
between the $(+/0)$ level and the valence band maximum (VBM) is at
$E_{+/0}^{\text{VBM}}\approx 1.1$ eV, which is smaller than that between the
$(0/-)$ level and the CBM at $E_{-/0}^{\text{CBM}}\approx 1.2-1.3$ eV for all
V${}_{\text{C}}$V${}_{\text{Si}}$ configurations in 4H SiC at $T=0$ K.
Therefore we concluded that the dark state of the divacancy is
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$, i.e., the singly negative charge state
of the defect, into which V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ is converted
upon photoexcitation below the reionization threshold energy because the
photon energy of the applied laser used to illuminate
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ ($>1.1$ eV) could immediately convert
V${}_{\text{C}}$V${}_{\text{Si}}^{+}$ into
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ with ejecting a hole in the valence
band. The experimental threshold energies recorded at cryogenic temperatures
and the calculated charge transition levels are listed in Table 3.

Table 3: Experimental threshold energies of photoionization recorded at cryogenic temperatures (Ref. Magnusson  _et al._ , 2018) which is approximated by the $(0/-)$ charge transition levels for all V${}_{\text{C}}$V${}_{\text{Si}}$ configurations referenced to the CBM. The corresponding HR ($S$) and DW ($w$) factors are also provided which are related to the phonon assisted ionization process. Center | Config. |  $E_{\text{rei}}$ (eV) |  $E^{\text{CBM}}_{0/-}$ (eV) | $S$ | $w$
---|---|---|---|---|---
PL1 | $hh$ | 1.310 | 1.245 | 2.691 | 0.068
PL2 | $kk$ | 1.310 | 1.209 | 3.038 | 0.048
PL3 | $hk$ | 1.321 | 1.307 | 2.490 | 0.083
PL4 | $kh$ | 1.281 | 1.174 | 3.267 | 0.038

The reionization threshold energies can depend on the temperature that we
divided to three contributors: (i) phonon-assisted photoionization, (ii)
entropy contribution to the charge transition levels, (iii) temperature shifts
of the band edges of the host crystal.

First, we consider here the phonon-assisted photoionization process. As we are
interested in the threshold energy, we consider only the CBM and the phonon-
assisted optical transition to CBM which is technically well represesented in
our simulation in terms of folding the M-point of the Brillouin-zone to the
$\Gamma$-point of the supercell. The calculated reionization threshold energy,
i.e., the $E^{\text{CBM}}_{0/-}$ in the particular case of divacancies in 4H
SiC, can be considered as an absorption of the electron from the highest
occupied defect level by the CBM states which corresponds to a no-phonon
process. At elevated temperatures, the vibronic excitated states of
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ may be occupied which effectively lowers
the threshold energy to reach CBM compared to that of the no-phonon process.
This is analogous to the ZPL energy vs. phonon sideband related energies in
the absorption process, so HR absorption spectrum calculations were applied
for all V${}_{\text{C}}$V${}_{\text{Si}}$ configurations as shown in Fig. 6.
We used the V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ phonons for calculating the
phonon sideband of the reionization spectra as we did for the PL spectra (see
Sec. III.1).

Figure 6: Calculated temperature dependent transition probability density
functions for HR absorption spectrum of V${}_{\text{C}}$V${}_{\text{Si}}^{-}$
$\rightarrow$ V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ transition for the four
configurations as indicated in the text boxes at $0$ K (black lines), $150$ K
(red lines), $300$ K (blue lines), $450$ K (green lines) and $600$ K (magenta
lines). The ZPL energies - associated with the reionization energies at
cryogenic temperature - were shifted to the experimental data recorded at
cryogenic temperatures (Ref. Magnusson  _et al._ , 2018 as listed in Table 3).

In Fig. 6, the ZPL peak corresponds to the no-phonon ionization energies and
coincides with the adiabatic $E^{\text{CBM}}_{0/-}$ charge transition levels.
By raising the temperature, peaks and broad features appear at lower energies
than the ZPL energy which are responsible to lower the reionization threshold
energies. There are some subtle differences between in the low-energy PSB
spectrum of the different configurations, e.g., sharp features are visible for
$hh$ and $hk$ configurations at $T=600$ K, while those for $kk$ and $kh$
configurations almost completely disappear. This implies a stronger electron-
phonon coupling in the reionization process for configurations with
V${}_{\text{C}}(h)$ which are associated with the vibration of the three
nearest neighbor Si atoms. The HR and DW factors of the absorption
(reionization) spectra for each configuration are listed in Table 3 which
characterize the general strength of the electron-phonon coupling associated
with the optical transition.

The calculation of the phonon-assisted reionization process can be expressed
as $E_{\text{rei}}=E_{\text{ex}}\pm E_{\text{ph}}$, where $E_{\text{ex}}$
represents the applied photoexcitation energy and $E_{\text{ph}}$ stands for
the energy of the absorbed ($-$) or emitted ($+$) phonon. The former process
will reduce the threshold energy for reionization. At croygenic temperatures
close to $T=0$ K, the no-phonon process is dominating for the threshold energy
of reionization, i.e., $E_{\text{ex}}^{\text{min}}=E_{\text{rei}}$ as listed
in Table 3. At elevated temperatures, Fig. 6 clearly shows that, indeed, the
threshold energy for reionization may be reduced but the absorption cross-
section has an exponential tail in the low-energy PSB of the optical
transition. However, it is difficult to define the exact condition at which
absorption cross-section the reionization process is observable in
experiments. In order to quantify the theoretical description of this process,
we setup a condition in which we assume that the reionization process is
definitely observable when the total reionization cross-section mediated by
the phonons at a given temperature is at least as large as the integrated
reionization cross-section associated with the half of the DW factor of the
reionization spectrum which corresponds to the reionization cross-section at
cryogenic (close to $T=0$ K) temperatures.

Since the calculated reionization spectra are normalized (those are PDFs),
exciting with the energy of $E_{\text{ex}}$ results in the total transition
probability , i.e., the cumulative transition probability (CTP) of

$\text{CTP}(E_{\text{ex}})=\int_{0}^{E_{\text{ex}}}\text{PDF}(\varepsilon)\mathrm{d}\varepsilon,$
(19)

where $\varepsilon$ is the variable for the excitation energy in the
reionization spectra. At cryogenic temperatures, we assume that
CTP($E_{\text{rei}}=E_{\text{ZPL}}$) = $w/2$ with assuming symmetrical
broadening of the ZPL peak of the reionization spectrum. By increasing the
temperature, the phonon excited states of the negative divacancy are occupied
according to the Bose-Einstein function, thus the effective reionization
energy decreases which results in a lower threshold energy of photoionization.
In other words, the phonon-assisted photoionization results in
PDF$(\varepsilon)>0$ for $\varepsilon<E_{\text{ZPL}}$ energy. One can find the
threshold reionization energy at a given temperature for each configuration
for which the reionization probability is equal to that of at $T=0$ K, i.e.,
PDF$(\varepsilon)=w/2$ is fulfilled.

Figure 7: (a) Numerical temperature dependent transition probability density
functions of the HR absorption (reionization) spectra for the
V${}_{\text{C}}$V${}_{\text{Si}}^{-}(hh)$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ transition. PSB below the ZPLs given
rise at elevated temperatures are enlarged and the corresponding reionization
threshold energies ($E_{\text{ex}}^{\text{min}}$) are indicated by vertical
dashed lines, whereas the cumulative transition probability (CTP) of $w/2$ are
represented by shaded areas. (b) CTP for the
V${}_{\text{C}}$V${}_{\text{Si}}^{-}(hh)$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ transition. Half of the DW factor,
i.e., $w/2$ is indicated in the expanded part of the CTP.

We calculated the CTP($\varepsilon$) at temperatures of
$T=0\leavevmode\nobreak\ \text{K},150\leavevmode\nobreak\
\text{K},300\leavevmode\nobreak\ \text{K},450\leavevmode\nobreak\
\text{K},600\leavevmode\nobreak\ \text{K}$ for each divacancy configuration in
4H SiC that we depict for $(hh)$ configuration in Fig. 7. CTP($\varepsilon$)
can be read out as shaded area under the curves in Fig. 7(a) whereas it is a
line plot in Fig. 7(b). CTP($\varepsilon$) at all investigated temperatures is
similar to the lineshape of the cumulative distribution function of the
Boltzmann distribution. The wavy lineshapes around the ZPL energy visible for
low-temperature simulations could be an artifact by neglecting the temperature
broadening of the ZPL emission and other size effects of the supercell
approach.

We also indicated $E_{\text{ex}}^{\text{min}}$ at each temperature in Fig.
7(a) and depicted for all V${}_{\text{C}}$V${}_{\text{Si}}$ configurations at
$T=0\dots 600$ K in $50$ K steps in Fig. 8(a). This approach results in a
plateau up to $T\approx 150-200$ K for $E_{\text{ex}}^{\text{min}}(T)$, and it
then drops for all the V${}_{\text{C}}$V${}_{\text{Si}}$ configurations at
higher temperatures. In other words, E${}_{\text{ph}}^{\text{min}}(T)$ is
about zero in our approximation for $T<150-200$ K, and it increases for
$T>150-200$ K. In order to gain further insight to the temperature dependence
of the phonon threshold energy, E${}_{\text{ph}}^{\text{min}}(T)$, we employed
the Boltzmann distribution [$f(\varepsilon,T)$] as the low concentration limit
of the Bose-Einstein distribution for the phonon bath appearing in the region
of [0, $E_{\text{rei}}$] as

$f(\varepsilon,T)=\frac{1}{kT}\cdot\text{exp}\Big{(}\frac{\varepsilon-
E_{\text{rei}}}{kT}\Big{)},$ (20)

where $k$=8.617 $\cdot$ 10-5 eV/K is the Boltzmann-constant. By setting the
value of the integral to $w/2$ in the region of [0, $E_{\text{rei}}$], we
arrive at

$\int_{0}^{E_{\text{ex}}^{\text{min}}}\frac{1}{kT}\cdot\text{exp}\Bigg{(}\frac{\varepsilon-
E_{\text{rei}}}{kT}\Bigg{)}\mathrm{d}\varepsilon=\frac{w}{2}$ (21)

By evaluating the integral of Eq. 21 and rearranging the equation, we arrive
at

$E_{\text{ph}}^{\text{min}}(T)=-\text{ln}\Bigg{[}\frac{w}{2}+\text{exp}\Bigg{(}\frac{E_{\text{rei}}}{kT}\Bigg{)}\Bigg{]}\cdot
kT\text{.}$ (22)

Since the values of $E_{\text{rei}}$ are around 1300 meV the typical values of
$\text{exp}\big{(}\frac{E_{\text{rei}}}{kT}\big{)}$ is about $10^{-12}$ at the
given temperature range, and hence negligible next to the calculated DW values
of $w/2\approx 0.03-0.09$ (see Table 3). Finally,
$E_{\text{ph}}^{\text{min}}(T)$ can be approximated as

$E_{\text{ph}}^{\text{min}}(T)\approx-\text{ln}\Bigg{(}\frac{w}{2}\Bigg{)}\cdot
kT$ (23)

yielding linear characteristics for $E_{\text{ph}}^{\text{min}}(T)$ with the
slope of $m=-\text{ln}(w/2)\cdot k$. Furthermore, since $w/2<1$,
$E_{\text{ph}}^{\text{min}}(T)$ increases with $T$ within this approximation.
As a result, reionization threshold energy can be expressed as

$E_{\text{ex}}^{\text{min}}(T)=E_{\text{rei}}-E_{\text{ph}}^{\text{min}}(T)\approx
E_{\text{rei}}+\text{ln}\Bigg{(}\frac{w}{2}\Bigg{)}\cdot kT\text{.}$ (24)

Linear behavior of $E_{\text{ex}}^{\text{min}}(T)$ is indeed manifested in our
simulations plotted in Fig. 8(a), for $T>200$ K. This approach predicts that
$E_{\text{rei}}(T=300\leavevmode\nobreak\ \text{K})$ reaches 1250-1320 meV
(990-940 nm) for V${}_{\text{C}}$V${}_{\text{Si}}$ configurations in 4H SiC
which is close to the typical photoexcitation energy of the neutral
divacancies in 4H SiC.

Figure 8: Reionization threshold energies [$E_{\text{ex}}^{\text{min}}(T)$]
(a) without thermal corrections on the charge transition levels, (b) corrected
with $\Delta F(T)$ and (c) corrected further with $\Delta E^{\text{CBM}}(T)$
as a function of temperature for all V${}_{\text{C}}$V${}_{\text{Si}}$ defect
configurations.

We note that Eq. 24 is valid for any point defects and their PL spectrum
within the employed approximations by substituting $E_{\text{rei}}$ with the
corresponding ZPL energy of (re)ionization. We emphasize that the actual
phonon-assisted reionization may be observed for $\varepsilon$ with
PDF$(\varepsilon)<w/2$ depending on the integration time and other factors in
the experiments that are discussed in Sec. III.4. Furthermore, the temperature
shifts of the charge transition level and the band edges, i.e., temperature
dependence of $E_{\text{rei}}$, should be also taken into for understanding
temperature dependence of the reionization process.

###  III.3  Temperature dependent reionization thershold energies including
contributions from entropy and temperature shifts of the conduction band
minimum

Temperature dependence of the charge transition levels of point defects was
discussed in Sec. II.2 which depends on the entropy effects in the formation
enthalpy of the point defects in their charge states ($\Delta F(T)$) and the
electron-phonon renormalization of the band edges. Thermal evolution of the
band edges in 4H SiC was already determined in a previous study Cannuccia and
Gali (2020). According to this work, thermal shifts of the VBM and CBM becomes
noticeble at around $T=300$ K, where the CBM/VBM is decreased/increased by
$\approx 5$ meV, respectively, with respect to their $T=0$ K value. In Table
4, we provide the values for $\Delta F(T)$ and
$E^{\text{CBM}}_{-/0}(\text{corr})$ for all V${}_{\text{C}}$V${}_{\text{Si}}$
defect configurations at room-temperature.

Table 4: Free energy correction term ($\Delta F$) to $E^{\text{CBM}}_{-/0}$ and the corrected charge transition levels referenced to the CBM [$E^{\text{CBM}}_{-/0}(\text{corr})$] at room temperature ($T=300$ K). We already included the $\Delta E_{\text{CBM}}(T)=5$ meV decrease of the CBM in the calculation of $E^{\text{CBM}}_{-/0}(\text{corr})$.  Config. | $hh$ | $kk$ | $hk$ | $kh$
---|---|---|---|---
$\Delta F$(eV) | 0.041 | -0.008 | -0.001 | 0.050
$E^{\text{CBM}}_{-/0}(\text{corr})$ (eV) | 1.281 | 1.196 | 1.301 | 1.219

We find that the $\Delta F$ is one order magnitude lower for the $kk$ and $hk$
configurations than those for the $hh$ and $kh$ configurations. Since
$F^{q}(T)$ depends on frequency of the phonon modes (cf. Eq. 16),
$\omega_{i}$, this trend may arise from the stronger distortion between the
neutral and negative charge states for the $hh$ and $kh$ configurations than
that for the $kk$ and $hk$ configurations. Indeed, the geometry distortion is
larger by about 10% for the $hh$ and $kh$ configurations than that for $kk$
and $hk$ configurations. As a result, correction in $E^{\text{CBM}}_{-/0}$
(cf. Table 3) for the $kk$ and $hk$ configurations is tiny and further reduces
the reionization threshold energies. On the other hand, the free energy
contribution rather inceases the reionization threshold energies for $hh$ and
$kh$ configurations and mostly compensates the thermal shift of the band edge
[cf. Figs. 8(b) and (c)]. These results clearly demonstrate that considering
only the free energy term in the temperature dependence of the ionization
threshold energies can lead to a false result as the contribution from the
phonon-assisted ionization process could dominate and even change the sign of
the thermal shift in the ionization threshold energy [cf. Fig. 8(a) and (c)].

###  III.4  Experimental aspects

The above consideration shows that the quenching behavior of the divacancy in
all configurations is distinctly different at low temperature (say, < 2 K) and
at elevated temperatures. While at low temperature there is a distinct
threshold for the excitation energy for each divacancy configuration, below
which quenching of the corresponding PL is observed but can be recovered using
higher-energy repump illumination. Above the threshold the PL of the
corresponding divacancy configuration does not quench even if no repump is
applied. At higher temperatures, however, there are no distinct thresholds
because phonon-assisted photoionization of the VV- charge state becomes
possible, whereas no phonons are available at low temperature for this
process. Given a fixed excitation energy $E_{\text{ex}}$, the minimum energy
of a phonon $E_{\text{ph}}^{\text{min}}$ that can combine with the incident
photon to produce phonon-assisted photoionization of VV- is given by the
energy balance quoted in Sec. III.1,

$E_{\text{rei}}=E_{\text{ex}}+E_{\text{ph}}^{\text{min}}$ (25)

Here we assume that $E_{\text{ex}}<E_{\text{rei}}$, hence
$E_{\text{ph}}^{\text{min}}>0$. All phonons with energies
$E_{\text{ph}}>E_{\text{ph}}^{\text{min}}$ also contribute to phonon-assisted
photoionization, hence the cumulative transition probability
CTP($E_{\text{ex}}$) is given by Eq. 19. From practical point of view, it may
be desirable to have the possibility for optical charge state control also at
elevated temperatures. This can be achieved by using lower excitation energies
$E_{\text{ex}}$, because these provide lower CTP($E_{\text{ex}}$), hence the
quenching behaviour can be preserved also at elevated temperatures. It is
clear that if $E_{\text{ph}}^{\text{min}}\gg k_{\text{B}}T$ for given
temperature $T$ (the case of lower-energy E${}_{\text{ex}}$) the amount of
phonons available for phonon-assisted reionization is negligible and
CTP($E_{\text{ex}}$) $\approx$ 0\. Note that $E_{\text{ex}}$ is bound from
below by the condition $E_{\text{ex}}>E_{\text{ex}}^{\text{ZPL}}$, where
$E_{\text{ex}}^{\text{ZPL}}$ refers to the corresponding divacancy
configuration. Our goal here is to examine the concept of phonon-assisted
reionization by comparing the quenching properties at higher temperatures for
two excitation energies, both below the thresholds, but one of which is close
to the threshold while the other is far enough so that
$E_{\text{ph}}^{\text{min}}=E_{\text{rei}}-E_{\text{ex}}\gg k_{\text{B}}T$.

Figure 9: Theoretical dependencies of the cumulative transition probability on
the temperature for the four divacancy configurations as denoted for each
curve taking all the temperature effects into account, and for the two
excitation wavelengths used in the experiment. The inset displays an energy
diagram of the phonon-assisted photoionization illustrated on the hk divacancy
configuration. The $(0/-)$ charge transition levels are labelled with the
corresponding configuration and the experimentally determined energy
separation from the conduction band. The two bold vertical arrows represent
the two experimental laser energies whereas the wavy arrows represent the
phonons with minimum energies needed for reionization of the divacancy from
negative to neutral charge state.

The energy diagram (inset in Fig. 9) depicts the two excitation cases for
$E_{\text{ex}}=1.270$ eV ($976$ nm), and $E_{\text{ex}}=1.177$ eV ($1053$ nm)
which have been used in the experiments presented below. Both excitations are
below the thresholds for all divacancy configurations, but the former is close
to the thresholds (only $\approx 11$ meV lower than the threshold for PL4, cf.
Table 2), whereas the $1.177$ eV excitation is about $100$ meV (up to $140$
meV for the different configurations) below all thresholds. As a consequence,
the condition $E_{\text{ph}}^{\text{min}}\gg k_{\text{B}}T$ is fulfilled only
for the second excitation at $1.177$ eV, as long as the temperature does not
exceed room temperature. The CTPs as a function of temperature for the two
used excitations and for the four divacancy configurations are calculated
using Eq. 19 and displayed in Fig. 9 in the temperature range $0–350$ K. Using
these dependencies we can make qualitative comparison with the experimental
data presented in Fig. 10 and obtained using these two fixed laser
excitations, $976$ and $1053$ nm. The repump laser beam used in these
experiments is at $532$ nm directed unfocused on the sample with power density
of the order of $1$ mW/cm-2. On the other hand, the excitation beam ($976$ or
$1053$ nm) is moderately focused on the sample to a spot of $\sim$1 mm
diameter, with estimated power density at the sample of the order of several
watts per cm-2. The experimental conditions are detailed in Ref. Magnusson
_et al._ , 2018.

Figure 10: Experimental spectra of the divacancy emission obtained at
different temperatures with the simultaneous application of the IR laser and
the $532$ nm repump laser (denoted NQ, or non-quenched) and with the IR laser
alone after a long term quenching ($Q$, or quenched). The panels to the left
(right) refer to the $976$ nm ($1053$ nm) IR laser, respectively. The two
bottom panels show also the anti-Stokes part of the spectra and illustrate the
upconversion in the emission of the silicon vacancy.

We compare now the quenching behaviour of the divacancy emission at different
temperatures when the two different excitations are used. The experimental
spectra obtained at different temperatures are presented in Fig. 10. The
spectra with the infrared laser excitation alone ($976$ or $1053$ nm) are
obtained after switching off the repump laser and quenching the divacancy PL
for about 1/2 hour, and are referred to further as the “quenched” spectra
(denoted by $Q$ in Fig. 10). The rest of the spectra are with the repump laser
applied in addition to the infrared excitation; under these conditions the
divacancy PL is stable and the spectra are denoted by NQ for “non-quenched” in
Fig. 10. One can see from the spectra that at liquid He temperature ($3.5$ K)
the divacancy PL quenches about $10–30$ times for each of the two excitations,
as expected because both excitation energies are below the thresholds for all
divacancy configurations. Also, at this temperature essentially no phonons are
available for phonon-assisted reionization, hence the accumulation of the
divacancy in the negative charge state is irreversible for both IR excitations
used. We notice also that in the particular sample used in these experiments
the quenched PL never reaches zero intensity, similar to sample 1 of Ref.
Magnusson  _et al._ , 2018.

At 80 K, however, we notice differences in the quenching behavior of one of
the divacancy configurations ($kh$). Namely, with the 976-nm excitation the
PL4 line corresponding to the $kh$ configuration shows nearly the same
intensity in the “quenched” spectrum as in the “non-quenched”, i.e., it does
not quench. We notice here that the relative contribution of the PL4 line in
the spectra at $80$ K is strongly reduced compared to the spectra at $3.5$ K,
but this behavior is not understood at present. The rest of the lines
(PL1-PL3) do quench exhibiting large intensity contrast (similar to that at
$3.5$ K) between the “quenched” and “non-quenched” spectra for both
excitations. We notice further that the PL4 also quenches in a way similar to
that at $3.5$ K if the 1053-nm excitation is used. Referring to the CTP plots
in Fig. 9 we notice that at temperature $T=80$ K only the $kh$ configuration
(corresponding to the PL4 line) exhibits CTP$\approx 0.0015$ for 976-nm
excitation, while for the rest of the divacancy configurations at $976$ nm and
all configurations at $1053$ nm the CTP is essentially negligible ($<3.5\times
10^{-4}$) for both excitation wavelengths. Thus, the experimental data for the
PL4 line implies that for CTP of only 0.0015 (0.15% probability of phonon-
assisted ionization of V${}_{\text{C}}$V${}_{\text{Si}}^{-}$) the reionization
process outcompetes the accumulation of the divacancy in the negative charge
state. Consequently, for CTPs in the low 10-3 range one may expect to notice
impact of the phonon-assisted reionization on the quenching process. This
provides a means to define semi-quantitatively the experimental conditions
when the rate of the phonon assisted reionization is comparable with the rate
of capture of electrons to neutral divacancies, thus creating negatively-
charged ones. We can define that this condition (similar rates of ionization
an reionization) is true if the quenched spectra exhibit roughly 50% of the
intensity of the non-quenched ones.

At temperatures higher than $80$ K the contribution of the $kh$ configuration
(PL4 line and phonon sideband) in the spectrum can be neglected owing to the
abovementioned rapid decrease of the PL4 intensity with increasing
temperature. Thus, we investigate the behavior of the PL1-PL3 lines, the
phonon sidebands of which coalesce into a broad band extending between $1100$
and $1320$ nm. Note that also the threshold energies of these three
configurations ($hk$, $hh$ and $kk$) are rather close according to the energy
diagram inset in Fig. 9. At $T=150$ K, the contributions from the these three
divacancy configurations are still discernible via their zero-phonon lines,
although the latter are much broader and weaker than at cryogenic temperature.
However, the dominant contribution to the spectrum at this and higher
temperatures comes from the PSBs of each configuration, which overlap to form
a single band concealing the individual contribution from the different
configurations. Therefore, in the following we shall compare the intensities
of the broad bands in the “quenched” and “non-quenched” spectra, looking for
each infared (IR) excitation for the temperature when the intensity in the
quenched spectrum is roughly half of the intensity in the non-quenched one.

We observe from Fig. 10 that our criterion is fulfilled for the 976-nm
excitation at $T=150$ K, when the intensity of the phonon sideband is about
1/2 of the intensity in the non-quenched one. The corresponding CPTs for the
976-nm excitation are $2.3\cdot 10^{-3}$, $1.2\cdot 10^{-3}$ and $6.5\cdot
10^{-3}$ for the $hk$, $hh$, and $kk$ configurations, respectively. These
values are in good agreement with our anticipation that the phonon-assisted
reionization rate is comparable with that of the divacancy ionization to
V${}_{\text{C}}$V${}_{\text{Si}}^{-}$. Further temperature increase to 200 K
completely inhibits the quenching as seen from the corresponding spectrum for
976 nm excitation in Fig. 10. On the other hand, the quenching is well
pronounced for the 1053 nm excitation at both temperatures, 150 K and 200 K.
For this excitation (1053 nm), one notices significant change in the quenching
contrast only at 250 K, when the intensity of the “quenched” spectrum is about
30% of that of the “non-quenched” one. The corresponding theoretical CPTs are
$1.2\cdot 10^{-4}$, $4.4\cdot 10^{-4}$ and $3.1\cdot 10^{-3}$ for the $hk$,
$hh$, and $kk$ configurations, respectively. These values imply that $kk$
configuration contributes the most or entirely to the “non-quenched” spectrum,
and we still expect impact on the quenching properties for such values of the
CPT for $hk$ and $hh$ configurations. Thus, the experimental data is in
reasonable agreement with the theory and we may conclude that the dominant
reason for the loss of the quenching properties above $T=150$ K quoted in Ref.
Magnusson  _et al._ , 2018 where 976-nm excitation has been used is the
phonon-assisted reionization activated at this and higher temperatures. In
contrast, the quenching property disappears completely only at room
temperature when the 1053 nm excitation is used. Indeed, the CPT values in
these conditions are $2.4\cdot 10^{-3}$, $1.3\cdot 10^{-3}$ and $1.2\cdot
10^{-2}$ for the $hk$, $hh$, and $kk$ configurations, respectively, which are
all above the quoted CPT threshold value.

We notice the appearance of the luminescence of the silicon vacancy
(V${}_{\text{Si}}^{-}$) at elevated temperatures when the 976-nm excitation is
used, via the upconversion mechanism discussed in Ref. Magnusson  _et al._ ,
2018. Indeed, at $T=250$ K the spectrum obtained with 976-nm excitation is
dominated by the contribution of V${}_{\text{Si}}^{-}$, as shown in the bottom
left panel of Fig. 10. However, the V${}_{\text{Si}}^{-}$ spectrum can be
discerned even below $T=200$ K when 976-nm excitation is used. For instance,
the Stokes part of the V${}_{\text{Si}}^{-}$ luminescence is clearly visible
in the spectrum at $T=200$ K (left panels). In contrast, weak contribution
from the upconverted V${}_{\text{Si}}^{-}$ luminescence can be observed with
1053-nm excitation only at room temperature ($293$ K, bottom right panel); at
$T=250$ K or lower temperatures its contribution to the spectrum is negligible
and the dominant luminescence is that of the divacancy. Since the upconverted
V${}_{\text{Si}}^{-}$ luminescence is associated with increased generation of
free electrons Magnusson  _et al._ (2018), its appearance serves as an
indicator of the rapid increase of the free-electron population with
temperature when 976-nm is used, owing to the rapid increase of the CTP with
increasing temperature for this excitation. On the contrary, the increase in
the CTP (and the electron population) is much smaller with the 1053-nm
excitation, so that the contrast between quenched and non-quenched spectra
disappears only at $T=293$ K and is associated with the weak appearance of
upconverted luminescence from V${}_{\text{Si}}^{-}$, as expected.

Thus, the experimental data suggests that the phonon-assisted divacancy
reionization apparently influences the quenching behavior at temperatures and
excitations for which the CTP is of the order of low $10^{-3}$, i.e., well
below 1%. We notice that this value is much smaller than the half of the
Debye-Waller factor ($w/2$) used as a reference during the discussion of the
probabilities in the previous section. We can understand the fact that low
reionization probability (of the order of 0.1%) is sufficient to counteract
the process of accumulation of the divacancy in the negative charge state, if
we consider the model describing the quenching dynamics presented in Ref.
Magnusson  _et al._ , 2018. In this simplest model two types of centers are
considered: traps capable to emit electrons to the conduction band as a result
of photoionization (via the IR excitation) and divacancies capable of
capturing electrons. Consequently, the traps can be in two states: with
captured electron (capable to emitted an electron to the conduction band via
photoionization), and with a missing electron (capable to recapture of free
electron). The divacancies in this model are also assumed to have two charge
states: neutral (capable of capturing electron) and negatively-charged (with
captured electron). The latter state cannot be photoionized using IR
excitation below the threshold, hence the accumulation of divacancies in the
negative charge state leads to quenching of the luminescence from the neutral
charge state.

Naturally, there exist various traps capable of emitting electrons via
photoionization, for instance, the nitrogen donors, the carbon vacancy, the
silicon vacancy in double (2$-$) and triple (3$-$) negative charge states,
just to name a few. The model in Ref. Magnusson  _et al._ , 2018 assumes only
one kind of traps which may be thought as having the dominant capture cross
section for electrons and the dominant concentration. If these were the
nitrogen donors, for instance, then the capture cross-section of ionized
donors for electrons is “giant” Lax (1960) owing to the possibility for
capture via the excited donor states. On the other hand capturing of electrons
occurs to the neutral divacancy which, therefore, is assumed to have much
smaller capture cross section for electrons. Indeed, in order to reproduce
correctly the shape of the time decay of the divacancy PL the model in Ref.
Magnusson  _et al._ , 2018 admits that the ratio of capture cross sections for
electrons of the divacancy and the trap is $2\cdot 10^{-3}$. In other words,
an electron is emitted from traps and recaptured from traps many times before
it eventually can be captured by the divacancy. On the other hand, one can
assume that the probability for absorption of a photon which is stipulated by
the optical-absorption cross-section is similar for the traps and the
divacancies. Recapturing of electrons is fast to the traps, but only about 0.1
% of the total free-electron population are captured to divacancies. Thus, the
divacancy population in the negative charge state grows slowly, and this slow
process can be counteracted by phonon-assisted reionization with a CTP of
about 0.1 %.

To conclude this section, we can state that excitations with longer wavelength
(lower energy) widen the temperature range in which the quenching phenomenon
can be observed. However, in many cases when the silicon vacancy
V${}_{\text{Si}}^{-}$ is also present, excitation with too high photon
energies (e.g., close to or even above the threshold of all divacancy
configurations) leads at elevated temperatures to upconverted PL of
V${}_{\text{Si}}^{-}$ instead of excitation of the divacancy PL, as
illustrated in the low left panel of Fig. 10. Thus, if the PL of the divacancy
is to be measured at higher temperatures or at room temperature, lower
excitation energies are also preferable. It should be noted, however, that it
might be possible to stabilize the divacancy PL by Fermi-level tuning even
when excitations above the thresholds is applied, but this is a subject of a
different study and will not be discussed here.

##  IV Conclusions

In this paper, we investigated the photoluminescence lineshape of the
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ qubits in 4H SiC at $T=0$ K by DFT
calculations. We obtained the corresponding PL spectra by using the HR theory
and found that the corresponding PSB exhibits a rich set of sharp features for
each configuration. The sharp features originate from quasilocalized phonon
modes. We showed via IPR analysis that phonons are localized on the three
neighbor Si atoms around the V${}_{\text{C}}$ in the phonon energy region of
$\approx 17-35$ meV while those falling in the $\approx 35-115$ meV are
localized on the three neighbor C atoms around the V${}_{\text{Si}}$. We also
calculated the HR and DW factors of the PL spectra yielding the smallest HR
factor (and thus the largest DW factor) for
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ implying
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ to be the most attractive
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ qubit candidate in terms of extracting
the brightest optical transition at the coherent ZPL emission among the
divacancy configurations in 4H SiC. These optical properties of
V${}_{\text{C}}$V${}_{\text{Si}}^{0}(hh)$ in 4H SiC are reminescent of the
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ in the cubic polytype of SiC, i.e., 3C
SiC.

We carried out experiments and ab initio calculations on the temperature
dependent V${}_{\text{C}}$V${}_{\text{Si}}^{-}$ $\rightarrow$
V${}_{\text{C}}$V${}_{\text{Si}}^{0}$ reionization process in 4H SiC. We found
that the reionization of divacancies in 4H SiC upon unfocused green repump
laser illumination is a complex process in which very efficient traps (most
likely nitrogen donors) and the silicon vacancy defects present in the sample
are involved. For efficient driving of divacancy qubits in 4H SiC at elevated
temperatures, lower excitation energies are preferable. We found that the PL
signals of divacancies can be fully recovered at room temperature with typical
infared photoexcitation energies without applying any green repump laser for
such divacancy configurations which have stable fluorescence up to room
temperature. This result strongly indicates that the reionization threshold
energies of divacancies decrease with increasing temperature. Our calculations
revealed that not the thermal shifts of the band edges nor the charge
transition levels of divacancies caused by entropy contributions to the
formation enthalpy are responsible for an effective lowering of the
reionization threshold energies at elevated temperatures but rather the
activation of phonon-assisted photoionization transitions. This result
highlights a general conclusion that phonon-assisted photoionization processes
should be always considered in understanding the thermally activated
photoionization of point defects in semiconductors and insulators.

###### Acknowledgements.

A. Cs. acknowledges the support from the ÚNKP-20-4 New National Excellence
Program of the Ministry for Innovation and Technology from the source of the
National Research, Development and Innovation Fund. A. G. acknowledges the
National Research, Development, and Innovation Office for Hungary Grant No.
KKP129866 of the National Excellence Program of Quantum-coherent materials
project, Grant No. 127902 of the EU QuantERA Nanospin project, Grant No.
2017-1.2.1-NKP-2017-00001 of the National Quantum Technology Program, and the
Quantum Information National Laboratory supported by the Ministry of
Innovation and Technology of Hungary, as well as the EU Commission for the
H2020 Quantum technology Flagship projects ASTERIQS (Grant No. 820394). We
acknowledge the EU H2020 project QuanTELCO (Grant No. 862721). The support
provided by the Swedish Research Council (Grant No. VR 2016-05362 for I. G. I.
and No. VR 2016-04068 for N. T. S.) and the Knut and Alice Wallenberg
Foundation Grant No. KAW 2018.0071 for N. T. S. and I. G. I. is acknowledged.
We acknowledge the computational sources provided by the Swedish National
Infrastructure for Computing (SNIC) at National Computation Centre (NSC)
partially funded by the Swedish Research Council through grant agreement No.
2018-05973 and the Hungarian Governmental Information Technology Development
Agency through the project “gallium”. We acknowledge that the results of this
research have been partially achieved using the DECI resource Eagle HPC based
in Poland at Poznan with support from the PRACE aisbl and resources provided
by the Hungarian Governmental Information Technology Development Agency.

## References

  * Doherty  _et al._ (2013) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. Hollenberg, Phys. Rep. 528, 1 (2013).
  * Gali (2019) A. Gali, Nanophotonics 8, 1907 (2019).
  * Koehl  _et al._ (2011) W. F. Koehl, B. B. Buckley, F. J. Heremans, G. Calusine, and D. D. Awschalom, Nature 479, 84 (2011).
  * Castelletto  _et al._ (2013) S. Castelletto, B. C. Johnson, V. IvĂˇdy, N. Stavrias, U. T., A. Gali, and T. Oshima, Nat. Mater. 13, 151 (2013).
  * Christle  _et al._ (2015) D. J. Christle, A. L. Falk, A. P., P. V. Klimov, J. Ul Hassan, N. T. Son, J. E., O. T., and D. D. Awschalom, Nat. Mater. 14, 164 (2015).
  * Widmann  _et al._ (2015) M. Widmann, S.-Y. Lee, T. Rendler, N. T. Son, H. Fedder, S. Paik, L.-P. Yang, N. Zhao, S. Yang, I. Booker, A. Denisenko, M. Jamali, S. A. Momenzadeh, I. Gerhardt, T. Ohshima, A. Gali, E. Janzén, and J. Wrachtrup, Nature Materials 14, 164 (2015).
  * Wolfowicz  _et al._ (2020) G. Wolfowicz, C. P. Anderson, B. Diler, O. G. Poluektov, F. J. Heremans, and D. D. Awschalom, Science Advances 6, eaaz1192 (2020), publisher: American Association for the Advancement of Science Section: Research Article.
  * Lohrmann  _et al._ (2015) A. Lohrmann, N. Iwamoto, Z. Bodrog, S. Castelletto, T. Ohshima, T. J. Karle, A. Gali, P. S., J. C. McCallum, and B. C. Johnson, Nat. Comm. 6, 7783 (2015).
  * Radulaski  _et al._ (2017) M. Radulaski, M. Widmann, M. Niethammer, J. L. Zhang, S.-Y. Lee, T. Rendler, K. G. Lagoudakis, N. T. Son, E. Janzén, T. Ohshima, J. Wrachtrup, and J. Vučkoviċ, Nano Letters 17, 1782 (2017), pMID: 28225630, https://doi.org/10.1021/acs.nanolett.6b05102 .
  * Spindlberger  _et al._ (2019) L. Spindlberger, A. Csóré, G. Thiering, S. Putz, R. Karhu, J. Hassan, N. Son, T. Fromherz, A. Gali, and M. Trupke, Phys. Rev. Applied 12, 014015 (2019).
  * Nagy  _et al._ (2019) R. Nagy, M. Niethammer, M. Widmann, Y.-C. Chen, P. Udvarhelyi, C. Bonato, J. U. Hassan, R. Karhu, I. G. Ivanov, N. T. Son, J. R. Maze, T. Ohshima, O. O. Soykal, A. Gali, S.-Y. Lee, F. Kaiser, and J. Wrachtrup, Nature Communications 10, 1954 (2019).
  * Babin  _et al._ (2022) C. Babin, R. Stöhr, N. Morioka, T. Linkewitz, T. Steidl, R. Wörnle, D. Liu, E. Hesselmeier, V. Vorobyov, A. Denisenko, M. Hentschel, C. Gobert, P. Berwian, G. V. Astakhov, W. Knolle, S. Majety, P. Saha, M. Radulaski, N. T. Son, J. Ul-Hassan, F. Kaiser, and J. Wrachtrup, Nature Materials 21, 67 (2022), number: 1 Publisher: Nature Publishing Group.
  * Lee  _et al._ (2015) S.-Y. Lee, M. Niethammer, and J. Wrachtrup, Phys. Rev. B 92, 115201 (2015).
  * Kraus  _et al._ (2014) H. Kraus, V. A. Soltamov, F. Fuchs, D. Simin, A. Sperlich, P. G. Baranov, G. V. Astakhov, and V. Dyakonov, Scientific Reports 4 (2014), 10.1038/srep05303.
  * Simin  _et al._ (2016) D. Simin, V. A. Soltamov, A. V. Poshakinskiy, A. N. Anisimov, R. A. Babunts, D. O. Tolmachev, E. N. Mokhov, M. Trupke, S. A. Tarasenko, A. Sperlich, P. G. Baranov, V. Dyakonov, and G. V. Astakhov, Physical Review X 6, 031014 (2016).
  * Simin  _et al._ (2015) D. Simin, F. Fuchs, H. Kraus, A. Sperlich, P. G. Baranov, G. V. Astakhov, and V. Dyakonov, Physical Review Applied 4, 014009 (2015).
  * Niethammer  _et al._ (2016) M. Niethammer, M. Widmann, S.-Y. Lee, P. Stenberg, O. Kordina, T. Ohshima, N. T. Son, E. Janzén, and J. Wrachtrup, Physical Review Applied 6, 034001 (2016).
  * Cochrane  _et al._ (2016) C. J. Cochrane, J. Blacksberg, M. A. Anders, and P. M. Lenahan, Scientific Reports 6, 37077 (2016).
  * Anisimov  _et al._ (2016) A. N. Anisimov, D. Simin, V. A. Soltamov, S. P. Lebedev, P. G. Baranov, G. V. Astakhov, and V. Dyakonov, Scientific Reports 6, 33301 (2016).
  * Gali (2011) A. Gali, physica status solidi (b) 248, 1337 (2011).
  * Magnusson and Janzén (2005) B. Magnusson and E. Janzén, Materials Science Forum 483-485, 341 (2005), conference Name: Silicon Carbide and Related Materials 2004 ISBN: 9780878499632 Publisher: Trans Tech Publications Ltd.
  * Baranov  _et al._ (2005) P. G. Baranov, I. V. Il’in, E. N. Mokhov, M. V. Muzafarova, S. B. Orlinskii, and J. Schmidt, Jounal of Experimental and Theoretical Physics Letters 82, 441 (2005).
  * Son  _et al._ (2006) N. T. Son, P. Carlsson, J. ul Hassan, E. Janzén, T. Umeda, J. Isoya, A. Gali, M. Bockstedte, N. Morishita, T. Ohshima, and H. Itoh, Phys. Rev. Lett. 96, 055501 (2006).
  * Somogyi and Gali (2014) B. Somogyi and A. Gali, Journal of Physics: Condensed Matter 26, 143202 (2014).
  * Somogyi  _et al._ (2012) B. Somogyi, V. Zolyomi, and A. Gali, Nanoscale 4, 7720 (2012).
  * Beke  _et al._ (2020) D. Beke, J. Valenta, G. KĂˇrolyhĂˇzy, S. Lenk, Z. CzigĂˇny, B. G. MĂˇrkus, K. KamarĂˇs, F. Simon, and A. Gali, The Journal of Physical Chemistry Letters 11, 1675 (2020), publisher: American Chemical Society.
  * Christle  _et al._ (2017) D. J. Christle, P. V. Klimov, C. F. de las Casas, K. Szász, V. Ivády, V. Jokubavicius, J. Ul Hassan, M. Syväjärvi, W. F. Koehl, T. Ohshima, N. T. Son, E. Janzén, A. Gali, and D. D. Awschalom, Phys. Rev. X 7, 021046 (2017).
  * Bourassa  _et al._ (2020) A. Bourassa, C. P. Anderson, K. C. Miao, M. Onizhuk, H. Ma, A. L. Crook, H. Abe, J. Ul-Hassan, T. Ohshima, N. T. Son, G. Galli, and D. D. Awschalom, Nature Materials 19, 1319 (2020).
  * Falk  _et al._ (2014) A. L. Falk, P. V. Klimov, B. B. Buckley, V. Ivády, I. A. Abrikosov, G. Calusine, W. F. Koehl, A. Gali, and D. D. Awschalom, Physical Review Letters 112, 187601 (2014).
  * Wolfowicz  _et al._ (2021) G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom, Nature Reviews Materials 6, 906 (2021), bandiera_abtest: a Cg_type: Nature Research Journals Number: 10 Primary_atype: Reviews Publisher: Nature Publishing Group Subject_term: Electronic devices;Quantum information;Qubits Subject_term_id: electronic-devices;quantum-information;qubits.
  * Gali  _et al._ (2006) A. Gali, M. Bockstedte, N. T. Son, T. Umeda, J. Isoya, and E. Janzén, in  _Silicon Carbide and Related Materials 2005_, Materials Science Forum, Vol. 527 (Trans Tech Publications Ltd, 2006) pp. 523–526.
  * Gali (2012) A. Gali, Journal of Materials Research 27, 897â€“909 (2012).
  * Gordon  _et al._ (2015) L. Gordon, A. Janotti, and C. G. Van de Walle, Phys. Rev. B 92, 045208 (2015).
  * Magnusson  _et al._ (2018) B. Magnusson, N. T. Son, A. Csóré, A. Gällström, T. Ohshima, A. Gali, and I. G. Ivanov, Phys. Rev. B 98, 195202 (2018).
  * Anderson  _et al._ (2019) C. P. Anderson, A. Bourassa, K. C. Miao, G. Wolfowicz, P. J. Mintun, A. L. Crook, H. Abe, J. Ul Hassan, N. T. Son, T. Ohshima, and D. D. Awschalom, Science 366, 1225 (2019), publisher: American Association for the Advancement of Science.
  * Zwier  _et al._ (2015) O. Zwier, D. O’Shea, A. R. Onur, and C. H. van der Wal, Sci. Rep. 5 (2015), 10.1038/srep10931.
  * Wolfowicz  _et al._ (2017) G. Wolfowicz, C. P. Anderson, A. L. Yeats, S. J. Whiteley, J. Niklas, O. G. Poluektov, F. J. Heremans, and D. D. Awschalom, Nat. Comm. 8 (2017), 10.1038/s41467-017-01993-4.
  * Golter and Lai (2017) D. A. Golter and C. W. Lai, Sci. Rep. 7 (2017), 10.1038/s41598-017-13813-2.
  * Siyushev  _et al._ (2013) P. Siyushev, H. Pinto, M. Vörös, A. Gali, F. Jelezko, and J. Wrachtrup, Phys. Rev. Lett. 110, 167402 (2013).
  * Huang  _et al._ (1950) K. Huang, A. Rhys, and N. F. Mott, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 204, 406 (1950), https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1950.0184 .
  * Alkauskas  _et al._ (2014) A. Alkauskas, B. B. Buckley, D. D. Awschalom, and C. G. V. de Walle, New Journal of Physics 16, 073026 (2014).
  * Gali  _et al._ (2016) A. Gali, D. T., M. Vörös, G. Thiering, E. Cannuccia, and A. Marini, Nat. Comm. 7 (2016), 10.1038/ncomms11327.
  * Debye (1913) P. Debye, Annalen der Physik 348, 49 (1913).
  * Waller (1923) I. Waller, Zeitschrift fur Physik 17, 398 (1923).
  * Csóré  _et al._ (2021) A. Csóré, N. T. Son, and A. Gali, Phys. Rev. B 104, 035207 (2021).
  * Freysoldt  _et al._ (2009) C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. 102, 016402 (2009).
  * Cannuccia and Gali (2020) E. Cannuccia and A. Gali, Phys. Rev. Materials 4, 014601 (2020).
  * Wickramaratne  _et al._ (2018) D. Wickramaratne, C. E. Dreyer, B. Monserrat, J.-X. Shen, J. L. Lyons, A. Alkauskas, and C. G. Van de Walle, Applied Physics Letters 113, 192106 (2018), https://doi.org/10.1063/1.5047808 .
  * Gali  _et al._ (2009) A. Gali, E. Janzén, P. Deák, G. Kresse, and E. Kaxiras, Phys. Rev. Lett. 103, 186404 (2009).
  * Bersuker (2006) I. B. Bersuker,  _The Jahn-Teller effect_ (Cambridge University Press, 2006).
  * Heyd  _et al._ (2003) J. Heyd, G. E. Scuseria, and M. Ernzerhof, The Journal of Chemical Physics 118, 8207 (2003).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Bell  _et al._ (1970) R. J. Bell, P. Dean, and D. C. Hibbins-Butler, Journal of Physics C: Solid State Physics 3, 2111 (1970).
  * Falk  _et al._ (2013) A. L. Falk, B. B. Buckley, G. Calusine, W. F. Koehl, V. V. Dobrovitski, A. Politi, C. A. Zorman, P. X.-L. Feng, and D. D. Awschalom, Nat Commun 4, 1819 (2013).
  * Hashemi  _et al._ (2021) A. Hashemi, C. Linderälv, A. V. Krasheninnikov, T. Ala-Nissila, P. Erhart, and H.-P. Komsa, Phys. Rev. B 103, 125203 (2021).
  * Csóré  _et al._ (2022) A. Csóré, N. Mukesh, G. Károlyházy, D. Beke, and A. Gali, Journal of Applied Physics 131, 071102 (2022), publisher: American Institute of Physics.
  * Lax (1960) M. Lax, Phys. Rev. 119, 1502 (1960).