# Knesset-DictaBERT:
A Hebrew Language Model for Parliamentary Proceedings

Gili Goldin
University of Haifa
<EMAIL_ADDRESS>
Shuly Wintner
University of Haifa
<EMAIL_ADDRESS>

###### Abstract

We present Knesset-DictaBERT, a large Hebrew language model fine-tuned on the
Knesset Corpus, which comprises Israeli parliamentary proceedings. The model
is based on the DictaBERT architecture and demonstrates significant
improvements in understanding parliamentary language according to the MLM
task. We provide a detailed evaluation of the model’s performance, showing
improvements in perplexity and accuracy over the baseline DictaBERT model.

Knesset-DictaBERT:
A Hebrew Language Model for Parliamentary Proceedings

Gili Goldin University of Haifa<EMAIL_ADDRESS>Shuly Wintner University
of Haifa<EMAIL_ADDRESS>

##  1 Introduction

The field of natural language processing (NLP) has seen remarkable
advancements in recent years, driven by the development of large language
models. These models have significantly enhanced the ability to understand and
generate human language. However, much of the effort of creating and training
NLP models has been focused on English, while fewer NLP models are available
for lower-resource languages such as Hebrew. We present one such model here.

Parliamentary proceedings are a valuable source of information for
understanding political discourse, legislative processes, and more. Analyzing
these texts requires models that can accurately capture the nuances of the
language used in such settings. Despite the importance of this task, there has
been a lack of specialized models trained on parliamentary corpora in Hebrew.
To address this gap, we fine-tuned the pre-trained Dicta-BERT model (Shmidman
et al., 2023) on the Knesset Corpus (Goldin et al., 2024), a dataset of
Israeli parliamentary proceedings. DictaBERT is a state-of-the-art Hebrew
language model based on the BERT architecture (Devlin et al., 2019), which was
trained on a diverse set of Hebrew texts. Knesset-DictaBERT is the resulting
fine-tuned model: it is specifically tailored for Hebrew parliamentary text.
In this paper we describe the training process of this model, provide a
detailed evaluation of its performance, and demonstrate its superiority over
the baseline DictaBERT model on the Knesset data. We believe that Knesset-
DictaBERT will be a valuable resource for researchers and other users working
on Hebrew language processing and political text analysis.

##  2 Methodology

We fine-tuned the pre-trained DictaBERT model, a state of the art language
model for Hebrew, specifically for the masked language modeling (MLM) task, on
the full Knesset Corpus dataset (Goldin et al., 2024). The model was
initialized using the pre-trained weights from the dicta-il/dictabert
checkpoint available on Hugging Face.

###  2.1 Dataset and Data preprocessing

We used the Knesset Corpus, which is is a large Hebrew dataset of Israeli
parliamentary proceedings, as the dataset for fine-tuning DictaBERT. The
corpus contains over 32M sentences, over 384M tokens and comprises texts from
both plenary and committee protocols.

The corpus was preprocessed to create text shards for efficient loading and
processing. We split the dataset into training, validation, and test sets with
ratios of 80%, 10%, and 10%, respectively.

First, we tokenized the input text using the AutoTokenizer from Hugging Face’s
transformers library. Each text sample was tokenized into sequences of tokens,
with each token represented by its corresponding ID from the tokenizer’s
vocabulary. To prepare the data for the MLM training, we grouped the tokenized
texts into fixed length chunks of 256 tokens. The last chunk was padded to
ensure consistent input sizes for the model. We used the
DataCollatorForLanguageModeling from the transformers library to dynamically
mask tokens during training, with a masking probability of 15%.

###  2.2 Training procedure and hyperparameters

We trained the model on a SLURM (Simple Linux Utility for Resource Management)
environment. We utilized a distributed training setup with the NCCL backend to
leverage multiple GPUs, ensuring efficient training and gradient
synchronization.

We used the following training configuration: We used a per-device batch size
of 32. In order to effectively double the batch size without increasing memory
usage, we set the gradient accumulation steps to 4. A learning rate of
$1e^{-4}$ was chosen for the AdamW optimizer and a weight decay of 0.01 was
applied to regularize the model. The model was trained for 2 full epochs. We
enabled ‘fp16’ (mixed precision) to speed up training and reduce memory usage.
Mixed precision training uses 16-bit floating-point numbers instead of the
standard 32-bit, which can significantly improve computational efficiency and
decrease memory usage Micikevicius et al. (2017). Periodic evaluations were
conducted on the validation set, and the best-performing model checkpoint was
identified based on the validation loss. Final evaluations were performed on
the test set to assess the model’s performance. The model is available on
Hugging Face hub at Knesset-DictaBERT.

##  3 Experiments and Results

The fine-tuned Knesset-DictaBERT model was evaluated on the test set, which
contained about 3.2 million sentences (about 38 million tokens), using
perplexity as the primary metric. The model achieved a perplexity of 6.60,
significantly outperforming the original DictaBERT model, which showed a
perplexity of 22.87. Evidently, the Knesset-DictaBERT model reflects the
language of Knesset proceedings much better than the baseline DictaBERT model.

We also evaluated the fine-tuned Knesset-DictaBERT model for its accuracy in
predicting masked tokens on a subset of the test set containing approximately
300K sentences (about 3.5 million tokens). Since MLM is considered a
challenging task, where multiple options may be equally plausible, the
evaluation focused on three accuracy metrics: top-1 accuracy, top-2 accuracy,
and top-5 accuracy. These metrics measure the model’s ability to correctly
predict the masked token as the first, within the first two, or within the
first five predictions, respectively.

The Knesset-DictaBERT model correctly identified the masked token in the top-1
prediction 52.55% of the time, compared to the original Dicta model, which
achieved a top-1 accuracy of 48.02%. Additionally, Knesset-DictaBERT succeeded
in cases where the original DictaBERT model did not a total of 52,464 times.
In contrast, the original DictaBERT model succeeded where Knesset-DictaBERT
did not in only 27,995 times. Furthermore, when considering the top-2
predictions, Knesset-DictaBERT correctly identified the masked token 63.07% of
the time, whereas the original DictaBERT model had a top-2 accuracy of 58.60%.
Moreover, Knesset-DictaBERT succeeded in 19,400 instances where the original
model failed to provide a correct prediction within the top-2, while the
original DictaBERT model, succeeded in only 13,862 instances where Knesset-
DictaBERT did not. On top of that, when extending the scope to the top-5
predictions, Knesset-DictaBERT demonstrated a significant improvement with a
73.59% accuracy, while the original DictaBERT model achieved a 68.98%
accuracy.

In all tested metrics the Knesset-DictaBERT model outperformed the original
DictaBERT model, indicating a more robust performance in predicting masked
tokens within parliamentary text. These results highlight the effectiveness of
fine-tuning on the specific parliamentary corpus. The results are presented in
Table 1.

Metric | Knesset-DictaBERT | Original DictaBERT
---|---|---
Perplexity | 6.60 | 22.87
Top-1 Accuracy | 52.55% | 48.02%
Top-2 Accuracy | 63.07% | 58.60%
Top-5 Accuracy | 73.59% | 68.98%
Table 1: Comparison of Knesset-DictaBERT and Original DictaBERT on perplexity
and accuracy metrics

##  4 Conclusion and Future Work

In this work, we successfully fine-tuned the DictaBERT model on the Knesset
Corpus to create Knesset-DictaBERT, a model proficient at understanding and
generating parliamentary language in Hebrew. The results indicate a robust
model performance, with substantial improvements over the baseline model.
Future work may involve evaluation on additional Hebrew datasets to enhance
the model’s generalization capabilities and fine-tuning other language models
on the Knesset corpus.

## Limitations

The model was fine-tuned specifically on the Knesset Corpus, which comprises
parliamentary proceedings. As a result, its performance on general Hebrew text
or other domains may not be as robust. However, the original DictaBERT model
was trained on a variety of resources in Hebrew, which probably allows the
Knesset-DictaBERT to still benefit from the diverse linguistic patterns and
vocabulary present in the broader training data of the original model.

## Ethical Considerations

The Knesset Corpus may contain inherent biases, reflecting the political and
social biases present in parliamentary discussions. Consequently, Knesset-
DictaBERT may inherit these biases.

## References

  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186, Minneapolis, Minnesota. Association for Computational Linguistics.
  * Goldin et al. (2024) Gili Goldin, Nick Howell, Noam Ordan, Ella Rabinovich, and Shuly Wintner. 2024.  The Knesset Corpus: An annotated corpus of Hebrew parliamentary proceedings.  _ArXiv_ , abs/2405.18115.
  * Micikevicius et al. (2017) Paulius Micikevicius, Sharan Narang, Jonah Alben, Gregory F. Diamos, Erich Elsen, David García, Boris Ginsburg, Michael Houston, Oleksii Kuchaiev, Ganesh Venkatesh, and Hao Wu. 2017.  Mixed precision training.  _CoRR_ , abs/1710.03740.
  * Shmidman et al. (2023) Shaltiel Shmidman, Avi Shmidman, and Moshe Koppel. 2023.  DictaBERT: A state-of-the-art BERT suite for Modern Hebrew.


# Bayesian Deep Learning Hyperparameter Search for Robust
Function Mapping to Polynomials with Noise

Nidhin Harilal  Indian Institute of Technology Gandhinagar, Gujarat, India
<EMAIL_ADDRESS>Udit Bhatia  Indian Institute of Technology
Gandhinagar, Gujarat, India
<EMAIL_ADDRESS>Auroop R. Ganguly  Northeastern University, Boston, MA,
USA
Pacific Northwest National Laboratory, Richland, WA, USA
<EMAIL_ADDRESS>
###### Abstract

Advances in neural architecture search, as well as explainability and
interpretability of connectionist architectures, have been reported in the
recent literature. However, our understanding of how to design Bayesian Deep
Learning (BDL) hyperparameters, specifically, the depth, width and ensemble
size, for robust function mapping with uncertainty quantification, is still
emerging. This paper attempts to further our understanding by mapping Bayesian
connectionist representations to polynomials of different orders with varying
noise types and ratios. We examine the noise-contaminated polynomials to
search for the combination of hyperparameters that can extract the underlying
polynomial signals while quantifying uncertainties based on the noise
attributes. Specifically, we attempt to study the question that an appropriate
neural architecture and ensemble configuration can be found to detect a signal
of any n-th (where n $\in$ N) order polynomial contaminated with noise having
different distributions and signal-to-noise (SNR) ratios and varying noise
attributes. Our results suggest the possible existence of an optimal network
depth as well as an optimal number of ensembles for prediction skills and
uncertainty quantification, respectively. However, optimality is not
discernible for width, even though the performance gain reduces with
increasing width at high values of width. Our experiments and insights can be
directional to understand theoretical properties of BDL representations and to
design practical solutions.

††preprint: APS/123-QED

##  I Introduction

