# Improving Self-supervised Learning for Out-of-distribution Task via
Auxiliary Classifier

Harshita Boonlia  CSE Department
Indian Institute of Technology Kharagpur
Kharagpur, India  Tanmoy Dam  SEIT
University of New South Wales
Canberra, Australia  Md Meftahul Ferdaus  ATMRI
Nanyang Technological University
Singapore  Sreenatha G. Anavatti  SEIT
University of New South Wales
Canberra, Australia  Hussein A. Abbass  SEIT
University of New South Wales
Canberra, Australia  Ankan Mullick  CSE Department
Indian Institute of Technology Kharagpur
Kharagpur, India

###### Abstract

In real world scenarios, the out-of-distribution (OOD) datasets may have a
large distributional shift from training datasets. This phenomena generally
occurs when a trained classifier is deployed on varying dynamic enviorments.
Hence, it causes a significance drop in performance. In this work, we are
trying to tackle this issue by end-to-end fashion. We observed that rotational
prediction (self-supervised) accuracy generally has a strong a relationship
with OOD tasks. Therefore, a multitasking learning framework is introduced
where an additional auxiliary classification head is included along with
semantic classification and rotational prediction head. The proposed learning
framework can be formulated into bi-level optimisation framework where the
upper-level is trained the parameters semantic classification and rotational
head. In the lower-level optimisation, the auxiliary classification head is
only updated through semantic classification head where semantic
classification head parameters are fixed. The proposed joint-learning
framework is maximised the rotational prediction head under given two semantic
classification head. The proposed method has been validated through three
unseen OOD datasets where it exhibits a significant improvement on performance
from other two baseline methods.

###### Index Terms:

out of distribution, self-supervised learning, auxiliary classifier

##  I Introduction

In machine learning community, benchmarks like ImageNet [deng2009imagenet],
CIFAR [krizhevsky2009learning] etc. are commonly used to know the
generalization ability of classifiers, where empirical results are computed
through held-out labeled test sets. In these benchmarks, classifiers are
evaluated through popular evaluation metrics such as top-1 error, mean average
precision, etc. However, when classifiers applied to real-world applications
like product recommendation, medical diagnosis, autonomous driving, their
performance evaluation is not straightforward. For instance, they regularly
experience novel distributions that are different from the original training
distribution. Thus, their generalization ability observed from existing
benchmarks is not reflected during deploying in real-world scenarios. To
estimate classifier’s accuracy in such scenarios, it needs to re-evaluated on
those real-world data sets. However, in many cases, annotations of test
samples are not provided. Under such environment, classifiers’ performance is
usually evaluated by collecting new labeled test sets. Nonetheless, labeling
adequate images in a novel scenario is very complex and highly expensive. To
minimize such labeling cost, researchers have investigated various approaches
for evaluating classifiers’ performance on unlabeled test sets. Some
researchers have developed complexity measurements on model parameters to
analyse generalization of the classifiers [arora2018stronger,
corneanu2020computing, jiang2018predicting]. However, these works are limited
to fixed test sets. Unlike them, varying test domains are considered in a
recent work [deng2021does], where unsupervised evaluation mechanism is
considered to calculate their classifier’s accuracy with unlabeled test sets.
To be specific, they have utilized self-supervised task to generate labels in
test sets. Thus, they can calculate the classifier’s accuracy easily without
any human annotations. Furthermore, they investigated the relationship between
a network with unlabeled test sets performing good at self-supervised task and
its successful object detection task under same settings. In their
investigation, they observed a proportional relationship between a
classifier’s performance and the self-supervised task’s accuracy. Such
relationship has been realized by measuring how semantic classification co-
relates with rotation prediction statistically under varying testing
environments. Since varying test domains are fitted with real-world settings,
we also have focused on such varying test domains in this work….

Learning effective representations for downstream tasks is an usual usage of
self-supervised tasks in machine learning community [he2020momentum,
noroozi2016unsupervised, gidaris2018unsupervised]. Besides, researchers
predicted classifier accuracy on novel test sets using the self-supervision
method [deng2021does]. They classified four different rotation angles of an
image at $\\{0^{\circ},~{}90^{\circ},~{}180^{\circ},~{}270^{\circ}\\}$ to pay
attention to the pretext task of rotation prediction. They jointly trained
their network for both classification task and pretext task using CIFAR-10,
MNIST, Tiny-ImageNet, and COCO, where they studied the correlation between
those two task’s accuracy. By considering many labeled testsets and plotting
classification verses rotation prediction accuracy, a strong correlation
(Pearson’s Correlation $r>0.88$) is witnessed between those accuracies. Based
on such finding, they learnt a linear regression model, which can predict
classification accuracy on unseen test sets. They obtained ground truths from
a given unlabeled test set by rotating images manually. It has been used to
calculate the rotation prediction accuracy on the test images using the multi-
task network. Afterward, this rotation prediction accuracy is used by the
linear regression model to predict the semantic classification accuracy.

Unlike them, in this work, we mainly focused on improving rotation prediction
accuracy (self-supervised learning) using an auxiliary classifier for out of
distribution task, which ultimately improves the semantic classification
accuracy. Therefore, our proposed approach can be formulated as a bi-level
optimization problem, where the rotation head and semantic classification head
are learnt in the upper level. On the other hand, auxiliary classification
head is learnt through the semantic classification head without updating the
semantic classifier parameters in the lower level as shown in Fig. 1. Since
better rotation prediction accuracy indicates a models higher ability to
capture OOD features, we designed our multi-task architecture by focusing to
maximize the performance of the rotation prediction head.

Major contributions in this work can be stated as follows:

  * •

We proposed a joint end-to-end multi-task framework called only self-
supervised learning (OOSL) for handling unseen OOD test sets.

  * •

We formulated the problem in a bi-level optimization fashion as such to
improve semantic classification performance by maximizing the rotation
prediction accuracy via an auxiliary classifier.

  * •

Our proposed framework has been validated using three unseen OOD data sets,
where a better semantic classification accuracy have been witnessed in
contrast with the baselines.

Figure 1: Proposed multi-task network structure for improving only the
accuracy of rotation prediction via auxiliary classifier

##  II Proposed method

In this paper, we have utilized a held-in training set and an unseen OOD test-
set. We define the given training set as
$D^{train}=\\{(\textbf{x}_{k},\textbf{y}_{k})\\}_{k=1}^{N}\in(X\times Y)$
where, $x_{k}\in X$ is the $k$-th training image and
$\textbf{y}_{k}\\{k=0,1,...,C-1\\}\in Y$ is corresponding to target variable
spanning over $C$ classes. $N$ is the total number of samples present in the
training dataset $D^{train}$. The OOD test set is defined in a similar
fashion, $D^{test}=\\{(\textbf{x}_{l},\textbf{y}_{l})\\}_{l=1}^{M}$ where
$\textbf{y}_{l}$ is target variable spanning over $C$ classes. Given a set of
observation $(\textbf{x}_{k},\textbf{y}_{k})\subseteq(X,Y)$ drawn from a joint
distribution $(\textbf{x}_{k},\textbf{y}_{k})\sim P_{XY}$, our objective is to
design a robust classifier that can maximise classification accuracy for
unseen OOD test-set.