Neural Networks (NNs) or connectionist representations were originally
inspired by the human brain [1], while feedforward NNs or MultiLayer
Perceptrons (MLPs) were later shown to act as universal function approximators
[2, 3, 4]. However, recent literature points to the imperfect nature of
biological analogies for NNs [5] and the “unreasonable effectiveness” of deep
learning [6], or deep NN representations. Bayesian methods for uncertainty
quantification (UQ) have been suggested for both shallow [7, 8, 9, 10] and
deep [11, 12, 13] NNs. Despite recent successes of connectionist architectures
[14, 15], especially deep learning NNs [16, 17, 6] including Bayesian Deep
Learning (BDL) [10, 18], major gaps remain in our theoretical understanding
and in the design of practical solutions. Deep learning representations, in
particular, appear at first glance to defy the principle of Occam’s Razor or
model parsimony, even though Bayesian [11, 10, 19, 20] or physics-guided [21,
22, 23] approaches may be viewed as constraining the plausible hypotheses
space. Given the simplifying assumptions often made to establish NN (including
deep learning) theory and the ad hoc nature of most engineering solutions, a
complementary approach may be rigorous design-of-experiments with simulated
data. Here we design a set of experiments to understand the function
approximation capability of NNs including deep representations. We map a set
of NNs, specifically (shallow and deep) MLPs, to a set of polynomials
contaminated with noise. The mapping is explored keeping in mind that both NNs
and polynomials are universal function approximators (UFAs) in principle. We
simulate data by varying the degree of the polynomials along with the type of
noise and the signal-to-noise ratios (SNR). The hyperparameters of the NN
(i.e., MLP) representations (i.e., depth and width for the connectionist
function representations and ensemble size for UQ) are examined to
characterize the robustness of the fit in terms of the ability to recover the
original polynomial (measured through prediction skill on test data) and the
noise attributes (measured through distributional distance metrics).

Cybenko [24], Funhashi [25] and Hornik et al. [2] proved that a finite linear
sum of continuous sigmoidal functions could approximate any function to a
desired degree of accuracy, or in other words, that MLPs act as UFAs.
Subsequently, numerous results have shown the UFA property of NNs for
different function classes [3, 26, 4, 27, 28]. Mathematically, this universal
approximation property of neural networks can be described as: Given a
continuous target goal function $f(x)$, there exists a output function $g(x)$
in a linear summation form:

$\displaystyle g(x)=\Sigma_{i=1}^{n_{h}}w_{2,i}\sigma(W_{1,i}x+b_{i})$ (1)
$\displaystyle\text{where}\;\sigma(x)=\frac{1}{1+e^{-x}}$ (2)

and $w_{1,i}$, $w_{2,i}$ are the weights between the first and the output
layers respectively and $b_{i}$ are the bias values associated with the layer
such that $|g(x)-f(x)|\leq\epsilon$ for all $x$ where $\epsilon$ is an
arbitrarily small number.

The UFA property of NNs show that, under various regularity assumptions and
given sufficient data, they can approximate a function $f(t)$ to any desired
degree of accuracy. This implies that for an adequate number of hidden layers
($l$) and nodes ($n_{h}$), there exists a set of corresponding weights and
biases to achieve function approximation to any desired level of accuracy.
However, the theory does not prescribe what the $l$ and $n_{h}$ should be,
thus imposing a significant practical challenge. Also, the non-unique trained
parameters of such NNs are difficult to interpret or explain. Given the
growing complexity of deep learning models and the associated computational
challenges, model parsimony and neural architecture search have become crucial
research areas.

The possible existence of an optimal depth in NNs has been explored in NN
hyperparameter search. This body of literature [29, 30, 31] shows the
existence of certain deep ReLU networks that cannot be realized form shallow
networks. Zhou et al. [32] on the other hand, explores how width affects the
expressiveness of neural networks and shows the existence of classes of wide
networks which cannot be realized by any narrow network whose depth is no more
than a polynomial bound. Poggio et al. [33] provide theorems and examples of a
class of compositional functions for which there is a gap between
approximation in shallow and deep networks.

Theoretical papers on ANNs including MLPs and DL, have focused on performance
guarantees often based on simplified assumptions [29, 30, 31, 28]. However,
what has received limited attention in the literature is the influence of
hyperparameter selection on robust out-of-sample generalization and
uncertainty characterization. One of the key challenges in the analysis of
generalization bounds in deep neural networks is that it may vary depending on
the data distribution on which networks is trained. Therefore, such an
exploration requires extensive empirical analysis. It is however, not
practically possible to experiment on all possible distributions that may
exist. Therefore, the analysis needs to be restricted to certain family of
distributions. Itay et al. [34] provides a theoritical and empirical analysis
to provide several new depth-based separation results on natural radial non-
linear functions such as balls and ellipses. However, their empirical analysis
is restricted to just two depth values.

Recent attempts at UQ on NNs rely on what have been called Bayesian approaches
[35] and have taken the form of so-called BDL [36, 11, 20]. One way of
incorporating Bayesian inferencing in neural networks relies on ensembles
developed through random selection of nodes via a “dropout” strategy (MC-
dropout based networks) [11]. The complexity of conventional neural network
representations for point predictions, along with the heuristic nature of BDL
for uncertainty quantification, implies that any specific DL (or BDL) based
function approximation need to be carefully examined by the ability to
delineate and distinguish between what may be viewed as the signal (with
repeatable or generalizable patterns that may be deterministic) versus noise
(which is not repeatable, indeed it is usually modeled as a stochastic
process).

In this paper, we examine the ability of MC-dropout based neural networks to
distinguish between signal and noise in polynomials of different orders
contaminated with different levels and types of uncorrelated random samples.
Our underlying hypothesis is that the DL-based point prediction can capture
the underlying order of the polynomial while the BDL-based uncertainty
quantification can delineate and capture the statistical attributes of the
noise. A second associated hypothesis - based on the UFA properties of both
NNs and polynomials - is that the behavior exhibited by polynomials
contaminated with noise can be expressed through NNs (specifically, MLPs) with
different hyperparameters. The simulation-based experimental design examine
the hypotheses via metrics for skills in prediction and UQ.

##  II Experimental Design

In this section, we discuss the overview of the experiments performed to
explore the performance of BDL in terms of modeling the underlying polynomial
and approximating noise attributes from data sets of noise-contaminated
polynomials.

###  II.1 Polynomial Dataset

We consider polynomials of different orders and coefficients and add different
types of noises with varied SNR with an understanding that polynomials are
Universal Function Approximators (UFAs) themselves [37]. We use 3 types of
noises: (a) Gaussian, (b) Exponential, ad (c) Rayleigh with SNR ranging from
10 to 30. Mathematically, this polynomial dataset ($P$) can be represented as
follows in Eq. (3):

$\displaystyle P=p(x,n)+\epsilon_{D(t,r)}\;\;\text{where,
}\,p(x,n)=\Sigma_{i=0}^{n}a_{i}x^{i}$ (3)

where $\epsilon_{D(t,snr)}$ represents the noise from distribution $t$ with
SNR level $r$. Table 1 shows the details of various attributes of attributed
dataset.

Table 1: Polynomial dataset description  Attribute | Types
---|---
Poly. Order ($n$) | 2, 3, 4, 5, 7, 10
Noise Types ($t$) | Gaussian111$f(x;\mu;\sigma)=\frac{1}{\sigma\sqrt{2\pi}}e^{-(x-\mu)^{2}/(2\sigma^{2})}$, Exponential222$f(x;\lambda)=\begin{cases}\lambda e^{-\lambda x}&x\geq 0\\\ 0&\text{otherwise}\end{cases}$ , Rayleigh333$f(x;\sigma)=\frac{x}{\sigma^{2}}e^{-x^{2}*(2\sigma^{2})}$
SNR Levels ($r$) | 10, 20 30

###  II.2 Neural Network Design

To understand how different levels of abstraction in neural networks affect
their representation power, we vary the following attributes in our
experiments: (a) number of hidden layers; (b) nodes in each layer; and (c)
number of ensembles.

We follow a bottom-up approach for designing the network configuration for our
experiments. That is, we start with a 1-layer neural network with six neurons
and gradually increase the width and depth of the network and each time we
initialize the network with random weights. Consider width $w\in\mathbb{N}$,
depth $d\in\mathbb{N}$ and number of ensembles $m\in\mathbb{N}$, we represent
$G(w,d,m)$ for a feedforward neural network with $d$-layers each having
$w$-neurons and the network having $m$-ensembles. Table 2 shows the the
considered set of values of different hyperparamters. We analyze the behaviour
of $G(w,d,m)$ on our synthetically generated noisy-polynomial data. This is
done by training each candidate network $G(w,d,m)$ on polynomial data with
different order and noise levels. Details regarding these polynomial orders
and noise types are given in Table 1. This experiment is repeated until all
possible hyperparamter combinations have been trained on all synthetically
generated data.

Table 2: Network Configurations  Parameter | Types
---|---
Depth | 1, 2, 3, 4, …, 15
Width | 6, 8, 10 16, 20, 30, …, 1000
No. of Ensembles | 1, 5, 10, 20, 30, …,40

###  II.3 Evaluation

This section discusses different criteria used for evaluating the degree of
approximation and robustness of the trained networks.

L2-norm: We use L2-norm to evaluate the differences between values predicted
by each of the models and the actual values. L2-norm represents the square
root of the second sample moment of the differences between predicted values
and observed values. For the predicted-vector $\hat{y}$ and actual output-
vector $y$, L2-norm is computed as follows in Eq. (4):

$\displaystyle||\hat{y}-y||_{2}=\left(\Sigma_{i=1}^{N}|\hat{y}-y|^{2}\right)$
(4)

Figure 1:  Experiment pipeline showing various steps involved to understand
relationship between hyper-parameters and model performance for various
network configurations.

Bhattacharyya distance: Bhattacharyya distance is one way of measuring the
similarity between two probability distributions. In our case, we use
Bhattacharyya distance to measure the similarity between the noise
distribution ($p_{1}$) present in the data and the residual ($p_{2}$) obtained
form the prediction. Bhattacharyya distance ($BD$) is calculated as follows in
Eq. (5):

$\begin{split}BD(p_{1},p_{2})=\frac{1}{4}ln\left(\frac{1}{4}\left(\frac{\sigma_{p1}^{2}}{\sigma_{p2}^{2}}+\frac{\sigma_{p2}^{2}}{\sigma_{p1}^{2}}+2\right)\right)\\\
+\frac{1}{4}\left(\frac{(\mu_{p2}-\mu_{p1})^{2}}{\sigma_{p2}^{2}+\sigma_{p1}^{2}}\right)\end{split}$
(5)

where $p_{2}=\hat{y}-y$ and $p_{1}=\epsilon_{D(t,r)}$ represents the noise
from distribution $t$ and SNR-level $r$.

While L2-norm gives an estimation of how close the predictions are to the
actual values, Bhattacharyya distance between the residual and noise, on the
other hand, gives an estimate of how much signal has been retained and what
amount of bias has been induced on the network due to the added noise.

##  III Results and Analysis

This section discusses the affect of varying the network size, the function
complexity, and the amount of noise added to the training data on performance
criteria including L2-norm and Bhattacharyya distance.

Figure 2 (A-C) shows the test-set error (L1-norm) as the number of hidden
units (width w) is varied from 6 to 140, and the depth ($d$) of the network
from 1 to 15. We observe better performance in networks with larger width.

Figure 2:  Top row: Test Loss (L1-norm) with increasing polynomial degree on
Gaussian noise with $SNR=20$. Bottom row: Test Loss (L1-norm) on networks with
width ($w$) ranging from 6 to 1000.

Increasing Width: We increase the range of width ($w$) to 1000 to investigate
the existence of ”optimality” (used loosely in this research) in terms of
width. Figure 2 (D) shows the test loss L1-norm on polynomial with degree 3
and right-side for degree 5. Both the polynomials are contaminated with
Gaussian noise having $SNR=20$. As we move in the direction of increasing
width (red arrow), we notice constant decrements in test loss values. Also,
notice the saturation of test loss on networks when increasing the width. We
note the high level of test-set loss in case of 1-layer network. A possible
explanation is that a single-layer is probably too small to accurately
characterize the target function.

Increasing Depth: For each value of width, we train networks increase with
depth values ranging from 1 to 15 to also assess the impact of depth
parameter. From the right heatmap in Figure 2 (E), as we move in the direction
of increasing width (blue arrow), we observe decreasing values of test loss
followed by a steady increase in test loss indicating deteriorating
performance while incrementing after a certain depth-value in the feedforward
neural network.

Figure 3:  A, B: Test Loss (L1-norm) and Bhattacharya distance criterion
landscapes on polynomial with degree 7 on exponential noise with $SNR=10$. C:
Test Loss (L1-norm) criterion landscape and individual loss breakdown along
width on $5^{th}$ degree polynomial data with exponential noise and $SNR=10$.

###  III.1 Loss Landscapes

Figure 3 (A, B) shows the discussed loss landscapes of L1-norm (A) and
Bhattacharya norm (B) criteria on out-of-distribution data points of the
polynomial with degree 7 on exponential noise with $SNR=10$. If we consider
the plot in Figure 3 (A), it can be observed that there is a constant downward
slope in the landscape as width increases (yellow arrow). On the other hand,
when we consider models with higher depth, we first observe a steep downward
slope followed by an upward slope leading to a valley-like structure. A
corresponding one-to-one scenario is observed in landscapes based on
Bhattacharya distance criterion. The only difference is that the surface
becomes upside down compared to the curvature observed in the case of L1-loss.
Figure 3 (C) shows a surface breakdown of the L1-loss values along with the
depth values, which allows for a more apparent observation of decreasing L-1
test values of candidates with higher width values.

Based on the results obtained, it can be inferred that there exists no
optimality in terms of width in neural network on both L-1 loss and
Bhattacharya distance criterion. Since, an increment in performance for the
task of polynomial approximation is observed on models with higher width,
therefore, a higher width is favourable. However, we do see that both L1-loss
values and Bhattacharya distance criterion values starts saturating, we can
consider it a performance gain vs model complexity tradeoff along the width.
Therefore, on the basis of a fixed upper-threshold of model complexity or
computational overheads, considering a model with maximum permissible width
value is expected to perform better than models with lower width values given
other hyperparameters are fixed.

While the analysis of candidates along increasing depth showed an increase
followed by decrease in performance values in both the criterions indicating
the existence of an optimal value of depth located in the region of changing
slope direction. However, this depth value is not universal across all the
width values, which means that the exact value of the optimal depth depends on
the value of width chosen. Mathematically, we can say that, for a given width
$w’\exists$ depth $d’$ such that $d’$ is optimal for all possible
$G(d,w’,m),d\geq 1$

###  III.2 Effect of Ensembles

Using the Monte-Carlo dropout method as a Bayesian approximation [11], we
consider our analysis on the third parameter, that is, the number of
ensembles. By allowing dropouts in the test time, we conduct multiple
inferences and consider the average of all the obtained models as the final
output. Furthermore, the dropout strategy ensures that each time a different
model is obtained with high probability, which helps in bringing diversity and
thus reducing variance in the final model output. Mathematically, this can be
shown as:

$\displaystyle\hat{y}=\frac{1}{N}\Sigma_{i}^{m}f_{i}{x}$ (6)

$\begin{split}E^{2}=(y-\hat{y})^{2}=\frac{1}{N}\Sigma_{i}^{m}(y-f_{i}(x))^{2}\\\
-\frac{1}{N}\Sigma_{i}^{m}(f_{i}(x)-\hat{y})^{2}\end{split}$ (7)

where $\hat{y}$, $y$ and $f_{i}$ represents the final model output, the ground
truth and the $i^{th}$ ensemble member respectively. The breakdown of
$(y-\hat{y})^{2}$ in two parts shows why ensembles work better than a single
model. The first term is the actual error whereas second term represents the
disagreement between ensembles.

In comparison to the Mean-squared-error (MSE) for a single model, ensembles
introduces another term (MSE between ensembles and the model output) as shown
in Equation 7 which contributes in reducing the overall MSE of the model [38].

Figure 4:  Bhattacharya criterion values of different ensemble sozes with
increasing network width on 7 degree polynomial and Gaussian noise with
$SNR=20$.

Figure 4 shows the Bhattacharya criterion values of models with different
number of ensembles. A sudden increase in performance is observed as we go
from a single model to an ensemble. But the performance gain on further
increasing the ensembles size is observed to decrease. The reason for this
observation is that the chances of finding uncorrelated models as we increase
the number of ensembles decreases. If a particular ensemble size already
incorporates the top-performing models then further adding other members will
not offer any benefit to the ensemble. From Figure 4, it can be seen that in
many networks, the highest ensemble size is not the best performer. Therefore,
an optimal ensemble size may exist beyond which an improvement in performance
isn’t expected as the ensemble will not able to be harness their contribution
effectively.

##  IV Conclusion and Discussion

We observed interesting relations between the choice of BDL hyperparameters
and corresponding skill metrics for predictions with uncertainty. While
directly relevant for the approximation of noisy polynomials, the insights may
be directional to explore BDL-based robust function approximation in wider
settings.

The surface curvature obtained in our results showed that as depth increases
past an optimal point, generalization performance tends to decrease. This
experimental determination of the existence of an optimal depth value in
polynomial function approximation tasks can allow for reduced number of trials
needed to find the best BDL model in practical settings.

The observation of a continuously decreasing but positive performance gain
with increasing width in our results indicated that an optimal width may not
exist. However, the gain in generalization performance beyond a certain high
value of width were observed to grow less statistically significant. These set
of results suggest that the width may be chosen based on a threshold related
to model complexity (e.g., a modified information criteria where model
complexity and performance gain on out-of-sample data may need to be balanced)
or based on the available computational resources.

We observed, through distributional distance metrics, that an optimal is
suggested for ensemble sizes in a dropout-based BDL. In other words, the best
approximations to the noise statistics (used to contaminate the underlying
polynomials) were obtained for certain optimal ensemble sizes. These set of
insights may be useful in the practical design and implementation of BDLs.

Our results point to the existence of an optimal depth and an optimal ensemble
size but no optimal width for BDL representations. The empirical insights
presented here may benefit from a closer relation to existing and potentially
new learning theory in these areas. Future experiments on simulated data need
to examine a broader set of simulations, including different types of
(potentially nonlinear) signals as well as (potentially correlated) noise
processes, in both time and space. Practical BDL guidelines must be developed
across multiple science, engineering, and business domains, by considering use
cases based on both realistic simulations as well as on real data.

###### Acknowledgements.

ARG was supported by four National Science Foundation Projects including NSF
BIG DATA under Grant No. 1447587, NSF Expedition in Computing under Grant No.
1029711, NSF CyberSEES under Grant No. 1442728, and NSF CRISP type II under
Grant No. 1735505. UB was supported by the Ministry of Education (India),
STARS under project Project ID: 367. ARG acknowledges his guest affiliation at
IIT Gandhinagar which made this work feasible.