###  II-A Multi-task learning

Researchers have observed in [gidaris2018unsupervised] that a higher accuracy
in rotation prediction indicates a model’s superiority in capturing
representation of learned features in lower dimensional manifold. Though they
followed a completely unsupervised learning strategy, they did not consider
the OOD task. In recent times, similar method has been developed to deal with
OOD problem [deng2021does], where a linearly proportional relationship between
rotation prediction and semantic classification has been observed. However,
their proposed two-stage method can not provide an end-to-end solution.
Besides, there are few questions still remains in the end-to-end method.

  * •

If rotational head estimates OOD classification accuracy perfectly, then how
do I maximise the rotational head accuracy under multi-tasking learning
framework?

  * •

How do I get an end-to-end solutions for predicting OOD downstream task?

###  II-B How to maximise the self-supervision(rotational-accuracy) task?

Network details: To attain this objective, a multi-task learning framework has
been developed for semantic classification along with self-supervised
(rotation prediction) task. We utilise the multi-head network along with the
same base network. In addition, the multitasking learning should not enforce
to more complexity while learning self-supervised task. To maximise the self-
supervised performance on base-network, we introduce an auxiliary classifier
along with semantic classifier head and self-supervised (rotation prediction)
head. That minimal changes can’t give huge burden on the base-network. For the
base(feature extraction) network, we have taken a convolutional neural network
(e.g. densenet)followed by three fully connected layers for three different
tasks. As depicted in Fig 1, the base feature extractor is parameterised by
$\theta_{b}$. All the remaining task-specific classification head are
described as follows,

  * •

semantic-classification prediction head is parameterised by $\theta_{sc}$.

  * •

rotation prediction head is parameterised by $\theta_{ss}$.

  * •

auxiliary semantic classification head is parameterised by $\theta_{sa}$.

Rotational prediction head: We follow the similar rotation transformation as
in [deng2021does, gidaris2018unsupervised]. The four geometrical rotational
transformations are applied to an train image(x), $F=\\{G_{r}(\textbf{x})\\}$,
where $G_{r}$ is the geometrical rotation function with four rotation angles
$r=\\{0^{\circ},90^{\circ},180^{\circ},270^{\circ}\\}$. This geometrical
transformation can not alter the invariant nature [dam2021improving].
Therefore, the rotational head can predict rotational accuracy by $4$-ways.

Loss functions: The proposed OSSL method is associated with three individual
classifications losses for three different tasks. The semantic classification
loss is defined as follows,

$\displaystyle L_{ch}=CE(\textbf{y}_{c},\theta_{sc}(\theta_{b}(\textbf{x})))$
(1)

where
$CE=-\frac{1}{N}\sum_{y_{c}=1}^{N}\textbf{y}_{c}\log(\theta_{sc}(\theta_{b}(\textbf{x})))$

The rotation prediction classification loss is defined as follows,

$\displaystyle
L_{rh}=\frac{1}{4}\sum_{r\in\\{0^{\circ},90^{\circ},180^{\circ},270^{\circ}\\}}CE(\textbf{y}_{r},\theta_{ss}(\theta_{b}(G_{r}(\textbf{x}))))$
(2)

where, $y_{r}$ is represented as one-hot-encode labels for all four
rotational.

The semantic auxiliary classification loss is defined as follows,

$\displaystyle L_{ah}=CE(\textbf{y}_{c},\theta_{sa}(\theta_{b}(\textbf{x})))$
(3)

We have utilized above three losses into a bi-level optimisation problem to
maximise the self-supervision performance. In upper-level optimisation, the
semantic classification head and rotation classification head parameters are
learnt simultaneously to update the base-network parameters as well as
corresponding task specific class parameters, where the objective can be
expressed as follows:

$\underset{\theta_{b},\theta_{sc},\theta_{ss}}{\text{min}}L_{upper}$ (4)

where $L_{upper}=(L_{ch}+L_{rh})$. This upper level optimization problem is
solvable with the stochastic gradient descent (SGD) method where it first
tunes the parameters of the both the task specific classifiers:

$\\{\theta_{sc},\theta_{ss},\theta_{b}\\}=\\{\theta_{sc},\theta_{ss},\theta_{b}\\}-l_{r}\sum_{D^{train}}\nabla_{\\{\theta_{sc},\theta_{ss},\theta_{b}\\}}L_{upper}$
(5)

where $l_{r}$ is the learning rate of the upper-level loop.

Similarly, in lower-level optimisation, the objective can be expressed as
follows:

$\underset{\theta_{b},\theta_{sa}}{\text{min}}L_{lower}$ (6)

where $L_{lower}=(L_{ch}-L_{rh})$. As similar to upper level, the SGD method
is used to optimize only the parameters of the auxiliary head and base
network:

$\\{\theta_{sa},\theta_{b}\\}=\\{\theta_{sa},\theta_{b}\\}-l_{r}\sum_{D^{train}}\nabla_{\\{\theta_{sa},\theta_{b}\\}}L_{lower}$
(7)

where, the same $l_{r}$ is being used to nullify the effects of semantic
classification head in the backward path. However, the semantic classification
head parameters ($\theta_{sc}$) remain fixed. The proposed OSSL learning
strategy is given in Algorithm 1.

Algorithm 1 Learning Strategy of OSSL

Input: training dataset $D^{train}$, testing dataset $D^{test}$, learning
rates $l_{r}$, iteration numbers $n_{epoch}$

Output: parameters of all the four networks
$\\{\theta_{b},\theta_{sc},\theta_{ss},\theta_{sa}\\}$,

for $p=1$ to $n_{epoch}$ do