## References

  * Amit and Amit [1992] D. J. Amit and D. J. Amit,  _Modeling brain function: The world of attractor neural networks_ (Cambridge university press, 1992).
  * Hornik  _et al._ [1989] K. Hornik, M. Stinchcombe, H. White, _et al._ , Multilayer feedforward networks are universal approximators., Neural networks 2, 359 (1989).
  * Hornik [1991] K. Hornik, Approximation capabilities of multilayer feedforward networks, Neural networks 4, 251 (1991).
  * Barron [1993] A. R. Barron, Universal approximation bounds for superpositions of a sigmoidal function, IEEE Transactions on Information theory 39, 930 (1993).
  * Lillicrap  _et al._ [2020] T. P. Lillicrap, A. Santoro, L. Marris, C. J. Akerman, and G. Hinton, Backpropagation and the brain, Nature Reviews Neuroscience , 1 (2020).
  * Sejnowski [2020] T. J. Sejnowski, The unreasonable effectiveness of deep learning in artificial intelligence, Proceedings of the National Academy of Sciences 117, 30033 (2020).
  * MacKay [1995a] D. J. MacKay, Bayesian neural networks and density networks, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 354, 73 (1995a).
  * Neal [1993] R. M. Neal, Bayesian learning via stochastic dynamics, in  _Advances in neural information processing systems_ (1993) pp. 475–482.
  * Neal [2012] R. M. Neal,  _Bayesian learning for neural networks_ , Vol. 118 (Springer Science & Business Media, 2012).
  * Kendall and Gal [2017] A. Kendall and Y. Gal, What uncertainties do we need in bayesian deep learning for computer vision?, arXiv preprint arXiv:1703.04977 (2017).
  * Gal and Ghahramani [2016] Y. Gal and Z. Ghahramani, Dropout as a bayesian approximation: Representing model uncertainty in deep learning, in  _international conference on machine learning_ (PMLR, 2016) pp. 1050–1059.
  * Gal  _et al._ [2017] Y. Gal, J. Hron, and A. Kendall, Concrete dropout, arXiv preprint arXiv:1705.07832 (2017).
  * Vandal  _et al._ [2018] T. Vandal, E. Kodra, J. Dy, S. Ganguly, R. Nemani, and A. R. Ganguly, Quantifying uncertainty in discrete-continuous and skewed data with bayesian deep learning, in  _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ (2018) pp. 2377–2386.
  * Norris [1994] D. Norris, Shortlist: A connectionist model of continuous speech recognition, Cognition 52, 189 (1994).
  * Christiansen and Chater [1999] M. H. Christiansen and N. Chater, Toward a connectionist model of recursion in human linguistic performance, Cognitive Science 23, 157 (1999).
  * Krizhevsky  _et al._ [2012] A. Krizhevsky, I. Sutskever, and G. E. Hinton, Imagenet classification with deep convolutional neural networks, Advances in neural information processing systems 25, 1097 (2012).
  * LeCun  _et al._ [2015] Y. LeCun, Y. Bengio, and G. Hinton, Deep learning, nature 521, 436 (2015).
  * Khan  _et al._ [2018] M. Khan, D. Nielsen, V. Tangkaratt, W. Lin, Y. Gal, and A. Srivastava, Fast and scalable bayesian deep learning by weight-perturbation in adam, in  _International Conference on Machine Learning_ (PMLR, 2018) pp. 2611–2620.
  * Vasudevan  _et al._ [2021] R. K. Vasudevan, M. Ziatdinov, L. Vlcek, and S. V. Kalinin, Off-the-shelf deep learning is not enough, and requires parsimony, bayesianity, and causality, npj Computational Materials 7, 1 (2021).
  * Luo and Kareem [2020] X. Luo and A. Kareem, Bayesian deep learning with hierarchical prior: Predictions from limited and noisy data, Structural Safety 84, 101918 (2020).
  * Jia  _et al._ [2021] X. Jia, J. Willard, A. Karpatne, J. S. Read, J. A. Zwart, M. Steinbach, and V. Kumar, Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles, ACM/IMS Transactions on Data Science 2, 1 (2021).
  * Jia  _et al._ [2019] X. Jia, J. Willard, A. Karpatne, J. Read, J. Zwart, M. Steinbach, and V. Kumar, Physics guided rnns for modeling dynamical systems: A case study in simulating lake temperature profiles, in  _Proceedings of the 2019 SIAM International Conference on Data Mining_ (SIAM, 2019) pp. 558–566.
  * Li and Saad [2018] B. Li and D. Saad, Exploring the function space of deep-learning machines, Physical review letters 120, 248301 (2018).
  * Cybenko [1989] G. Cybenko, Approximation by superpositions of a sigmoidal function, Mathematics of control, signals and systems 2, 303 (1989).
  * Funahashi [1989] K.-I. Funahashi, On the approximate realization of continuous mappings by neural networks, Neural networks 2, 183 (1989).
  * Park and Sandberg [1991] J. Park and I. W. Sandberg, Universal approximation using radial-basis-function networks, Neural computation 3, 246 (1991).
  * Leshno  _et al._ [1993] M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, Multilayer feedforward networks with a nonpolynomial activation function can approximate any function, Neural networks 6, 861 (1993).
  * Bengio and Delalleau [2011] Y. Bengio and O. Delalleau, On the expressive power of deep architectures, in  _International conference on algorithmic learning theory_ (Springer, 2011) pp. 18–36.
  * Eldan and Shamir [2016] R. Eldan and O. Shamir, The power of depth for feedforward neural networks, in  _Conference on learning theory_ (PMLR, 2016) pp. 907–940.
  * Cohen  _et al._ [2016] N. Cohen, O. Sharir, and A. Shashua, On the expressive power of deep learning: A tensor analysis, in  _Conference on learning theory_ (PMLR, 2016) pp. 698–728.
  * Liang and Srikant [2016] S. Liang and R. Srikant, Why deep neural networks for function approximation?, arXiv preprint arXiv:1610.04161 (2016).
  * Lu  _et al._ [2017] Z. Lu, H. Pu, F. Wang, Z. Hu, and L. Wang, The expressive power of neural networks: A view from the width, arXiv preprint arXiv:1709.02540 (2017).
  * Poggio  _et al._ [2017] T. Poggio, H. Mhaskar, L. Rosasco, B. Miranda, and Q. Liao, Why and when can deep-but not shallow-networks avoid the curse of dimensionality: a review, International Journal of Automation and Computing 14, 503 (2017).
  * Safran and Shamir [2017] I. Safran and O. Shamir, Depth-width tradeoffs in approximating natural functions with neural networks, in  _International Conference on Machine Learning_ (PMLR, 2017) pp. 2979–2987.
  * Izmailov  _et al._ [2020] P. Izmailov, W. J. Maddox, P. Kirichenko, T. Garipov, D. Vetrov, and A. G. Wilson, Subspace inference for bayesian deep learning, in  _Uncertainty in Artificial Intelligence_ (PMLR, 2020) pp. 1169–1179.
  * Wang and Yeung [2016] H. Wang and D.-Y. Yeung, Towards bayesian deep learning: A framework and some existing methods, IEEE Transactions on Knowledge and Data Engineering 28, 3395 (2016).
  * Pinkus [2000] A. Pinkus, Weierstrass and approximation theory, Journal of Approximation Theory 107, 1 (2000).
  * Hansen and Salamon [1990] L. K. Hansen and P. Salamon, Neural network ensembles, IEEE transactions on pattern analysis and machine intelligence 12, 993 (1990).
  * Blumer  _et al._ [1987] A. Blumer, A. Ehrenfeucht, D. Haussler, and M. K. Warmuth, Occam’s razor, Information processing letters 24, 377 (1987).
  * Wang  _et al._ [2020] S. Wang, T. H. McCormick, and J. T. Leek, Methods for correcting inference based on outcomes predicted by machine learning, Proceedings of the National Academy of Sciences 117, 30266 (2020).
  * Rasmussen and Ghahramani [2000] C. Rasmussen and Z. Ghahramani, Occam’s razor, Advances in neural information processing systems 13, 294 (2000).
  * MacKay [1995b] D. J. MacKay, Probable networks and plausible predictions—a review of practical bayesian methods for supervised neural networks, Network: computation in neural systems 6, 469 (1995b).
  * Karnin [1990] E. D. Karnin, A simple procedure for pruning back-propagation trained neural networks, IEEE transactions on neural networks 1, 239 (1990).
  * Fahlman and Lebiere [1989] S. Fahlman and C. Lebiere, The cascade-correlation learning architecture, Advances in neural information processing systems 2, 524 (1989).
  * Kossaifi  _et al._ [2017] J. Kossaifi, A. Khanna, Z. Lipton, T. Furlanello, and A. Anandkumar, Tensor contraction layers for parsimonious deep nets, in  _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition Workshops_ (2017) pp. 26–32.
  * Lawrence  _et al._ [1998] S. Lawrence, C. L. Giles, and A. C. Tsoi,  _What size neural network gives optimal generalization? Convergence properties of backpropagation_ , Tech. Rep. (1998).
  * Malach and Shalev-Shwartz [2019] E. Malach and S. Shalev-Shwartz, Is deeper better only when shallow is good?, in  _Advances in Neural Information Processing Systems_ (2019) pp. 6429–6438.
  * Bekele and Lawson [2019] E. Bekele and W. Lawson, The deeper, the better: analysis of person attributes recognition, in  _2019 14th IEEE International Conference on Automatic Face & Gesture Recognition (FG 2019)_ (IEEE, 2019) pp. 1–8.
  * Becker  _et al._ [2020] S. Becker, Y. Zhang, _et al._ , Geometry of energy landscapes and the optimizability of deep neural networks, Physical review letters 124, 108301 (2020).
  * Geman  _et al._ [1992] S. Geman, E. Bienenstock, and R. Doursat, Neural networks and the bias/variance dilemma, Neural computation 4, 1 (1992).
  * Belkin  _et al._ [2019] M. Belkin, D. Hsu, S. Ma, and S. Mandal, Reconciling modern machine-learning practice and the classical bias–variance trade-off, Proceedings of the National Academy of Sciences 116, 15849 (2019).
  * Molchanov  _et al._ [2019] P. Molchanov, A. Mallya, S. Tyree, I. Frosio, and J. Kautz, Importance estimation for neural network pruning, in  _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ (2019) pp. 11264–11272.
  * Moody [1994] J. Moody, Prediction risk and architecture selection for neural networks, in  _From statistics to neural networks_ (Springer, 1994) pp. 147–165.
  * Elsken  _et al._ [2018] T. Elsken, J. H. Metzen, and F. Hutter, Neural architecture search: A survey, arXiv preprint arXiv:1808.05377 (2018).
  * Lampinen and Vehtari [2001] J. Lampinen and A. Vehtari, Bayesian approach for neural networks—review and case studies, Neural networks 14, 257 (2001).
  * Loquercio  _et al._ [2020] A. Loquercio, M. Segu, and D. Scaramuzza, A general framework for uncertainty estimation in deep learning, IEEE Robotics and Automation Letters 5, 3153 (2020).
  * Valentini and Masulli [2002] G. Valentini and F. Masulli, Ensembles of learning machines, in  _Italian workshop on neural nets_ (Springer, 2002) pp. 3–20.


# Spin-phonon decoherence in solid-state paramagnetic defects from first
principles

Sourav Mondal  Alessandro Lunghi<EMAIL_ADDRESS>School of Physics, AMBER and
CRANN Institute, Trinity College, Dublin 2, Ireland

###### Abstract

Paramagnetic defects in diamond and hexagonal boron nitride possess a unique
combination of spin and optical properties that make them prototypical solid-
state qubits. Despite the coherence of these spin qubits being critically
limited by spin-phonon relaxation, a full understanding of this process is not
yet available. Here we apply ab initio spin dynamics simulations to this
problem and quantitatively reproduce the experimental temperature dependence
of spin relaxation time and spin coherence time. We demonstrate that low-
frequency two-phonon modulations of the zero-field splitting are responsible
for spin relaxation and decoherence, and point to the nature of vibrations in
2-dimensional materials as the culprit for their shorter coherence time. These
results provide a novel interpretation to spin-phonon decoherence in solid-
state paramagnetic defects, offer a new strategy to correctly interpret
experimental results, and pave the way for the accelerated design of new spin
qubits.

## Introduction

Defects in solid-state semiconductors often introduce additional electronic
states with energy lower than the band-gap, leading to color centers. Tens of
different color centers are known for only diamondalkahtani2018fluorescent
and Silicon Carbidecastelletto2020silicon  and their presence often enriches
the original material’s optical and magnetic properties, enabling interesting
applications in the fields of sensing, photonics, and more. In particular, the
negative nitrogen-vacancy (NV-) center in diamond consists, represented in
Fig. 1, of a spin-1 system, with spin density localized at the defect site.
NV- centers hold a special role in the quantum technology ecosystem thanks to
some ideal properties, including long relaxation and coherence times ($T_{1}$
and $T_{2}$, respectively) at ambient temperaturebalasubramanian2009ultralong
; bar2013solid , chemical and mechanical robustnessschirhagl2014nitrogen , and
the possibility to optically address and initialize the qubit
statesgruber1997scanning ; jelezko2004observation . NV- centers have thus
found immediate application in quantum technologies, going from
sensinggrotz2011sensing ; tetienne2013spin ; schirhagl2014nitrogen  to quantum
communicationstogan2010quantum . On the verge of NV- success, other materials
with similar properties have been individuated leading to a fast-growing
family of solid-state spin qubitswolfowicz2021quantum . This includes negative
boron vacancies in hexagonal boron-nitride ($\mathrm{V_{B}^{-}}$), reported in
Fig. 1, which has emerged as a promising 2D-material-based qubit with an
optimal spin-optical interface and mechanical
propertiesgottscholl2020initialization ; gottscholl2021room ; stern2022room .

Figure 1: Defects’ Structure. Left panel: The structure of the defects NV- in
the unit cell of diamond. Right panel: the structure of
V${}_{\mathbf{B}}^{-}1$ embedded on a 6x6 supercell of h-BN. Carbon atoms are
ported in grey, nitrogen atoms in blue, and boron atoms in pink. The vacancy
is highlighted by the presence of a fictitious yellow atom.

Spin coherence is key for any quantum application and ultimately sets the
limit of sensing accuracy or computation fidelity. Despite its central role in
the physics of paramagnetic defects, the contribution of spin-phonon
interaction to relaxation and decoherence has not yet been fully understood
and experiments are invariably interpreted by means of phenomenological models
based on a simplistic Debye picture of phononswalker19685 ;
wolfowicz2021quantum . Such a state of affairs effectively prevents the
establishment of a rigorous understanding of spin dynamics in solid-state
qubits and several outstanding questions are left unanswered. For instance, it
is yet not clear what mechanism of spin relaxation is the driving one at non-
cryogenic temperatures, with one- (Direct and Orbach) and two-phonon (Raman)
mechanisms often interchangeably cited as possible relaxation
routesjarmola2012temperature ; gottscholl2021room . Moreover, spin-phonon
relaxation has been shown to limit spin coherence at ambient temperature, but
with significant deviations from the commonly expected limit
$T_{2}=2T_{1}$bar2013solid , suggesting that additional vibronic contributions
to decoherence might play an important role. Besides the importance of such
questions from a fundamental science point of view, the current limits to our
understanding of the microscopic processes determining spin coherence in
paramagnetic solid-state materials pose serious challenges for their use as
quantum sensors, namely one of the main applications for such compounds. In a
typical experiment, the information on the spin’s environment is obtained
through the study of spin-environment fingerprints present in the spin
dynamicsdegen2017quantum . Such protocol is inherently connected to our
ability to build reliable and quantitative models for spin-environment
interactions and dynamics, of which spin-phonon coupling and relaxation are
excellent examples. Last but not least, from a materials science point of
view, the lack of robust structure-property relations for these compounds
makes the identification of new paramagnetic defects with desired properties
an outstanding challenge, so far mostly left to serendipity.

Figure 2: Phonon density of state and spin-phonon coupling. Results for NV-
and V${}_{\mathbf{B}}^{-}1$ are reported in the left and right panel,
respectively. Red curves report the phonon density of states, while the blue
curves report the spin-phonon coupling density $V(\omega)$. Star symbols point
to the phonon modes displayed as insets.

Ab initio methods have been applied to the problem of spin relaxation in the
context of molecular systems in the last few yearslunghi2019phonons ;
escalera2019exploring ; lunghi2020limit ; reta2021ab ; lunghi2022toward ;
mondal2022unraveling , leading to unprecedented insights into the physics of
spin systems. On the ohter hand, only very recent attempts to apply ab initio
strategies to solid-state defects are availableastner2018solid ; tang2022first
; cambria2022temperature , and a full picture of what is the leading
relaxation mechanism at ambient temperature in both 3D and 2D semiconductors,
the role of phonons on coherence, and what type of vibrations are responsible
for such processes have not yet been fully elucidated. Here we extend the
state-of-the-art in this nascent fieldlunghi2022spin ; lunghi2022toward  to
solid-state compounds and obtain a quantitative microscopic picture of spin
relaxation and decoherence in NV- and $\mathrm{V_{B}^{-}}$. Our results can be
exported to the rest of the increasingly large family of $S=1$ spin
qubitswolfowicz2021quantum ; bayliss2020optically  and have the potential to
provide unprecedented guidance toward the design of new materials with optimal
features for quantum technologies.

## Results

Spin-Phonon Coupling. The ground state of both NV- and V${}_{\mathrm{B}}^{-}$
is well described by a $S=1$ Hamiltonian

$\hat{H}_{S}=\vec{\mathbf{S}}\cdot\mathbf{D}\cdot\vec{\mathbf{S}}+\mu_{\mathrm{B}}\>\vec{\mathbf{S}}\cdot\mathbf{g}\cdot\vec{\mathbf{B}}\>,$
(1)

where the Bohr’s magneton, $\mu_{\mathrm{B}}$, and the Lande’s tensor,
$\mathbf{g}$, mediate the interaction with the external field
$\vec{\mathbf{B}}$, and $\vec{\mathbf{D}}$ is the zero-field splitting matrix.
By expressing this Hamiltonian in a reference frame with the $z$-axis parallel
to the defect’s $C_{3}$ symmetry axis, the first term of Eq. 1 reduces to
$\hat{H}_{S}=DS_{z}^{2}$, with $D\sim 0.118$ cm-1 and $\sim 0.092$ cm-1 for
V${}_{\mathrm{B}}^{-}$, and NV-, respectively, whereas $\mathbf{g}$ closely
matches the free electron value for both compounds. The zero-field splitting
stabilizes the $\hat{S}_{z}$ eigenstate $|0\rangle$ with respect to the pair
$|\pm 1\rangle$ by an energy amount of $D$. Finally, the degeneracy of $|\pm
1\rangle$ is broken by the Zeeman interaction (see Fig. S6). Density
functional theory (DFT) nicely reproduces these experimental values when
applied to a small cluster of atoms around the defect (see SM). This approach
is used to describe the defects’ spin properties throughout this work.

The interaction between the defect’s spin and the phonons is modeled by the
spin-phonon Hamiltonianlunghi2022toward

$\hat{H}_{\mathrm{S-Ph}}=\sum_{\alpha}\left(\frac{\partial\hat{H}_{S}}{\partial
q_{\alpha}}\right)q_{\alpha}+\frac{1}{2}\sum_{\alpha,\beta}\left(\frac{\partial^{2}\hat{H}_{S}}{\partial
q_{\alpha}\partial q_{\beta}}\right)q_{\alpha}q_{\beta}\>,$ (2)

where we assume the phonons $q_{\alpha}$ to be quantum harmonic oscillators of
frequency $\omega_{\alpha}/2\pi$. Here we consider the modulation of the zero-
field splitting $\mathbf{D}$ as the source of spin-phonon coupling and
numerically compute its first and second-order derivatives. In order to reduce
the computational overheads associated with a second-order numerical
differentiation, we trained a neural network on a set of $\mathbf{D}$ values
computed with DFT over 1600 random distortions of the considered
clusterlunghi2020limit ; garlatti2022critical . The neural network is tested
over 400 residual distortions not used during the training and shows a
remarkable ability to predict the values of $\mathbf{D}$ for unseen geometries
with a root-mean-squared error of only 1.8 x 10-4 cm-1 for NV- and 2.6 x 10-4
cm-1 for V${}_{\mathrm{B}}^{-}$, see Figs. S8-S10. The net is then used to
perform a 6- and 36-point numerical differentiation of $\mathbf{D}$ for first-
and second-order, respectively.

Phonons were computed with periodic DFT, and the vibrational density of states
is shown in Fig. 2. In order to individuate the coupling strength of each
phonon we studied the function

$V(\omega)=\sum_{ij}^{3}\left(\frac{\partial D_{ij}}{\partial
q_{\alpha}}\right)^{2}\delta(\omega-\omega_{\alpha})\>.$ (3)

$V(\omega)$ is plotted in Fig. 2, together the phonon density of states, and
shows that strongly coupled vibrations are present over the entire spectrum. A
visual inspection of the strongest coupled mode for NV- shows that spin-phonon
coupling is promoted by either compressive waves of diamonds (mode at
$\sim$350 cm-1 in Fig. 2) and distortions localized at the defect site (mode
at $\sim$500 cm-1 in Fig. 2). The analysis of modes for V${}_{\mathrm{B}}^{-}$
shows that spin-phonon coupling is effectively originated by out-of-plane
modes that locally distort the defect.

Spin-Phonon Relaxation. Spin-phonon relaxation is described by means of
density-matrix perturbation theory, where phonons, namely the perturbation,
act as a thermal reservoir with infinite specific heat and instantaneous
relaxation time. In this framework, the Markovian time-evolution of the
reduced spin density matrix in the interaction picture, $\hat{\rho}_{s}(t)$,
is described bylunghi2022toward

$\frac{d\hat{\rho_{s}}}{dt}=\left(\hat{\hat{\mathbf{R2}}}^{1-\mathrm{ph}}+\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}+\hat{\hat{\mathbf{R4}}}^{2-\mathrm{ph}}\right)\>\hat{\rho}_{s}(t)\>,$
(4)

where $\hat{\hat{\mathbf{R2}}}^{1-\mathrm{ph}}$ describes resonant single-
phonon processes arising from the first term of the Hamiltonian in Eq. 2 and
second-order density matrix perturbation theory. As shown in Fig. 2, phonons’
energies largely exceed the zero-field splitting of the spin qubits, resulting
in the absence of one-phonon processes, namely direct and Orbach relaxation
mechanisms, which would require degeneracy among spin and phonon states. In
principle, low-energy phonons would appear by increasing the size of the cell
used for phonons simulations, thus recovering the long-wavelength acoustic
modes of the lattice. However, as shown for molecular qubitslunghi2019phonons
; garlatti2020unveiling , these modes have vanishingly small spin-phonon
coupling and become relevant to spin dynamics only in the low-temperature
relaxation regimeastner2018solid ; gugler2018ab . Low-temperature spin
relaxation is also strongly affected by cross-relaxation, as shown by
experiments conducted on samples with different defect
concentrationsjarmola2012temperature , and does not represent an ideal regime
to compare simulations and experiments. Therefore, we neglect the contribution
$\hat{\hat{\mathbf{R2}}}^{1-ph}$ to spin dynamics and focus on the more
application-relevant high-temperature limit dominated by two-phonon processes.

The term $\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ accounts for two-phonon
processes leading to Raman relaxation and originates from the quadratic spin-
phonon coupling terms of Eq. 2 at the second-order of perturbation theory. In
this approximation, the dynamics of the diagonal elements of $\hat{\rho}_{s}$
reads

$R2^{2-\mathrm{ph}}_{bb,aa}=\frac{\pi}{4\hbar^{2}}\sum_{\alpha\beta}V^{\alpha\beta}_{ba}V^{\alpha\beta}_{ab}G^{2-\mathrm{ph}}(\omega_{ba},\omega_{\alpha},\omega_{\beta})$
(5)