Update $\\{\theta_{sc},\theta_{ss},\theta_{b}\\}$ parameters by using equation
(5) /* upper level optimisations/*

Update $\\{\theta_{sa},\theta_{b}\\}$ parameters by using equation (7) when
$\theta_{sc}$ is fixed /* lower level optimisations/*

if $p\geq 50\,\&\,p\geq 1$ then

calculate testing accuracy for $D^{test}$

end if

end for

##  III Experiments and Validation

In this paper, our proposed OSSL method has been compared with two other
baseline methods associated with two different losses, where the parameters of
the first baseline is updated through semantic classification loss($L_{ch}$)
and the second one is updated through semantic classification with rotational
head losses ($L_{ch}+L_{rh}$). For experimental validations, two popular
classification benchmarks data sets have been considered to train the model
such as: digits (MNIST) and natural image (CIFAR-10) data set. For both data
sets, three different unseen OOD test-sets have been utilised to evaluate the
model prediction accuracy.

The LeNet-5 [lecun2015lenet] model is a popular architecture for classifying
the digit datasets (MNIST). Therefore, we consider LeNet-5 as a base feature
extractor along with three classification heads. Original MNIST dataset is
applied to train the model parameters, but, two different unseen OOD data sets
namely USPS [hull1994database] and SVHN [netzer2011reading] are used to test
it. Besides, both the unseen test-sets are having same number of classes($10$)
as in the training set. Therefore, it is practical to use these data sets as
unseen OOD test-sets. On top of the backbone feature extractor i.e. LeNet-5,
three tasks specific fully connected layers are being used. In addition, to
analyse the effectiveness of the proposed method in a complex dataset,
$Densenet-40$ ($40$ layers) architecture [huang2017densely] is applied as a
backbone feature extractor. In this case, CIFAR-10 is used to train the model,
whereas $CIFAR-10.1$ is utilized as an unseen test-set to evaluate the model
performance. CIFAR-10.1 is a popular benchmark dataset for the OOD
classification task where collected test-set samples distributional shift
cannot vary too much compared with the original CIFAR −10 samples
[recht2018cifar10.1, torralba2008tinyimages]. Moreover, CIFAR-10.1 dataset
samples are subset of the Imagenet [krizhevsky2012imagenet] dataset samples.

TABLE I: The quantitative classification performance for unseen OOD test-set. While the base model train with held-in digit datasets, USPS and SVHN are being used for unseen(held-out) test sets. For CIFAR-10 datast, CIFAR 10.1 is used as unseen test-set. The reported classification accuracy is given in $\%$ Trainset | MNIST | CIFAR 10
---|---|---
unseen OOD test-set | USPS | SVHN | CIFAR 10.1
$L_{ch}$ | 60.85 | 22.25 | 83.30
$L_{ch}+L_{rh}$ | 64.82 | 19.74 | 85.05
OSSL | 65.22 | 21.09 | 88.00

Table I represents the quantitative performance comparison for different
unseen test-sets. For a fair comparison among these methods, the same model
architecture is considered for all cases. For digit classification problem,
original MNIST train-set is used to train the model. The obtained baseline
classification performance (considering $L_{ch}$) for USPS testset is
$60.85\%$. A significance improvement on performance has been observed when
rotation prediction loss is considered along with semantic classification
loss. However, our proposed OSSL method has obtained best classification
accuracy compared with the above two methods. The obtained classification
accuracy is $65.22\%$. However, SVHN test-set, the best obtained model
accuracy is $22.25\%$ for using only classification head. A significant
performance drop is observed while considering the rational prediction head
along with classification head. However, OSSL has obtained better model
accuracy than the rotational head prediction, but it cannot outperform the
baseline performance. The main reason could be that the unseen SVHN test set
has large distributional sift samples compared to the held-in train set
[hsu2020generalized]. Therefore, large distributional shift data can’t be
represented by geometrical rotation with given held-in train-set. As a
results, the classification head performance is declined under multi-tasking
learning framework.

In addition, we have also considered the a complex CIFAR-10 dataset where more
complex network architecture has been utilised as a baseline. It is clearly
observed from the Table I, the proposed OSSL has outperformed with other two
methods. While considering only classification head, the obtained accuracy is
$83.30\%$. A significant performance improvement has been observed by
considering both classification head and rotation head prediction accuracy.
Therefore, the obtained classification accuracy is $85.05\%$. However, the
proposed method , OSSL has attained based model accuracy among other two
methods. The obtained accuracy is $88.00\%$ and a significant improvement on
performance has been observed.

##  IV limitations

The proposed OSSL method is based on maximisation of rotational prediction
performance so that the improvement of semantic classification has been
observed. However, for some OOD datasets, rotational head can’t guarantee to
estimate unseen test-set accurately. As a consequence, rotational prediction
head has influenced to make classification head performance deteriorated and
it is due to sharing the common features by the two classification heads.
Moreover, the under-laying assumption is that rotational head prediction based
OOD task has to be well-defined and significance [sun2020test]. Otherwise,
rotational prediction head can’t capture the significance features. As it is
observed from the SVHN datasets, the training and testing dataset has large
distributional shift. However, a classifier is called robust when it can
detect and reject such large distributional shift from held-out dataset
[lee2018simple, hendrycks2016baseline, hsu2020generalized] .

##  V Conclusion & Future Works

This paper presents a joint learning strategy to improve classification for
unseen OOD downstream task. If OOD datasets has higher correlation with
rotation prediction accuracy, then we can maximise the rotation accuracy to
get better semantic classification performance. To attain this objective, we
formulated a bi-level optimisation framework where an additional auxiliary
classifier was introduced to nullify the semantic classification head impacts
on base-feature extractors. The proposed OSSL method has been validated
through three unseen OOD datasets. A significant improvement in classification
performance is observed than the two other baselines.Some of the future
directions, proposed learning strategy can be directly applied as follows,

  * •

put a set of penalty for different dynamic test enviorments through invariant
risk minimisation principle for handling large distributional shifts
[arjovsky2019invariant].

  * •

for handling adversarial robustness through maximising the rotational accuracy
is a one of the major key directions of OSSL method [chen2019improving]

  * •

For handling class imbalance in the training set, a latent preserving GAN
[dam2021does] can be generated for minority class for dynamic tests
environments.


# High temperature spectroscopy of ensembles of nitrogen vacancy centers in
diamond

Mohammed Attrash  Oleg Shtempluck  Eyal Buks<EMAIL_ADDRESS>Andrew
and Erna Viterbi Department of Electrical Engineering, Technion, Haifa 32000
Israel

###### Abstract

We study the spectroscopy of an ensemble of negatively charged nitrogen-
vacancy (NV-) centers in diamond at high temperatures between room temperature
and $700\operatorname{K}$ under high vacuum conditions. Spin resonances are
studied using optical detection of magnetic resonance (ODMR), and optical
spectroscopy is employed to study radiative transitions. Upon increasing the
temperature the intensity of radiative decay in visible and infra-red
decreased. In addition, the ODMR resonance frequencies were decreased, and the
phonon line emission shifted to higher wavelengths. Density functional theory
calculation of the zero-field splitting parameter ($D$) revealed that thermal
expansion is not enough to explain the shift in the ODMR frequencies. Fitting
the measured intensity of photo-luminescence with the theoretical predictions
of the Mott-Seitz model yields the value of $0.22\operatorname{eV}$ for the
energy barrier associated with nonradiative decay.

##  I Introduction

Negatively charged nitrogen vacancy (NV-) center in diamond is considered as a
candidate for quantum technologies due to its unique physical properties, such
as long coherence time at room temperature. The technique of optical detection
of magnetic resonance (ODMR) allows monitoring spin resonances Gruber_2012 ;
Wolf_041001 ; Juraschek_127601 ; Pinto_1 , whereas radiative transitions can
be probed using optical spectroscopy. Many researchers investigated the
temperature dependence of luminescence of NV Toyli_031001 ; Collins_2177 ;
Chen_161903 ; Doherty_041201 ; Acosta_070801 . It was found that the magnetic
resonance frequencies were decreased after increasing the temperature. The
shift of the frequency lines is attributed to the temperature dependence of
the NV- zero-field splitting parameter $D$. This dependency was explained by a
combination of two mechanisms, thermal expansion Ivady_235205 and electron-
phonon interactions Doherty_041201 . It is acclaimed that the non-radiative
processes shorten the excited-state lifetime at high temperatures
Plakhotnik_9751 ; Toyli_031001 ; Acosta_201202 . Beside shift in frequencies,
the photo-luminescence (PL) intensity was found to decrease at high
temperatures. Two models were proposed to explain PL intensity decay: Mott-
Seitz Gurney_69 ; Seitz_74 and Schon-Klasens Schon_463 ; Wise_226 ;
Klasens_306 . In the Mott-Seitz model, the PL quenching at high temperature is
due to the enhancement of the non-radiative process. On the other hand, in the
Schon-Klasens model, the PL decreases due to radiationless recombination of
holes and electrons at non-radiative recombination centers. Despite the
extensive studies of the NV at high temperatures, some remaining issues were
not thoroughly investigated.

Figure 1: Experimental setup. The diamond wafer is attached to a heater stage
inside the HV chamber. Both a MW LA and a multi-mode optical fiber are
attached to a copper-made finger cooled by liquid nitrogen. Vacuum
feedthroughs are employed to connect both MW LA and optical fiber to room
temperature instruments. The bare end of the optical fiber is positioned in
front of the diamond wafer. The fiber is split to 3 components (using fiber
couplers) , which are connected to the PD, VS and IRS.  The NV- is excited by
a $532\operatorname{nm}$ laser delivered via free space. The optical
propagation direction is parallel to the [110] plane. The magnetic field is
applied by a cylindrical neodymium magnet. Figure 2: (a) Schematics of the NV-
electronic structure, the ${}^{3}\mathrm{A}_{2}$ and ${}^{3}\mathrm{E}$ are
the triplet ground and excited states, respectively, while 1A1 and 1E are
singlet states, which involved in the intersystem crossing. The solid lines
represent the optical transitions, while the dashed lines represent the non-
radiative transitions. (b) NV unit cell structure, where carbon atoms are in
gray, nitrogen in blue, and vacancy in white. (c) ODMR plot at room
temperature. The color bar indicates the normalized ODMR intensity. The dashed
lines represent the resonances related to spin transitions from
$m_{\mathrm{s}}=\pm 1$ to $m_{\mathrm{s}}=0$ for the four different NV- center
directions (the resonances corresponding to two directions are nearly
degenerate). (d) PL measurement of HPHT diamond [110] before and after
electron irradiation and annealing. For HPHT diamond, two features are
observed: at 573 nm, which is related to diamond Raman peak Mildren_18950 ;
Chandran_221602 , and at 637 nm, which is related to NV- zero-phonon line.
After electron irradiation and annealing, a very low peak intensity at 575 nm
can be observed which is related to NV0 centers in diamond. In addition, the
phonon line intensity (600 - 850 nm) increased.

In this work we employ both ODMR in the microwave (MW) band, and optical
spectroscopy in the optical band of $500-1500\operatorname{nm}$, to study
ensembles of NV- centers in a high vacuum (HV) system at different
temperatures ($300-550\operatorname{K}$). In addition, we investigate the PL
as a function of temperature and calculate the energy barrier according to
Mott-Seitz model of ensembles of NV- centers.  The comparison between our
experimental results and theory yields partial agreement.

##  II Experimental setup

Type Ib high pressure high temperature (HPHT) single crystal diamond with a
nitrogen concentration lower than $200$ ppm was laser-cut along the [110]
plane, polished and irradiated with $2.8\operatorname{MeV}$ electrons at a
doze of $8\times 10^{18}\operatorname{cm}^{-2}$, and annealed at
$900^{\circ}$C for 2 hours. The estimated NV- concentration is $3.3\times
10^{17}\operatorname{cm}^{-3}$ Alfasi_214111 . The diamond was cleaned in an
acid mixture of Perchloric, Sulfuric, and Fuming Nitric acid for 1 hour.

The experimental setup is schematically shown in Fig. 1. The vacuum system,
which has a base pressure of $1\times 10^{-6}\operatorname{Torr}$, includes an
ultra HV heater (HeatWave Labs Inc., model 101491) connected to a thermocouple
(type K) and a temperature controller (model 101303), MW loop antenna (LA),
and a bare multi-mode optical fiber pointing to the diamond to collect PL
emitted from the diamond. Liquid nitrogen was flowed in a metal finger
attached to the LA and the optical fiber to avoid annealing damage. The
optical fiber is split (1x2 Step-Index Multimode Fiber Optic Couplers -
Thorlabs) to 3 components, which are connected to: photo-diode (PD) (Thorlabs,
PDA100A), visible spectrometer (VS) (Thorlabs, CCS175,
$500-1000\operatorname{nm}$), and infrared spectrometer (IRS) (Ibsen
photonics, ROCK NIR $900-1700\operatorname{nm}$). The relative optical power
delivered to the PD, VS and IRS is $0.5$, $0.25$ and $0.25$, respectively. The
spectrometers allow monitoring both the triplet-triplet transition at
wavelength of $637\operatorname{nm}$ and the singlet-singlet transition at
$1042\operatorname{nm}$ [see the schematic energy level diagram shown in Fig.
2(a)].

The NV- centers were excited by a $532\operatorname{nm}$ green laser delivered
via free space. The LA is connected to a radio frequency signal source (RFSG),
and the MW amplitude was modulated by a $151\operatorname{Hz}$ sine wave. The
PD is attached to a $600\operatorname{nm}$ long-pass filter, and its signal
was demodulated by a lock-in amplifier. To apply magnetic field on the
diamond, a cylindrical neodymium magnet is positioned outside the vacuum
chamber using a motorized stage. To reduce laser intensity fluctuations, we
used proportional integral derivative (PID) controller (SIM 960 - Stanford
Research Systems). A beam splitter in front of the laser together with a PD
are used to generate the PID input signal. The longtime optical intensity
stability is improved by more that a factor of 10 (compared with the stability
obtained without the PID).

##  III ODMR

Ignoring the hyper-fine interaction, the NV- spin triplet ground state
Hamiltonian $\mathcal{H}_{\mathrm{NV}}$ is given by

$\frac{\mathcal{H}_{\mathrm{NV}}}{\hbar}=\frac{D\mathcal{S}_{z}^{2}}{\hbar^{2}}-\frac{\gamma_{\mathrm{e}}\vec{B}\cdot\mathcal{S}}{\hbar}+\frac{E_{\mathrm{NV}}\left(S_{+}^{2}+S_{-}^{2}\right)}{2\hbar^{2}}\;,$
(1)

where $D$ and $E_{\mathrm{NV}}$ are the zero field splitting parameters, which
equale to $2\pi\times 2.88\operatorname{GHz}$ and $2\pi\times
10\operatorname{MHz}$ at $3.6\operatorname{K}$, respectively Alfasi_214111 ,
$\mathcal{S}$ is the total spin operator
$\mathcal{S}=\mathcal{S}_{x}\hat{x}+\mathcal{S}_{y}\hat{y}+\mathcal{S}_{z}\hat{z}$
, $\gamma_{\mathrm{e}}$ $=2\pi\times
28.03\operatorname{GHz}\operatorname{T}^{-1}$ is the electron spin
gyromagnetic ratio, $\vec{B}$ is the magnetic field, and
$\mathcal{S}_{\pm}=\mathcal{S}_{x}\pm i\mathcal{S}_{y}$. In this Hamiltonian
we neglect also the on-axis electric field and some strain components which
contribute to the $\mathcal{S}_{z}^{2}$ term and the off-axis electric field
and other strain components which contribute to $\mathcal{S}_{x}$ and
$\mathcal{S}_{y}$ terms Udvarhelyi_075201 . Under continuous laser excitation,
the NV- is polarized to the spin state $m_{\mathrm{s}}=0$, which has a
brighter PL Doherty_1 .

Figure 2(a) shows a schematic energy level diagram of an NV- center. The
structure of nitrogen vacancy center in diamond is shown in Fig. 2(b).
Measured ODMR of our sample as a function of magnetic field is shown in Fig.
2(c). The eight lines that are observed in this figure represent the
resonances related to spin transitions from $m_{\mathrm{s}}=\pm 1$ to
$m_{\mathrm{s}}=0$ for the four different NV- center directions. The dashed
lines represent the theoretically calculated frequencies of the same
resonances, which are obtained by numerically diagonlizing the Hamiltonian
$\mathcal{H}_{\mathrm{NV}}$ (1). Note that two pairs of resonances, which are
represented by the pink dash lines, are nearly degenerate (note that their
ODMR signal is nearly double the signal of the other NV directions). Figure
2(d) reveals the PL measurement of the HPHT diamond sample before and after
electron irradiation and annealing. The measurement was recorded by Micro-
Raman spectrometer (LABRAM HR, HORIBA, Jobin Yvon) using
$532\operatorname{nm}$ wavelength laser at room temperature (laser power was
lower than 10 mW). The HPHT diamond reveals a peak at $573\operatorname{nm}$
related to diamond Raman peak Mildren_18950 ; Chandran_221602 , and another
low intensity peak at $637\operatorname{nm}$ that is related to the zero
phonon line of NV- [see Fig. 2(a)] Prawer_2537 ; Rondin_115449 . After
electron irradiation and annealing process, the intensity of NV- zero phonon
line increased, beside the PL intensity increase in the entire range of
$650-850\operatorname{nm}$. The intensity of the NV0 zero phonon line is very
low in the treated diamond sample. The reason for relatively low NV0 peak
intensity may be related to low laser power, as was found in Ref. Manson_1705
.

##  IV Zero field splitting $D$

Figure 3: (a) and (b) ODMR as a function of temperature recorded at magnetic
field of $0$ and $10\operatorname{mT}$, respectively. The ODMR frequencies
decrease as temperature increases due to the decrease of $D$. Fitting of $D$
was done according to the parameters presented in Ref. Toyli_031001 . (c) The
shift $\Delta D$ in $D$ as a function of temperature. Note that $D$ is assumed
to be independent on the externally applied magnetic field for all fitting
lines in (a) and (b). The previous experimental data [see blue stars in (c)]
are from Feshchenko_1609_06519 ; Doherty_041201 .

To numerically investigate the shift in the zero field splitting $D$, density
functional theory (DFT) Hohenberg_B864 ; Kohn_A1133 was used as implemented in
the QUANTUM ESPRESSO software Giannozzi_395502 . The defect structure was
modeled with $3\times 3\times 2$ and $2\times 2\times 2$ super-cell of diamond
(144 and 64 atoms, respectively) and $\Gamma$ point sampling of Brillouin
zone. The calculation were performed with projector augmented-wave (PAW)
method Blochl_17953 and Perdew-Burke-Ernzerhof (PBE) exchange-correlation
functional Perdew_3865 with a plane-wave cutoff of $300$ Ry. To investigate
the thermal expansion effect on $D$, the lattice length of the structure was
chosen to suit the temperature as found in previous experiments Skinner_39 ;
Sato_092102 , and the structure was relaxed at constant lattice length until
the total force on the atoms was lower than $26\times
10^{-3}\operatorname{eV}\operatorname{\text{\AA}}^{-1}$. The zero field
splitting $D$ is believed to be related to spin-spin interactions of unpaired
electrons Lenef_13441 . Its tensor components $D_{ab}$ are given by
harriman2013theoretical

$\displaystyle D_{ab}$
$\displaystyle=\frac{\left(\gamma_{\mathrm{e}}\hbar\right)^{2}}{2S(2S-1)}\frac{\mu_{0}}{4\pi}$
$\displaystyle\times\sum_{i<j}\chi_{ij}\left\langle\psi_{ij}\right|\frac{r^{2}\delta_{ab}-3r_{a}r_{b}}{r^{5}}\left|\psi_{ij}\right\rangle\;,$

where $a,b\in\left\\{x,y,z\right\\}$ are Cartesian indices, $\mu_{0}$ is the
magnetic permeability of free space, the summation includes all pairs of
occupied Kohn-Sham orbitals, $\chi_{ij}=\pm 1$ for parallel and anti-parallel
spins, respectively, $S$ is the spin number, and $\psi_{ij}(r,r^{\prime})$ is
$2\times 2$ determinant which can be written as
$\psi_{ij}(r,r^{\prime})=2^{-1/2}[\psi_{i}(r)\psi_{j}(r^{\prime})-\psi_{i}(r^{\prime})\psi_{j}(r)]$.
A Python package named PyZFS was used to calculate $D_{ab}$ Ma_2160 .

According to our DFT calculations, the $D$ at zero temperature equals to
$2.84\operatorname{GHz}$ for $3\times 3\times 2$ super-cell
($2.77\operatorname{GHz}$ for $2\times 2\times 2$ super-cell) while the
experimental value at $3.6\operatorname{K}$ equals to
$2.88\operatorname{GHz}$. For comparison to other DFT results, $D$ was found
in the range of $\left[2.7,3\right]\operatorname{GHz}$ Bhandari_014115 ;
Whiteley_490 ; Ghosh_043801 . The shift of $D$ at a given temperature,
relative to its value at zero temperature, is denoted by $\Delta D$. Some
studies calculated the shift $\Delta D$ according to the thermal expansion at
temperatures below $300\operatorname{K}$ Doherty_041201 ; Ivady_235205 . To
account for our experimental results, we calculate $\Delta D$ based on the
thermal expansion in the range from absolute zero to $500\operatorname{K}$.

Figure 4: (a) Optical emission spectra of NV- at different temperatures. The
peak intensity shifted to higher wavelength at increased temperatures. (b)
Peak shift $\Delta E_{\mathrm{V}}$ as a function of the temperature $T$.

In order to investigate the temperature effect on the ODMR, the diamond was
measured at different temperatures under constant magnetic field. The ODMR at
temperatures below $520\operatorname{K}$ [see Fig 3(a)] reveals resonance
lines related to NV-. Upon increasing the temperature, the resonance
frequencies and their intensity were decreased, and above
$520\operatorname{K}$ the resonance lines could not be resolved. It is
expected that the shift in resonance lines is due to the decrease of $D$. The
measured shift $\Delta D$ in $D$ [see Fig. 3(b)] at $500\operatorname{K}$
equals to $-32\operatorname{MHz}$ and the slope $\mathrm{d}D/\mathrm{d}T$ is
of order of $-100\operatorname{kHz}\operatorname{K}^{-1}$, in agreement with
the value reported in Toyli_031001 . The normalized slope
$\mathrm{d}D/(D\mathrm{d}T)$ is $-4.3\times 10^{-5}\operatorname{K}^{-1}$,
which is higher than what was found in experiments carried out at lower
temperatures, for which the value of $-2.61\times
10^{-5}\operatorname{K}^{-1}$ has been measured Acosta_070801 ; Chen_161903 .

The experimental results of $D$ at temperatures lower than
$300\operatorname{K}$ were taken from other reports Chen_161903 ;
Doherty_041201 . As can be observed from simulation results, $D$ decreases as
the temperature increases. The shift in $D$ at $500\operatorname{K}$ compare
to its value at $300\operatorname{K}$ according to our calculations is
$-10\operatorname{MHz}$ ($-4\operatorname{MHz}$ for the $2\times 2\times 2$
super-cell). As can be seen from the plot in Fig. 3(b), the thermal expansion
is not enough to explain the shift in $D$. The unit cell length $a$ of diamond
at $298\operatorname{K}$ and $506\operatorname{K}$ equals to
$3.5668\operatorname{\text{\AA}}$ and $3.5680\operatorname{\text{\AA}}$,
respectively Skinner_39 , therefore,  $\mathrm{d}a/(a\mathrm{d}T)\simeq
1.6\times 10^{-6}\operatorname{K}^{-1}$ in this range. This value is too small
to account for the above-mentioned experimental finding that
$\mathrm{d}D/(D\mathrm{d}T)$ is $-4.3\times 10^{-5}\operatorname{K}^{-1}$ in
the same range.

According to a first-principles calculation reported in Tang_2205_02791 , the
shift in $D$ at $500\operatorname{K}$ is $-18\operatorname{MHz}$, while the
shift in our data is $-32\operatorname{MHz}$. Good agreement is found between
our experimental results, and the ones reported in Ref. Toyli_031001 , in
which the zero-field splitting parameter $D$ is expressed as a power series of
the temperature $T$.

Figure 5: IR emission as a function of temperature $T$ and wavelength
$\lambda$. The peak at 1042 nm, which represent the radiative decay in the
intermediate singlet states 1A${}_{1}\xrightarrow{}^{1}$E, can be observed at
temperatures below $450\operatorname{K}$. Figure 6: PL intensity
$I_{\mathrm{PL}}$ as a function of temperature $T$. The red solid line is
calculated using Eq. (4).

The ODMR signal intensity can be derived from the rate equation that governs
the time evolution of the spin polarization (see appendix A of Ref.
Alfasi_214111 ). As can be seen from Fig. 3, ODMR signals could not be
experimentally resolved above a temperature of about $520\operatorname{K}$. On
the other hand, PL signals were measured at higher temperatures up to about
$700\operatorname{K}$ (see Fig. 6). This observation suggests that the
dominant mechanism responsible for the observed ODMR signal drop with
temperature [see Fig. 3(a)], is dephasing, which enhances spin transverse
relaxation, and consequently disables ODMR above $520\operatorname{K}$.

##  V Optical spectroscopy

The fluorescence emission spectra of NV- (measured by the VS) at different
temperatures is shown in Fig. 4. The fluorescence can be resolved at
temperatures up to $520\operatorname{K}$. The spectra is shifted to longer
wavelengths at higher temperatures. The shift in wavelength at which the
emission intensity is maximized, which is denoted by $\lambda_{\mathrm{max}}$,
equals to $hc/60\operatorname{meV}$ at $500\operatorname{K}$, where $h$ is the
Planck constant and $c$ is the speed of light, and thus
$\mathrm{d}\lambda_{\mathrm{max}}/(\lambda_{\mathrm{max}}\mathrm{d}T)\simeq
3.1\times 10^{-4}\operatorname{K}^{-1}$. Fitting of the phonon shift was
performed according to Liu_3391

$\Delta
E_{v}(T)=\frac{-2A}{\exp\left({\frac{\hbar\omega_{0}}{2k_{\mathrm{B}}T}}\right)-1}\;,$
(3)

where $A$ depends on the diamond dispersion curves details, $k_{\mathrm{B}}$
is the Boltzmann’s constant, and $\omega_{0}$ is the zone-center phonon
frequency. The fitting yields that $\hbar\omega_{0}=173\pm
10\operatorname{meV}$ and $A=246\pm 34\operatorname{meV}$. Note that the zone-
center optical phonon Raman energy of diamond equals to
$165\operatorname{meV}$ ($1332\operatorname{cm^{-1}}$) Solin_1687 .

The IR spectra of the NV- ensemble (measured by the IRS) at different
temperatures is shown in Fig. 5. The peak at $1042\operatorname{nm}$, which is
associated with radiation decay from the ${}^{1}\mathrm{A}_{1}$ to
$\mathrm{E}_{1}$ levels [see Fig.2(a)], is observed at room temperatures. Upon
heating the diamond, the $1042\operatorname{nm}$ peak intensity decreased and
became unresolved at temperatures above $450\operatorname{K}$. Temperature
dependence of the ${}^{1}\mathrm{A}_{1}$ to $\mathrm{E}_{1}$ transition
wavelength has been observed in Doherty_041201 . However, this dependency
cannot be resolved in our measurements due to limited IRS resolution. Another
feature can be observed at $1350\operatorname{nm}$. We find that this feature
disappears when a 808 nm long pass filter is added, i.e. it originates from
second order diffraction (note that $1350/2=675$). The singlet NV- transition
(from ${}^{1}\mathrm{A}_{1}$ to $\mathrm{E}_{1}$ states) wavelength is found
to be less sensitive to temperature compare to the triplet transition (from
${}^{3}\mathrm{A}_{2}$ to ${}^{3}\mathrm{E}$ states) wavelength.

Finally, we discuss the PL intensity and energy barrier for the nonradiative
process, which is denoted by $U_{\mathrm{b}}$. The PL intensity
$I_{\mathrm{PL}}$ as a function of the temperature $T$ is shown in Fig. 6. In
the Mott-Seitz model this dependency of $I_{\mathrm{PL}}$ on $T$ is given by
[see Eq. (2) of Ref. Akselrod_3364 , and Eq. (5) of Ref. Chithambo_1880 ]

$I_{\mathrm{PL}}=\frac{I_{0}}{1+C\exp\left(-\frac{U_{\mathrm{b}}}{k_{\mathrm{B}}T}\right)}\;,$
(4)

where $I_{0}$ is the zero temperature intensity, and $C$ is a constant.
Fitting the measured PL intensity $I_{\mathrm{PL}}$ with Eq. (4) yields
$U_{\mathrm{b}}$ $=0.22\pm 0.05\operatorname{eV}$ and $C=200$. Note that
$U_{\mathrm{b}}$ was found in ref. Toyli_031001 to be $0.48\operatorname{eV}$.
The deviation between this result Toyli_031001 , which was obtained with a
single NV- center, and our extracted value for the energetic barrier, can
perhaps be attributed to impurities near the single NV center that was studied
in Ref. Toyli_031001 . The effect of nearby impurities on coherence lifetime
of NV- centers has been studied in Maze_644 .

##  VI Summary

In this study we investigated the optical properties of an ensemble of NV-
centers in diamond as a function of temperature. The NV- emission in the
visible ($\sim 700\operatorname{nm}$) and IR ($\sim 1042\operatorname{nm}$)
regions decreased upon increasing the temperature, and the phonon line
wavelength at maximum intensity of the NV- ($650-800\operatorname{nm}$)
increased. The $D$ parameter decreased upon increasing the temperature, and
the thermal expansion model is not enough to explain this behaviour. By
applying the Mott-Seitz model we find that the energy barrier for nonradiative
decay $U_{\mathrm{b}}$ is in the range
$\left[0.17,0.27\right]\operatorname{eV}$.

This work was supported by the Israeli science foundation, the Israeli
ministry of science, and by the Technion security research foundation. We
would like to thank Joan Adler for providing computational facilities.

## References

  * (1) A Gruber, A Dräbenstedt, C Tietz, L Fleury, J Wrachtrup, and C Von Borczyskowski,  “Scanning confocal optical microscopy and magnetic resonance on single defect centers”,  Science, vol. 276, no. 5321, pp. 2012–2014, 1997.
  * (2) Thomas Wolf, Philipp Neumann, Kazuo Nakamura, Hitoshi Sumiya, Takeshi Ohshima, Junichi Isoya, and Jörg Wrachtrup,  “Subpicotesla diamond magnetometry”,  Physical Review X, vol. 5, no. 4, pp. 041001, 2015.
  * (3) Dominik M. Juraschek, Quintin N. Meier, Morgan Trassin, Susan E. Trolier-McKinstry, Christian L. Degen, and Nicola A. Spaldin,  “Dynamical magnetic field accompanying the motion of ferroelectric domain walls”,  Phys. Rev. Lett., vol. 123, pp. 127601, Sep 2019.
  * (4) Dinesh Pinto, Domenico Paone, Bastian Kern, Tim Dierker, René Wieczorek, Aparajita Singha, Durga Dasari, Amit Finkler, Wolfgang Harneit, Jörg Wrachtrup, et al.,  “Readout and control of an endofullerene electronic spin”,  Nature communications, vol. 11, no. 1, pp. 1–6, 2020.
  * (5) D. M. Toyli, D. J. Christle, A. Alkauskas, B. B. Buckley, C. G. Van de Walle, and D. D. Awschalom,  “Measurement and control of single nitrogen-vacancy center spins above 600 k”,  Phys. Rev. X, vol. 2, pp. 031001, Jul 2012.
  * (6) A. T. Collins, M. F. Thomaz, and M. I. B. Jorge,  “Luminescence decay time of the 1.945 eV centre in type Ib diamond”,  Journal of Physics C Solid State Physics, vol. 16, no. 11, pp. 2177–2181, April 1983.
  * (7) X.-D. Chen, C.-H. Dong, F.-W. Sun, C.-L. Zou, J.-M. Cui, Z.-F. Han, and G.-C. Guo,  “Temperature dependent energy level shifts of nitrogen-vacancy centers in diamond”,  Applied Physics Letters, vol. 99, no. 16, pp. 161903, 2011.
  * (8) M. W. Doherty, V. M. Acosta, A. Jarmola, M. S. J. Barson, N. B. Manson, D. Budker, and L. C. L. Hollenberg,  “Temperature shifts of the resonances of the ${\mathrm{nv}}^{-}$ center in diamond”,  Phys. Rev. B, vol. 90, pp. 041201, Jul 2014.
  * (9) V. M. Acosta, E. Bauch, M. P. Ledbetter, A. Waxman, L.-S. Bouchard, and D. Budker,  “Temperature dependence of the nitrogen-vacancy magnetic resonance in diamond”,  Phys. Rev. Lett., vol. 104, pp. 070801, Feb 2010.
  * (10) Viktor Ivády, Tamás Simon, Jeronimo R. Maze, I. A. Abrikosov, and Adam Gali,  “Pressure and temperature dependence of the zero-field splitting in the ground state of nv centers in diamond: A first-principles study”,  Phys. Rev. B, vol. 90, pp. 235205, Dec 2014.
  * (11) Taras Plakhotnik and Daniel Gruber,  “Luminescence of nitrogen-vacancy centers in nanodiamonds at temperatures between 300 and 700 k: perspectives on nanothermometry”,  Phys. Chem. Chem. Phys., vol. 12, pp. 9751–9756, 2010.
  * (12) V. M. Acosta, A. Jarmola, E. Bauch, and D. Budker,  “Optical properties of the nitrogen-vacancy singlet levels in diamond”,  Phys. Rev. B, vol. 82, pp. 201202, Nov 2010.
  * (13) R. W. Gurney and N. F. Mott,  “Luminescence in solids”,  Trans. Faraday Soc., vol. 35, pp. 69–73, 1939.
  * (14) Frederick Seitz,  “An interpretation of crystal luminescence”,  Trans. Faraday Soc., vol. 35, pp. 74–85, 1939.
  * (15) Michael Schön,  “Zum leuchtmechanismus der kristallphosphore”,  Zeitschrift für Physik, vol. 119, no. 7, pp. 463–471, 1942\.
  * (16) M. E. Wise and H. A. Klasens,  “Fluorescence efficiency and hole migration in zinc sulphides”,  J. Opt. Soc. Am., vol. 38, no. 3, pp. 226–231, Mar 1948.
  * (17) H. Klasens,  “Transfer of energy between centres in zinc sulphide phosphors”,  Nature, vol. 158, pp. 306–307, 1946.
  * (18) Richard P Mildren, James E Butler, and James R Rabeau,  “Cvd-diamond external cavity raman laser at 573 nm”,  Optics express, vol. 16, no. 23, pp. 18950–18955, 2008.
  * (19) Maneesh Chandran, Shaul Michaelson, Cecile Saguy, and Alon Hoffman,  “Fabrication of a nanometer thick nitrogen delta doped layer at the sub-surface region of (100) diamond”,  Applied Physics Letters, vol. 109, no. 22, pp. 221602, 2016.
  * (20) Nir Alfasi, Sergei Masis, Oleg Shtempluck, and Eyal Buks,  “Detection of paramagnetic defects in diamond using off-resonance excitation of nv centers”,  Phys. Rev. B, vol. 99, pp. 214111, Jun 2019.
  * (21) Péter Udvarhelyi, V. O. Shkolnikov, Adam Gali, Guido Burkard, and András Pályi,  “Spin-strain interaction in nitrogen-vacancy centers in diamond”,  Phys. Rev. B, vol. 98, pp. 075201, Aug 2018.
  * (22) Marcus W Doherty, Neil B Manson, Paul Delaney, Fedor Jelezko, Jörg Wrachtrup, and Lloyd CL Hollenberg,  “The nitrogen-vacancy colour centre in diamond”,  Physics Reports, vol. 528, no. 1, pp. 1–45, 2013.
  * (23) Steven Prawer and Robert J Nemanich,  “Raman spectroscopy of diamond and doped diamond”,  Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, vol. 362, no. 1824, pp. 2537–2565, 2004.
  * (24) L. Rondin, G. Dantelle, A. Slablab, F. Grosshans, F. Treussart, P. Bergonzo, S. Perruchas, T. Gacoin, M. Chaigneau, H.-C. Chang, V. Jacques, and J.-F. Roch,  “Surface-induced charge state conversion of nitrogen-vacancy defects in nanodiamonds”,  Phys. Rev. B, vol. 82, pp. 115449, Sep 2010.
  * (25) NB Manson and JP Harrison,  “Photo-ionization of the nitrogen-vacancy center in diamond”,  Diamond and related materials, vol. 14, no. 10, pp. 1705–1710, 2005\.
  * (26) AV Feshchenko, O-P Saira, and JP Pekola,  “Thermal conductance of nb thin films at sub-kelvin temperatures”,  arXiv:1609.06519, 2016.
  * (27) P. Hohenberg and W. Kohn,  “Inhomogeneous electron gas”,  Phys. Rev., vol. 136, pp. B864–B871, Nov 1964.
  * (28) W. Kohn and L. J. Sham,  “Self-consistent equations including exchange and correlation effects”,  Phys. Rev., vol. 140, pp. A1133–A1138, Nov 1965.
  * (29) Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, Andrea Dal Corso, Stefano de Gironcoli, Stefano Fabris, Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, Anton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitsonen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch,  “QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials”,  Journal of Physics: Condensed Matter, vol. 21, no. 39, pp. 395502, sep 2009.
  * (30) P. E. Blöchl,  “Projector augmented-wave method”,  Phys. Rev. B, vol. 50, pp. 17953–17979, Dec 1994.
  * (31) John P. Perdew, Kieron Burke, and Matthias Ernzerhof,  “Generalized gradient approximation made simple”,  Phys. Rev. Lett., vol. 77, pp. 3865–3868, Oct 1996.
  * (32) Brian J. Skinner,  “The Thermal Expansions of Thoria, Periclase and Diamond*”,  American Mineralogist, vol. 42, no. 1-2, pp. 39–55, 02 1957.
  * (33) Toshimaro Sato, Kazutoshi Ohashi, Tomoko Sudoh, Katsuji Haruna, and Hiroshi Maeta,  “Thermal expansion of a high purity synthetic diamond single crystal at low temperatures”,  Phys. Rev. B, vol. 65, pp. 092102, Feb 2002.
  * (34) A. Lenef and S. C. Rand,  “Electronic structure of the n-v center in diamond: Theory”,  Phys. Rev. B, vol. 53, pp. 13441–13455, May 1996.
  * (35) John E Harriman,  Theoretical foundations of electron spin resonance: physical chemistry: a series of monographs,  Academic press, 1978.
  * (36) He Ma, Marco Govoni, and Giulia Galli,  “Pyzfs: A python package for first-principles calculations of zero-field splitting tensors”,  Journal of Open Source Software, vol. 5, no. 47, pp. 2160, 2020\.
  * (37) Churna Bhandari, Aleksander L. Wysocki, Sophia E. Economou, Pratibha Dev, and Kyungwha Park,  “Multiconfigurational study of the negatively charged nitrogen-vacancy center in diamond”,  Phys. Rev. B, vol. 103, pp. 014115, Jan 2021.
  * (38) Samuel J Whiteley, Gary Wolfowicz, Christopher P Anderson, Alexandre Bourassa, He Ma, Meng Ye, Gerwin Koolstra, Kevin J Satzinger, Martin V Holt, F Joseph Heremans, et al.,  “Spin–phonon interactions in silicon carbide addressed by gaussian acoustics”,  Nature Physics, vol. 15, no. 5, pp. 490–495, 2019.
  * (39) Krishnendu Ghosh, He Ma, Vikram Gavini, and Giulia Galli,  “All-electron density functional calculations for electron and nuclear spin interactions in molecules and solids”,  Phys. Rev. Materials, vol. 3, pp. 043801, Apr 2019.
  * (40) Hao Tang, Ariel Rebekah Barr, Guoqing Wang, Paola Cappellaro, and Ju Li,  “First-principles calculation of the temperature-dependent transition energies in spin defects”,  2022\.
  * (41) Ming S. Liu, Les A. Bursill, S. Prawer, and R. Beserman,  “Temperature dependence of the first-order raman phonon line of diamond”,  Phys. Rev. B, vol. 61, pp. 3391–3395, Feb 2000.
  * (42) S. A. Solin and A. K. Ramdas,  “Raman spectrum of diamond”,  Phys. Rev. B, vol. 1, pp. 1687–1698, Feb 1970.
  * (43) M. S. Akselrod, N. Agersnap Larsen, V. Whitley, and S. W. S. McKeever,  “Thermal quenching of f-center luminescence in al2o3:c”,  Journal of Applied Physics, vol. 84, no. 6, pp. 3364–3373, 1998\.
  * (44) M L Chithambo,  “The analysis of time-resolved optically stimulated luminescence: II. computer simulations and experimental results”,  Journal of Physics D: Applied Physics, vol. 40, no. 7, pp. 1880–1889, mar 2007.
  * (45) JR Maze, PL Stanwix, JS Hodges, S Hong, JM Taylor, P Cappellaro, L Jiang, MV Gurudev Dutt, E Togan, AS Zibrov, et al.,  “Nanoscale magnetic sensing with an individual electronic spin in diamond”,  Nature, vol. 455, no. 7213, pp. 644–647, 2008.