where $V^{\alpha\beta}_{ab}$ stands for $\langle
a|(\partial^{2}\hat{H}_{\mathrm{s}}/\partial q_{\alpha}\partial
q_{\beta})|b\rangle$. The function $G^{2-\mathrm{ph}}$ accounts for three
possible processes involving two phonons: absorption of two phonons, emission
of two phonons, and simultaneous emission of one phonon and absorption of a
second one. Among all the possible two-phonon processes, two-phonon
absorption/emission is negligible as they require phonons with energy lower
than the zero-field splitting, absent in our framework. Therefore, the
simultaneous emission/absorption of two phonons is the only process that
enables energy-conserving two-phonon spin transitions, for which
$G^{2-\mathrm{ph}}$ reads

$G^{2-\mathrm{ph}}(\omega_{ba},\omega_{\alpha},\omega_{\beta})=\delta(\omega_{ba}-\omega_{\alpha}+\omega_{\beta})\bar{n}_{\alpha}(\bar{n}_{\beta}+1)\>.$
(6)

Finally, the term $\hat{\hat{\mathbf{R4}}}^{2-\mathrm{ph}}$ accounts for those
two-phonon processes due to linear spin-phonon coupling at the fourth-order
perturbation theory and reads

$R4_{bb,aa}^{2-ph}=\frac{\pi}{2\hbar^{2}}\sum_{\alpha\beta}\left|T^{\alpha\beta,+}_{ba}+T^{\beta\alpha,-}_{ba}\right|^{2}G^{2-ph}(\omega_{ba},\omega_{\alpha},\omega_{\beta})\>,$
(7)

where the terms

$T^{\alpha\beta,\pm}_{ba}=\sum_{c}\frac{\langle
b|(\partial\hat{H}_{s}/\partial q_{{\alpha}})|c\rangle\langle
c|(\partial\hat{H}_{s}/\partial
q_{{\beta}})|a\rangle}{E_{c}-E_{a}\pm\hbar\omega_{\beta}}\>,$ (8)

include the contributions of virtual excited states $|c\rangle$.

Figure 3:  Spin-phonon relaxation. The left panel reports the simulated rate
$(T_{1})^{-1}$ (black solid line) as a function of $T$ for NV- defect in
diamond together with the experimental rates (red squares)gottscholl2021room .
The dashed yellow line represents the fitting with Eq. 9. Dashed blue and
green lines are the two separate contributions of Eq. 9. The right panel
reports the simulated rate as a function of $T$ for V${}_{\mathrm{B}}^{-}$
(solid black line) together with experimental rates (red squares). The dashed
yellow line represents the fitting with the function $A\>T^{2}$.

We compute $\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ and
$\hat{\hat{\mathbf{R4}}}^{2-\mathrm{ph}}$ with the spin-phonon coupling
coefficients and phonon frequencies determined by electronic structure
calculations and simulate the relaxation dynamics of both defects with the
software MolForgelunghi2022toward ; MolForge . In accordance with the
experimental setup, we initialize the system in $|1\rangle$ and then monitor
the population recovery of $|0\rangle$ in time in the presence of a field of 1
mTjarmola2012temperature ; gottscholl2021room . As shown in Fig. S5 we observe
a mono-exponential profile that leads to a single relaxation rate constant
$(T_{1})^{-1}$. Fig. 3 reports the simulated $(T_{1})^{-1}$ as function of
temperature obtained by applying $\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ to
both defects. The comparison with experimental data reveals a striking
agreement for both materials, and simulations nicely reproduce both the
$T$-dependence and the trend in relaxation rates. On the other hand,
simulations carried out with $\hat{\hat{\mathbf{R4}}}^{2-\mathrm{ph}}$ lead to
much slower relaxation rates (see Fig. S4), thus leading to the conclusion
that two-phonon relaxation due to quadratic spin-phonon coupling is the most
relevant pathway in solid-state defects. The robustness of our conclusions is
demonstrated by performing simulations with different phonons cell sizes, as
shown in Fig. S2.

The excellent accuracy of simulations validates the computational approach and
allows us to revisit the interpretation of spin relaxation in these systems.
Firstly, our simulations reveal that no one-phonon Orbach process is involved
in the spin relaxation of either NV- and V${}_{\mathrm{B}}^{-}$ defects at
high-temperature, as previously hypothesized in the fitting of experimental
datajarmola2012temperature . Our finding is in agreement with spin energy-
conservation arguments and the low-energy cut-off frequency of the phonons’
density of states. Moreover, the fitting of experimental data accounted for
the Raman contribution to spin relaxation by means of a power-law expression
$(T_{1})^{-1}=AT^{5}$jarmola2012temperature . Such expression was derived over
60 years ago by assuming the term $\hat{\hat{\mathbf{R4}}}^{2-\mathrm{ph}}$ as
the dominating driving mechanism for Raman relaxation, a Debye-like structure
of the phonons’ density of states, and the absence of low-energy spin excited
stateswalker19685 . Our simulations show that contrary to this model’s
assumptions, Raman relaxation does instead proceed through the term
$\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ and that the phonons relevant to
spin relaxation are not at all well described by a Debye model. The analysis
of our simulations indeed shows that a strict $(T_{1})^{-1}\propto T^{5}$
relation does not hold and that the $T$-dependence of spin relaxation has a
more complex profile. As suggested by the $T$-dependence of the function
$G^{2-\mathrm{ph}}$, we successfully modelled the simulations of NV- by means
of the expression

$(T_{1})^{-1}=A\frac{\>e^{-B/k_{B}T}}{(e^{-B/k_{B}T}-1)^{2}}+C\frac{\>e^{-D/k_{B}T}}{(e^{-D/k_{B}T}-1)^{2}}\>,$
(9)

where the fitted coefficients $C=326$ cm-1 and $D=576$ cm-1 nicely correspond
to the first few peaks in the spin-phonon coupling density $V(\omega)$
reported in Fig. 2, in qualitative agreement with a recent
workcambria2022temperature . The same type of expression would apply to
V${}_{B}^{-}$, but differently from NV-, the spin relaxation of the
V${}_{B}^{-}$ is determined in the high-$T$ regime ($k_{B}T>\hbar\omega$) with
respect to the first available modes, thus leading to a simple
$(T_{1})^{-1}=A\>T^{2}$ contribution, also in nice agreement with the
experimentally determined $T^{2.5}$ power-lawgottscholl2021room . To
corroborate the interpretation that the first available phonons are indeed the
prominent ones for spin relaxation, we simulate $T_{1}$ by imposing a variable
cut-off on the energy of phonons included in the computation of
$\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ and show that even the inclusion of
the first few available phonons is enough to obtain the correct order of
magnitude of relaxation rate, with the inclusion of phonons up to 800 cm-1 and
80 cm-1 to entirely recover the converged value of $T_{1}$ for NV- and
$V_{\mathrm{B}}^{-}$, respectively (see Fig. S7). This is in agreement with
previous findings for magnetic moleculeslunghi2022toward  and provides a
simple rationale for explaining the striking difference among the relaxation
rate of NV- and V${}_{\mathrm{B}}^{-}$. Indeed, despite being both materials
with stiff chemical bonds, the 2D nature of h-BN enables low-energy out-of-
plane flexural vibrations that are sensibly populated even at low-$T$.

$\bar{n}_{\alpha}$$\bar{n}_{\alpha}$ \+ 1 $\bar{n}_{\beta}$$\bar{n}_{\beta}$
\+ 1 $\hat{V}_{sph}$$\hat{V}_{ss}$$\gamma_{\alpha}$$\gamma_{\beta}$

Figure 4: Pure dephasing mechanism and correlations between spin-phonon
relaxation time (T1) and coherence time (T2). Top panel: schematic
representation of two-phonon pure dephasing mechanism. Two degenerate phonons
exchange population and lead to dephasing through the spin-phonon coupling
interaction $\hat{V}_{sph}$, similarly to the pure dephasing due to nuclear
spin flip-flop processes among spins with identical gyromagnetic factors
$\gamma_{\alpha}=\gamma_{\beta}$. Middle panel: Black and red symbols
represent the simulated and experimentalbar2013solid  results for NV-. The
dashed line reports $T_{2}=0.56T_{1}$. Bottom panel: Black symbols represent
the simulated results for V${}_{\mathrm{B}}^{-}$. The dashed line reports
$T_{2}=0.1T_{1}$.

Spin-Phonon Decoherence. Finally, we turn to the study of the limits imposed
to spin coherence by spin-phonon relaxation by studying the dynamics of the
out-of-diagonal elements of $\hat{\rho}_{s}$. Considering the effect of the
super-operator $\hat{\hat{\mathbf{R2}}}^{2-\mathrm{ph}}$ we obtain the
analytical result

$(T_{2})^{-1}=(2T_{1})^{-1}+(T_{2}^{*})^{-1}\>,$ (10)

where the first term on the r.h.s. is the commonly reported upper bound of
coherence due to a finite spin lifetime. The term $T_{2}^{*}$ represents the
pure-dephasing term and reads

$\displaystyle\frac{1}{T_{2}^{*}}\propto$
$\displaystyle\sum_{\alpha,\beta}\left(-2V^{\alpha\beta}_{aa}V^{\alpha\beta}_{bb}+V^{\alpha\beta}_{aa}V^{\alpha,\beta}_{aa}+V^{\alpha\beta}_{bb}V^{\alpha\beta}_{bb}\right)\cdot$
(11) $\displaystyle\cdot
G^{2-\mathrm{ph}}(0,\omega_{\alpha},\omega_{\beta})\>,$ (12)

where $a,b=\\{|0\rangle,|1\rangle\\}$. Such a mechanism thus involves an
energy-conserving simultaneous emission and absorption of two degenerate
phonons. If the two phonons have different spin-phonon coupling, this process
leads to decoherence, equivalently to the flip-flop nuclear spin transitions
limiting the coherence of a central electronic spinyao2006theory ;
witzel2006quantum ; seo2016quantum . These two processes are schematically
shown in the top panel of Fig. 4 Notably, the corresponding expression of pure
dephasing for one-phonon processes leads to a null contribution, as energy
conservation would require zero-frequency phonons. Recently, it has been shown
that pure spin dephasing is an active decoherence mechanism in molecular
$S=1/2$ spin qubitsgarlatti2022critical . Here we test the possibility for
this process to also be influencing the values of $T_{2}$ for NV- observed in
dynamical decoupling experimentsbar2013solid . In such experimental
conditions, the contribution of spin-spin dipolar interactions to dephasing is
systematically removed by applying several refocusing $\pi$ pulses, revealing
the contribution of phonons to decoherence and allowing us to compare
simulations and experimental results. We study the evolution of the state
$|\Psi\rangle=(|0\rangle+|1\rangle)/\sqrt{2}$, as generated in experiments
through a $\pi/2$ pulse resonant with the transition
$|0\rangle\rightarrow|1\rangle$bar2013solid . The middle panel of Fig. 4
reports the ratio between the values of $T_{2}$ and $T_{1}$ at different
values of temperature for NV-. The comparison between simulations and
experimental results shows a remarkable agreement between the two, with
simulations correctly capturing the experimentally observed trend $T_{2}\sim
0.5T_{1}$, demonstrating the origin of such relation in the presence of pure
dephasing. Interestingly, the simulations for V${}_{\mathrm{B}}^{-}$ reported
in the bottom panel of Fig. 4 show an even more drastic reduction of $T_{2}$,
namely $T_{2}\sim 0.1T_{1}$.

## Discussion

The theory of spin relaxation in solid-state paramagnetic defects or
impurities has a long-standing tradition with roots in the early days of
magnetic resonance. Seminal contributions to the theory of spin
relaxationvan1940paramagnetic  were rapidly followed by a proliferation of
phenomenological models based on varied assumptions about the details of spin-
phonon coupling and materials’ vibrational density of
statesshrivastava1983theory ; lunghi2022spin . Until very recently, such
models have remained the only tool available for the interpretation of
experimental relaxation data. However, experimental data has proved to be far
from easy to interpret and multiple models can often be equally good
explanations of experimental evidence. The field of paramagnetic defects is no
exception and no clear understanding of spin relaxation in this class of
materials has yet been established. Providing a clear interpretation of the
mechanism of spin relaxation is one of the main results of this work. Ab
initio simulations, free from any adjustable parameter, enabled the
identification of two-phonon Raman relaxation due to the quadratic coupling of
low-energy optical phonons and zero-field splitting as the main relaxation
pathway for NV- and V${}_{\mathrm{B}}^{-}$. Although alternative relaxation
pathways might be operative in spin qubits with different spin multiplicities
or low-lying excited electronic stateslunghi2022toward , these results are
likely exportable to all spin qubits exhibiting zero-field splitting, as
suggested by their validity for defects in both 3D and 2D materials. As a
consequence of the correct reinterpretation of the spin relaxation mechanism,
we were also able to unravel the role of phonons on spin decoherence. The
latter is invariably found to be bound from above by $T_{1}$ bar2013solid ;
simin2017locking  instead of the way around, as suggested by one-phonon
relaxation theories. As shown, two-phonon processes introduce a pure dephasing
mechanism and naturally explain such apparent deviations, leading to a
coherent picture of spin dynamics.

The present results have far-reaching consequences and suggest several
possible pathways for future development. For instance, the relevance of low-
energy phonons for the coherence dynamics of $S=1$ spin qubits has several
important consequences for the design of novel spin qubits. Although bulk
diamond probably already offers the best possible environment possible in
terms of the (high) energy of optical vibrations, new spin qubits can be
searched among materials with high-energy optical vibrations and a natural low
density of nuclear spinskanai2022generalized  in order to simultaneously
minimize the effect of both spin-phonon coupling and spin-spin dipolar
interactions.

Moreover, our results suggest that low-energy vibrational states of shallow
defects and nano-diamonds might play a role in reducing the coherence time of
electron spins in these substrates naydenov2011dynamical . Indeed, it has
recently been shown that the drastically reduced $T_{1}$ observed for NV-
nano-diamonds has a temperature profile consistent with phonon-induced
relaxationde2020temperature . DeGuillebon et al. invoked the presence of
magnetic impurities at the surface of the nano-diamonds to interpret their
resultsde2020temperature , but our simulations suggest that surface
vibrational states, qualitatively similar to those of V${}_{\mathrm{B}}^{-}$,
might be enough to explain the drastic reduction in $T_{1}$. 2D materials
hetero-structures and tailored surface coating might also offer a rich
playground to tune the vibrational properties of spin qubits’ host
environment.

Overall, our simulations open a window on the fine details of spin-phonon
coupling in solid-state paramagnetic defects that can eventually inform
several applications. Indeed, the coupling of vibrations and electronic states
are responsible for the temperature-dependent shift in spin energy transitions
often used as a sensing strategy of the spin environmentacosta2010temperature
, as well as for the optical polarization of defects’ spin
statesthiering2018theory . Moreover, whilst minimizing the effect of spin-
phonon coupling is a key future challenge to increase coherence time at high-
temperature, a tailored strong coupling with surface or bulk acoustic waves
holds great promise towards the control and coupling of multiple spin systems
in solid-state architecturesalbrecht2013coupling ; golter2016optomechanical ;
whiteley2019spin . The proposed computational framework has a strong potential
to aid the development of all these areas.

In conclusion, we have here used ab initio spin dynamics simulations to
determine the spin-phonon coupling coefficients in two prototypical solid-
state qubits and demonstrated how this interaction leads to spin relaxation
and decoherence. By using a fully parameter-free computational approach we
have provided novel insights into the nature of two-phonon Raman relaxation
and revealed low-energy phonons as the main source of decoherence in solid-
state qubits at high temperatures. We anticipate that our results will
facilitate the interpretation of experimental relaxation data in spin qubits
and speed up the identification of novel quantum sensors.

## Materials and Methods

Crystal Structures. The defects’ structure has been generated starting from a
supercell of the pristine materials, diamond and h-BN, respectively. Cell
optimizations and simulation of phonons for both the defects have been
performed for two sizes of the supercell, i.e 3x3x3 and 4x4x4 diamond
supercells for NV- and of 9x9 and 12x12 h-BN supercells for
V${}_{\mathrm{B}}^{-}$. Supercells of 3x3x3 and 4x4x4 for NV- contain 215
atoms and 511 atoms, respectively, whereas, 9x9 and 12x12 supercells for
V${}_{\mathrm{B}}^{-}$ contain 161 atoms and 287 atoms. The supercells are
shown in Fig. S1.

Electronic structure simulations. Cell optimizations and phonon calculations
for all the supercells of NV- defect in diamond and V${}_{\mathrm{B}}^{-}$
defect in h-BN were performed within the framework of periodic density
functional theory using CP2K.kuhne2020cp2k  All the optimization and phonon
calculations were done at $\Gamma$-point. A very tight force convergence
criteria of 10−7 a.u. and SCF convergence criteria of 10-10 a.u. for the
energy were employed for the cell optimizations. A plane wave cutoff of 1000
Ry with DZVP-MOLOPT Gaussian basis sets and Goedecker-Tetter Hutter
pseudopotentialsgoedecker1996separable  were used for both systems. The
Perdew-Burke-Ernzerhof (PBE) functional and DFT-D3 methods for dispersion
corrections had been used.perdew1996generalized ; grimme2010consistent

Density Functional Theory, as implemented in the software ORCAneese2020orca ,
had been used to compute the anisotropy tensor $\mathbf{D}$. Clusters
(passivated with H) created from atoms around the defect had been used for
this purpose. A convergence test for different sizes of clusters had been
performed to select the smallest cluster size which can reproduce the
experimental $D$ value, see Table S1. A similar strategy had been taken for
the choice of the basis set and the DFT functional. Convergence tests had been
performed to select the best ones in comparison to the experiments. A detailed
comparison of different basis sets and different functionals has been
presented in the table below, see Table. S2. Based upon the comparison, we
have used 5 Å cluster and BHANDHLY becke1993new  functional for NV- defect and
7.5 Å cluster and PBE0 adamo1999toward  functional for V${}_{\mathrm{B}}^{-}$
defect, respectively. The basis set of DKH-def2-SVP had been used for all the
atoms present in both the defects NV- and V${}_{\mathrm{B}}^{-}$. The
hydrogenated clusters of NV- defect (5 Å) and V${}_{\mathrm{B}}^{-}$ (7.5 Å)
defect have been shown in Fig. S1.

Training of Neural Network and Test set generation. The training sets for the
neural network have been generated by calculating $D$ for over 2000 distorted
structures. The distorted structures were obtained by displacing all the
Cartesian coordinates of all the atoms from the optimized geometry. The
displacement was given by a random value in the range from $-$0.05 Å to 0.05
Å. Among the total of 2000 distorted geometries, 1600 were used for training
the ML models, while 400 were used for validation. The neural network we used
has 3 hidden layers consisting of 128, 64, and 16 nodes, respectively. The
input layer has N nodes, where N is the number of atoms present in the
cluster, and the output layer has 1 node. Each hidden layer has a sigmoid as
the activation function and L2-regularization.

Calculation of spin-phonon coupling coefficients. All the relaxation rates
presented in the main manuscript have been obtained using the spin-phonon
coupling coefficients calculated from the machine learning model which was
trained with $D$ values obtained from DFT calculations. The numerical
differentiation method has been used to calculate the spin-phonon coupling
coefficients. In the numerical differentiation method, one-dimensional grids
of 6 points and two-dimensional grids of 36 points have been used to determine
the first-order and second-order coefficients of spin-phonon coupling. Machine
learning models have been used to calculate the $D$ values for all the points.
Each one of these 2-dimensional grids of points was interpolated using a
third-order polynomial,
$P(x,y)=a_{30}x^{3}+a_{03}y^{3}+a_{21}x^{2}y+a_{12}xy^{2}+a_{20}x^{2}+a_{02}y^{2}+a_{11}xy+a_{10}x+a_{01}y+a_{00}.$
The coefficients $a_{01}$ and $a_{10}$ correspond to first-order spin-phonon
coupling coefficients, $(\partial\hat{H}_{s}/\partial X_{i})$, and coefficient
corresponding to $a_{11}$, $a_{20}$ and $a_{02}$ represents the mixed and pure
second-order Cartesian derivatives $(\partial^{2}\hat{H}_{s}/\partial
X_{i}\partial X_{j})$, where the indexes $i,j$ runs over all the 3N degrees of
freedom of the cell. Once these parameters have been computed, the spin-phonon
coupling coefficients are obtained by means of the relations

$\left(\frac{\partial\hat{H}_{s}}{\partial
q_{{\alpha}}}\right)=\sum_{i}^{3N}\sqrt{\frac{\hbar}{\omega_{\alpha}m_{i}}}L_{i\alpha}\left(\frac{\partial\hat{H}_{s}}{\partial
X_{i}}\right)\>,$ (13)

and

$\left(\frac{\partial^{2}\hat{H}_{s}}{\partial q_{{\alpha}}\partial
q_{{\beta}}}\right)=\sum_{ij}^{3N}\sqrt{\frac{\hbar}{\omega_{\alpha}m_{i}}}\sqrt{\frac{\hbar}{\omega_{\beta}m_{j}}}L_{i\alpha}L_{j\beta}\left(\frac{\partial^{2}\hat{H}_{s}}{\partial
X_{i}\partial X_{j}}\right)\>,$ (14)

where $L_{i\alpha}$ are the eigenvectors of the Hessian matrix.

Supplementary Materials
The full expression of the spin-phonon transition rates. Tables S1-S2. Figures
S1-S10.

Acknowledgements and Funding
This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. [948493]). Computational resources were provided by the
Trinity College Research IT and the Irish Centre for High-End Computing
(ICHEC).

Authors Contributions
A.L. conceived the project. S.M. performed all simulations and analyzed the
data. All authors contributed to the discussion of the results and to writing
the manuscript.

Data Availability
All the spin dynamics simulations discussed in the present manuscript have
been carried out with the software MolForge v1.0.0lunghi2022toward ; MolForge
. The latest version of MolForge is available at
github.com/LunghiGroup/MolForge. All data needed to evaluate the conclusions
in the paper are present in the paper and/or the Supplementary Materials.
Additional data related to this paper may be requested from the authors.

Conflict of interests
The authors declare no competing interests.

## References

  * (1) Alkahtani, M. H. _et al._ Fluorescent nanodiamonds: past, present, and future.  _Nanophotonics_ 7, 1423–1453 (2018).
  * (2) Castelletto, S. & Boretti, A.  Silicon carbide color centers for quantum applications.  _Journal of Physics: Photonics_ 2, 022001 (2020).
  * (3) Balasubramanian, G. _et al._ Ultralong spin coherence time in isotopically engineered diamond.  _Nature materials_ 8, 383–387 (2009).
  * (4) Bar-Gill, N., Pham, L. M., Jarmola, A., Budker, D. & Walsworth, R. L.  Solid-state electronic spin coherence time approaching one second.  _Nature communications_ 4, 1–6 (2013).
  * (5) Schirhagl, R., Chang, K., Loretz, M. & Degen, C. L.  Nitrogen-vacancy centers in diamond: nanoscale sensors for physics and biology.  _Annual review of physical chemistry_ 65, 83–105 (2014).
  * (6) Gruber, A. _et al._ Scanning confocal optical microscopy and magnetic resonance on single defect centers.  _Science_ 276, 2012–2014 (1997).
  * (7) Jelezko, F., Gaebel, T., Popa, I., Gruber, A. & Wrachtrup, J.  Observation of coherent oscillations in a single electron spin.  _Physical review letters_ 92, 076401 (2004).
  * (8) Grotz, B. _et al._ Sensing external spins with nitrogen-vacancy diamond.  _New Journal of Physics_ 13, 055004 (2011).
  * (9) Tetienne, J.-P. _et al._ Spin relaxometry of single nitrogen-vacancy defects in diamond nanocrystals for magnetic noise sensing.  _Physical Review B_ 87, 235436 (2013).
  * (10) Togan, E. _et al._ Quantum entanglement between an optical photon and a solid-state spin qubit.  _Nature_ 466, 730–734 (2010).
  * (11) Wolfowicz, G. _et al._ Quantum guidelines for solid-state spin defects.  _Nature Reviews Materials_ 6, 906–925 (2021).
  * (12) Gottscholl, A. _et al._ Initialization and read-out of intrinsic spin defects in a van der waals crystal at room temperature.  _Nature materials_ 19, 540–545 (2020).
  * (13) Gottscholl, A. _et al._ Room temperature coherent control of spin defects in hexagonal boron nitride.  _Science Advances_ 7, eabf3630 (2021).
  * (14) Stern, H. L. _et al._ Room-temperature optically detected magnetic resonance of single defects in hexagonal boron nitride.  _Nature communications_ 13, 618 (2022).
  * (15) Walker, M.  At 5 spin–lattice relaxation rate for non-kramers ions.  _Canadian Journal of Physics_ 46, 1347–1353 (1968).
  * (16) Jarmola, A., Acosta, V., Jensen, K., Chemerisov, S. & Budker, D.  Temperature-and magnetic-field-dependent longitudinal spin relaxation in nitrogen-vacancy ensembles in diamond.  _Physical review letters_ 108, 197601 (2012).
  * (17) Degen, C. L., Reinhard, F. & Cappellaro, P.  Quantum sensing.  _Reviews of modern physics_ 89, 035002 (2017).
  * (18) Lunghi, A. & Sanvito, S.  How do phonons relax molecular spins?  _Science advances_ 5, eaax7163 (2019).
  * (19) Escalera-Moreno, L., Baldoví, J. J., Gaita-Ariño, A. & Coronado, E.  Exploring the high-temperature frontier in molecular nanomagnets: from lanthanides to actinides.  _Inorganic Chemistry_ 58, 11883–11892 (2019).
  * (20) Lunghi, A. & Sanvito, S.  The limit of spin lifetime in solid-state electronic spins.  _The Journal of Physical Chemistry Letters_ 11, 6273–6278 (2020).
  * (21) Reta, D., Kragskow, J. G. & Chilton, N. F.  Ab initio prediction of high-temperature magnetic relaxation rates in single-molecule magnets.  _Journal of the American Chemical Society_ 143, 5943–5950 (2021).
  * (22) Lunghi, A.  Toward exact predictions of spin-phonon relaxation times: An ab initio implementation of open quantum systems theory.  _Science Advances_ 8, eabn7880 (2022).
  * (23) Mondal, S. & Lunghi, A.  Unraveling the contributions to spin–lattice relaxation in kramers single-molecule magnets.  _Journal of the American Chemical Society_ (2022).
  * (24) Astner, T. _et al._ Solid-state electron spin lifetime limited by phononic vacuum modes.  _Nature materials_ 17, 313–317 (2018).
  * (25) Tang, H., Barr, A. R., Wang, G., Cappellaro, P. & Li, J.  First-principles calculation of the temperature-dependent transition energies in spin defects.  _arXiv preprint arXiv:2205.02791_ (2022).
  * (26) Cambria, M. _et al._ Temperature-dependent phonon-induced relaxation of the nitrogen-vacancy spin triplet in diamond.  _arXiv preprint arXiv:2209.14446_ (2022).
  * (27) Lunghi, A.  Spin-phonon relaxation in magnetic molecules: Theory, predictions and insights.  _arXiv preprint arXiv:2202.03776_ (2022).
  * (28) Bayliss, S. _et al._ Optically addressable molecular spins for quantum information processing.  _Science_ 370, 1309–1312 (2020).
  * (29) Garlatti, E. _et al._ The critical role of ultra-low energy vibrations in the relaxation dynamics of molecular qubits.  _Nature Communications_ In Press (2023).
  * (30) Garlatti, E. _et al._ Unveiling phonons in a molecular qubit with four-dimensional inelastic neutron scattering and density functional theory.  _Nature communications_ 11, 1–10 (2020).
  * (31) Gugler, J. _et al._ Ab initio calculation of the spin lattice relaxation time t 1 for nitrogen-vacancy centers in diamond.  _Physical Review B_ 98, 214442 (2018).
  * (32) Lunghi, A.  Molforge v.1.0.0 DOI:10.5281/zenodo.7596042 (2023).
  * (33) Yao, W., Liu, R.-B. & Sham, L.  Theory of electron spin decoherence by interacting nuclear spins in a quantum dot.  _Physical Review B_ 74, 195301 (2006).
  * (34) Witzel, W. & Sarma, S. D.  Quantum theory for electron spin decoherence induced by nuclear spin dynamics in semiconductor quantum computer architectures: Spectral diffusion of localized electron spins in the nuclear solid-state environment.  _Physical Review B_ 74, 035322 (2006).
  * (35) Seo, H. _et al._ Quantum decoherence dynamics of divacancy spins in silicon carbide.  _Nature communications_ 7, 1–9 (2016).
  * (36) Van Vleck, J.  Paramagnetic relaxation times for titanium and chrome alum.  _Physical Review_ 57, 426 (1940).
  * (37) Shrivastava, K.  Theory of spin–lattice relaxation.  _physica status solidi (b)_ 117, 437–458 (1983).
  * (38) Simin, D. _et al._ Locking of electron spin coherence above 20 ms in natural silicon carbide.  _Physical Review B_ 95, 161201 (2017).
  * (39) Kanai, S. _et al._ Generalized scaling of spin qubit coherence in over 12,000 host materials.  _Proceedings of the National Academy of Sciences_ 119, e2121808119 (2022).
  * (40) Naydenov, B. _et al._ Dynamical decoupling of a single-electron spin at room temperature.  _Physical Review B_ 83, 081201 (2011).
  * (41) de Guillebon, T., Vindolet, B., Roch, J.-F., Jacques, V. & Rondin, L.  Temperature dependence of the longitudinal spin relaxation time t 1 of single nitrogen-vacancy centers in nanodiamonds.  _Physical Review B_ 102, 165427 (2020).
  * (42) Acosta, V. M. _et al._ Temperature dependence of the nitrogen-vacancy magnetic resonance in diamond.  _Physical review letters_ 104, 070801 (2010).
  * (43) Thiering, G. & Gali, A.  Theory of the optical spin-polarization loop of the nitrogen-vacancy center in diamond.  _Physical Review B_ 98, 085207 (2018).
  * (44) Albrecht, A., Retzker, A., Jelezko, F. & Plenio, M. B.  Coupling of nitrogen vacancy centres in nanodiamonds by means of phonons.  _New Journal of Physics_ 15, 083014 (2013).
  * (45) Golter, D. A., Oo, T., Amezcua, M., Stewart, K. A. & Wang, H.  Optomechanical quantum control of a nitrogen-vacancy center in diamond.  _Physical review letters_ 116, 143602 (2016).
  * (46) Whiteley, S. J. _et al._ Spin–phonon interactions in silicon carbide addressed by gaussian acoustics.  _Nature Physics_ 15, 490–495 (2019).
  * (47) Kühne, T. D. _et al._ Cp2k: An electronic structure and molecular dynamics software package-quickstep: Efficient and accurate electronic structure calculations.  _The Journal of Chemical Physics_ 152, 194103 (2020).
  * (48) Goedecker, S., Teter, M. & Hutter, J.  Separable dual-space gaussian pseudopotentials.  _Physical Review B_ 54, 1703 (1996).
  * (49) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Physical review letters_ 77, 3865 (1996).
  * (50) Grimme, S., Antony, J., Ehrlich, S. & Krieg, H.  A consistent and accurate ab initio parametrization of density functional dispersion correction (dft-d) for the 94 elements h-pu.  _The Journal of chemical physics_ 132, 154104 (2010).
  * (51) Neese, F., Wennmohs, F., Becker, U. & Riplinger, C.  The orca quantum chemistry program package.  _The Journal of chemical physics_ 152, 224108 (2020).
  * (52) Becke, A. D.  A new mixing of hartree–fock and local density-functional theories.  _The Journal of chemical physics_ 98, 1372–1377 (1993).
  * (53) Adamo, C. & Barone, V.  Toward reliable density functional methods without adjustable parameters: The pbe0 model.  _The Journal of chemical physics_ 110, 6158–6170 (1999).