* Powell, (1978) Powell, M. J. (1978).  Algorithms for nonlinear constraints that use lagrangian functions.  Mathematical programming, 14:224–248.
  * Præstgaard and Wellner, (1993) Præstgaard, J. and Wellner, J. A. (1993).  Exchangeably weighted bootstraps of the general empirical process.  The Annals of Probability, pages 2053–2086.
  * Quintos, (1998) Quintos, C. E. (1998).  Analysis of cointegration vectors using the gmm approach.  Journal of econometrics, 85(1):155–188.
  * Rothenberg, (1971) Rothenberg, T. J. (1971).  Identification in parametric models.  Econometrica: Journal of the Econometric Society, pages 577–591.
  * Rothenberg and Stock, (1997) Rothenberg, T. J. and Stock, J. H. (1997).  Inference in a nearly integrated autoregressive model with nonnormal innovations.  Journal of Econometrics, 80(2):269–286.
  * Rubinov et al., (2002) Rubinov, A. M., Huang, X., and Yang, X. (2002).  The zero duality gap property and lower semicontinuity of the perturbation function.  Mathematics of Operations Research, 27(4):775–791.
  * Ruszczyński and Shapiro, (2006) Ruszczyński, A. and Shapiro, A. (2006).  Conditional risk mappings.  Mathematics of operations research, 31(3):544–561.
  * Saikkonen, (1991) Saikkonen, P. (1991).  Asymptotically efficient estimation of cointegration regressions.  Econometric theory, 7(1):1–21.
  * Sargan, (1976) Sargan, J. D. (1976).  Econometric estimators and the edgeworth approximation.  Econometrica: Journal of the Econometric Society, pages 421–448.
  * Seo and Shin, (2016) Seo, M. H. and Shin, Y. (2016).  Dynamic panels with threshold effect and endogeneity.  Journal of econometrics, 195(2):169–186.
  * Spini, (2021) Spini, P. E. (2021).  Robustness, heterogeneous treatment effects and covariate shifts.  arXiv preprint arXiv:2112.09259.
  * Stock and Watson, (1993) Stock, J. H. and Watson, M. W. (1993).  A simple estimator of cointegrating vectors in higher order integrated systems.  Econometrica: journal of the Econometric Society, pages 783–820.
  * Su and Chen, (2013) Su, L. and Chen, Q. (2013).  Testing homogeneity in panel data models with interactive fixed effects.  Econometric Theory, 29(6):1079–1135.
  * Su et al., (2015) Su, L., Jin, S., and Zhang, Y. (2015).  Specification test for panel data models with interactive fixed effects.  Journal of Econometrics, 186(1):222–244.
  * Su and Lu, (2013) Su, L. and Lu, X. (2013).  Nonparametric dynamic panel data models: kernel estimation and specification testing.  Journal of Econometrics, 176(2):112–133.
  * Su et al., (2016) Su, L., Shi, Z., and Phillips, P. C. (2016).  Identifying latent structures in panel data.  Econometrica, 84(6):2215–2264.
  * Tabri, (2021) Tabri, R. V. (2021).  The information projection in moment inequality models: Existence, dual representation, and approximation.  arXiv preprint arXiv:2107.07140.
  * Van Garderen, (1997) Van Garderen, K. J. (1997).  Curved exponential models in econometrics.  Econometric Theory, 13(6):771–790.
  * Verdier, (2016) Verdier, V. (2016).  Estimation of dynamic panel data models with cross-sectional dependence: Using cluster dependence for efficiency.  Journal of Applied Econometrics, 31(1):85–105.
  * Vesely et al., (2023) Vesely, A., Finos, L., and Goeman, J. J. (2023).  Permutation-based true discovery guarantee by sum tests.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 85(3):664–683.
  * Wagner et al., (2020) Wagner, M., Grabarczyk, P., and Hong, S. H. (2020).  Fully modified ols estimation and inference for seemingly unrelated cointegrating polynomial regressions and the environmental kuznets curve for carbon dioxide emissions.  Journal of Econometrics, 214(1):216–255.
  * Wald, (1950) Wald, A. (1950).  Statistical decision functions.
  * Westerlund, (2019) Westerlund, J. (2019).  On estimation and inference in heterogeneous panel regressions with interactive effects.  Journal of Time Series Analysis, 40(5):852–857.
  * Westerlund et al., (2022) Westerlund, J., Norkutė, M., and Stauskas, O. (2022).  The factor analytical approach in trending near unit root panels.  Journal of Time Series Analysis, 43(3):501–508.
  * White, (1996) White, H. (1996).  Estimation, inference and specification analysis.  Number 22. Cambridge university press.
  * White, (2014) White, H. (2014).  Asymptotic theory for econometricians.  Academic press.
  * Wintoki et al., (2012) Wintoki, M. B., Linck, J. S., and Netter, J. M. (2012).  Endogeneity and the dynamics of internal corporate governance.  Journal of financial economics, 105(3):581–606.
  * Wu et al., (2023) Wu, J., Song, X., and Xiao, Z. (2023).  Testing for trend specifications in panel data models.  Journal of Business & Economic Statistics, 41(2):453–466.
  * Yalcin and Kasimbeyli, (2020) Yalcin, G. D. and Kasimbeyli, R. (2020).  On weak conjugacy, augmented lagrangians and duality in nonconvex optimization.  Mathematical Methods of Operations Research, 92:199–228.
  * Yan and Cheng, (2022) Yan, Y. and Cheng, T. (2022).  Factor-augmented forecasting regressions with threshold effects.  The Econometrics Journal, 25(1):134–154.
  * Yu and Fan, (2021) Yu, P. and Fan, X. (2021).  Threshold regression with a threshold boundary.  Journal of Business & Economic Statistics, 39(4):953–971.
  * Zellner et al., (1974) Zellner, A., Palm, F., et al. (1974).  Time series analysis and simultaneous equation econometric models.  Journal of econometrics, 2(1):17–54.


# Dictionary Learning for Clustering on Hyperspectral Images

Joshua Bruton  School of Computer Science and
Applied Mathematics
University of the Witwatersrand
Johannesburg, Gauteng  Hairong Wang  School of Computer Science and
Applied Mathematics
University of the Witwatersrand
Johannesburg, Gauteng

###### Abstract

Dictionary learning and sparse coding have been widely studied as mechanisms
for unsupervised feature learning. Unsupervised learning could bring enormous
benefit to the processing of hyperspectral images and to other remote sensing
data analysis because labelled data is often scarce in this field. We propose
a method for clustering the pixels of hyperspectral images using sparse
coefficients computed from a representative dictionary as features. We show
empirically that the proposed method works more effectively than clustering on
the original pixels. We also demonstrate that our approach, in certain
circumstances, outperforms the clustering results of features extracted using
principal component analysis and non-negative matrix factorisation.
Furthermore, our method is suitable for applications in repetitively
clustering an ever-growing amount of high-dimensional data, which is the case
when working with hyperspectral satellite imagery.

###### Index Terms:

dictionary learning, hyperspectral images, signal processing, spectral
clustering, unsupervised learning, feature extraction

##  I Introduction

When clustering high-dimensional data we often seek an effective
representation in a latent space in the pursuit of more discriminative
features. Dictionary learning provides a way to learn this representation of
the latent structure of the data. After a dictionary has been acquired, sparse
coding allows one to reconstruct the original data using a linear combination
of the dictionary atoms with a set of weight vectors called sparse
coefficients. A sparse coefficient vector is higher-dimensional than the
original samples in the data however they are sparse, i.e. most of the
components have zero values. To accurately reconstruct an original datum, a
sparse coefficient must prioritise specific atoms (columns) in the dictionary
in the correct proportion. Given that the dictionary provides a general view
of the data in a lower-dimensional space, the atoms that sparse coding selects
to represent a particular input must provide for any variance within the mean
structure of the data. Thus, as has been shown in recent literature, a sparse
coefficient matrix computed using a representative dictionary learned from a
specific data set will provide useful discriminative information [12, 11, 9,
6].

In this paper, we propose a dictionary-based clustering method that takes
advantage of that discriminative information. There exists a trend across many
fields towards favouring unsupervised and semi-supervised techniques. The
amount of data being generated in the modern world cannot be managed solely
through the paradigm of supervised learning, in which a label has to be
assigned to every sample. Other techniques which require less a-priori
knowledge of the data are needed. Clustering, or automatically assigning
labels to samples in a data set based on some similarity criteria, is one such
technique. The fields of signal processing and remote sensing are no exception
in this trend towards less supervised approaches. We chose to apply our method
to hyperspectral images (HSIs), in particular three widely used hyperspectral
data sets. Modern hyperspectral sensors have high spectral resolution covering
up to several hundred bands. Hyperspectral images acquired using such sensors
are high dimensional not only in the spatial domain, but also in the spectral
domain. This justifies the use of feature extraction techniques that can
represent this high-dimensional information in a lower-dimensional space.

Before outlining the method in detail, we summarise some of the related work
in the field and explain some important background information in section II.
We provide insight into parameter-tuning and implementation details in section
III. In section IV, we demonstrate the efficacy of the method empirically and
show that the method outperforms two simple approaches to feature extraction,
principal component analysis (PCA) and non-negative matrix factorisation
(NMF), in certain circumstances. We end the paper in section V with concluding
remarks and describe some of the future work that could add to the results
shown here.

##  II Background and Related Work

###  II-A Online Dictionary Learning and Sparse Coding

In this subsection we aim to provide a working understanding of dictionary
learning and of the algorithm we use to perform it. In the context of this
paper, dictionaries are dense matrices of base elements that best describe the
latent structure of the data from which they are learned. Dictionary learning
is a ubiquitous technique with which to learn these data representations. More
formally, when creating a dictionary, we aim to represent the input data in a
matrix X with two matrices: a dictionary, D, and a matrix of sparse
coefficients, A. Here, $\mathbf{D}\textbf{A}$ is called the sparse
approximation of X. Basis pursuit is not a working solution, as we also
require the dictionary to be overcomplete so as not to be limited by
orthogonality. By “overcomplete” we mean that if
$\mathbf{X}\in\mathbb{R}^{m\times n}$ and $\mathbf{D}\in\mathbb{R}^{m\times
k}$ then $k>m$. The dictionary being overcomplete guarantees that, for
$\textbf{x}\in\textbf{X}$, $\textbf{x}=\mathbf{D}\boldsymbol{\alpha}$ will
have infinitely many solutions for $\boldsymbol{\alpha}$, we evaluate these
solutions based on their sparsity, this is called “sparse coding” or “atom
decomposition”[2]. Methods for sparse coding are usually called “pursuit
algorithms” and they aim to find approximate solutions to the problem because
exact computation is intractable. We perform sparse coding with a greedy
pursuit algorithm called orthogonal matching pursuit (OMP) [4].

In the problem of dictionary learning we are given X and aim to find D and A.
If we frame the problem around the solution for A, we call this process
“simultaneous sparse approximation”. We cannot approach this problem
analytically, and so we turn to optimisation techniques. Optimisation
techniques effectively split simultaneous sparse approximation into two
manageable steps. In the first step, we compute the sparse coefficient matrix
A given the dictionary D; we can do this with any pursuit algorithm. In the
second step we are given the sparse coefficient matrix A and must update D
sensibly. An intuitive equation for performing the dictionary update is
provided by [1]:

$\textbf{D}^{\prime}=\underset{\textbf{D}}{\mathrm{argmin}}{1\over
t}\sum_{i=1}^{t}{{1\over
2}||\textbf{x}_{i}-\textbf{D}\boldsymbol{\alpha}_{i}||_{2}^{2}+\lambda||\boldsymbol{\alpha}_{i}||_{1}}$
(1)

where $\textbf{D}^{\prime}$ is the new dictionary,
$\textbf{x}_{i}\in\textbf{X}$, $\boldsymbol{\alpha}_{i}\in\textbf{A}$ and
$\lambda$ is a regularisation parameter. This equation endeavours to find the
dictionary D that minimises the reconstruction error under some constraints.
The first term, minimised in isolation from the second, would find the closest
approximation of $\textbf{x}_{i}$ without considering the sparsity of
$\boldsymbol{\alpha}_{i}$. We want $\boldsymbol{\alpha}_{i}$ to be sparse, and
so the second term penalises $\boldsymbol{\alpha}_{i}$ in proportion to the
number of nonzero elements it contains[2]. This equation can be altered
algebraically to instead operate on each column of the dictionary separately,
this is done in Algorithm 1. We also have to ensure that the columns of D are
normalised at the end of each update step, in practice this is enforced by
normalising the columns in D as they are updated. The algorithm described here
is called online dictionary learning (ODL) and was originally proposed in [1].
A pseudo-code implementation of ODL can be seen in Algorithm 1.

Algorithm 1 Online dictionary learning algorithm [1]

0: Number of iterations $T$, dictionary $\textbf{D}_{0}$
$\in\mathbb{R}^{m\times k}$, an algorithm, p(x), to draw random variable x
from data set X, regularisation term $\lambda$

1: BEGIN

2: Initialise $\textbf{A}_{0}\leftarrow 0$ and $\textbf{B}_{0}\leftarrow 0$

3: for $t=1$ to $T$:

4: draw $\textbf{x}_{t}$ from data set using $p(\textbf{x})$

5: sparse coding step using OMP to compute $\boldsymbol{\alpha}$

6:
$\textbf{A}_{t}\leftarrow\textbf{A}_{t-1}+\boldsymbol{\alpha}_{t}\boldsymbol{\alpha}_{t}^{T}$

7:
$\textbf{B}_{t}\leftarrow\textbf{B}_{t-1}+\textbf{x}_{t}\boldsymbol{\alpha}_{t}^{T}$

8: for $i=1$ to $k$:

9: Update $i^{\text{th}}$ column of dictionary using (1)

10: END

ODL is distinct from its predecessors, K-SVD[2] and the method of optimal
directions (MOD)[3], in many ways, but the most significant distinction is
that at each iteration of ODL we process a single datum and update the
dictionary accordingly. ODL allows one to routinely update an already learned
dictionary with new data, without reconsidering the data that has already been
seen. This significantly reduces the computational cost of maintaining a
representative dictionary for a regularly updated data set. This quality also
allows ODL to scale more effectively to larger data sets because it allows it
to take advantage of mini-batches.

With the techniques of dictionary learning and sparse coding we can compute
the sparse coefficient matrix A. $\mathbf{D}\textbf{A}$ is a highly compressed
version of the original data, and this approximation can provide high-fidelity
reconstructions. Dictionary learning and sparse coding are often used for
noise-removal, in-painting and signal compression. There has also been much
research conducted into the discriminative power of the sparse coefficient
matrix A. Typically, this experimentation is conducted with supervised
learning tasks in mind. As mentioned earlier, this work has a different focus.

###  II-B Spectral Clustering and KMeans

Clustering is the task of automatically labelling samples such that samples
with the same label are similar by some metric, each group of samples is
called a cluster [7]. The samples defining the centre of each cluster are
called centroids. The KMeans algorithm is a popular approach to clustering.
Each sample in the data is assigned to the centroid it is most similar to.
After all the assignments are complete, the centroids are updated to the
centre (mean) of all the points in their respective clusters. This process is
repeated until there is no change, or minimal change, in the position of each
centroid. In KMeans, similarity is typically defined in terms of Euclidean
distance. The simplicity of KMeans is appealing however the algorithm makes a
number of limiting assumptions about the structure of the data [2].

KMeans defines a circle around each centroid, if a sample falls within that
circle it is considered a member of the cluster that that centroid defines.
This means that KMeans can only define clusters that are both convex and
circular. One method for tempering this assumption is applying expectation
maximisation (EM) to a Gaussian mixture model. This allows for ellipsoidal
clusters as well as circular ones. This technique does not scale as
efficiently to high-dimensional data and so also requires the application of a
feature extraction algorithm like PCA. Another limitation of this technique is
that it does not address the assumption of convexity [15].

Spectral clustering (SC) addresses the problem of convexity by clustering with
respect to a mathematical graph. To achieve this we first construct a graph
from our data, the graph can be constructed using any technique but KNN is
certainly the most common [13]. A graph containing $N$ vertices, corresponding
to $N$ datum $\\{S_{i}\\}_{i=1}^{N}$, can be represented as an $N\times N$
affinity matrix. In this affinity matrix each entry $(i,j)$ represents the
weight of the edge between vertex $i$ and vertex $j$, we refer to each of
these weights as $W_{i,j}$. From an $N\times N$ affinity matrix we can obtain
$N$ elements $d_{i}$ with the equation $d_{i}=\sum_{j=1}^{N}{W_{i,j}}$ for
$i\in\\{1,2,...,N\\}$, these define the diagonal entries of the degree matrix.
We then construct the $N\times N$ Laplacian matrix, L, using the following
rule [14]:

$\textbf{L}_{i,j}=\left\\{\begin{array}[]{lr}d_{i}&:i=j\\\
-W_{i,j}&:\text{otherwise}\end{array}\right.$ (2)

Our goal in constructing the Laplacian matrix is to perform graph
partitioning, the Laplacian matrix has useful properties to facilitate this.
We first decompose the Laplacian matrix into eigenvalues and eigenvectors, for
an $N\times N$ Laplacian matrix we will obtain $N$ eigenvalues,
$\\{\lambda_{1},\lambda_{2},...,\lambda_{N}\\}$, and $N$ corresponding
eigenvectors, $\\{\textbf{x}_{1},\textbf{x}_{2},...,\textbf{x}_{N}\\}$. The
eigenvalues and eigenvectors are sorted in ascending order of eigenvalue, so
$\lambda_{1}<\lambda_{2}<...<\lambda_{N}$, this is called the “spectrum” of
the Laplacian matrix[13]. Eigenvalues appearing earlier in the spectrum
correspond to points of lower connectivity in the affinity graph. Each of
those eigenvalues corresponds to an $N$-dimensional eigenvector, and each of
the $N$ entries in each eigenvector correspond to a datum in our original
data. We cluster the data using these eigenvectors. Let $\textbf{x}_{i,j}$
represent the $j$-th entry in the $i$-th eigenvector. We can then represent
each of our datum $S_{i}$ as
$(\textbf{x}_{1,i},\textbf{x}_{2,i},\textbf{x}_{3,i},...,\textbf{x}_{N,i})$.
When deciding which of the eigenvectors to use in this representation, two
factors must be considered. First, if the original affinity graph was fully-
connected, then the first eigenvector (corresponding to the lowest eigenvalue)
need not be considered as its components are constant. Second, the number of
subsequent eigenvectors (not counting the first) should be one less than the
number of desired clusters[14]. Using this final representation of the data as
input into the KMeans algorithm will provide a better clustering.

###  II-C Joint Sparse Representation

As mentioned in subsection II-A and shown in Algorithm 1, ODL updates columns
in the dictionary iteratively with respect to randomly selected samples. Joint
sparse representation (JSR) alters this update technique to operate over a
sliding window instead of a single pixel, an illustration of this can be seen
in Fig. 1. This requires the dictionary to use specific columns to represent
pixels that are close together in the spatial domain of the image; thus
enforcing the idea of spatial locality.

Altering a standard ODL implementation to include JSR involves constructing a
set of overlapping tiles, $\\{{\textbf{T}_{i}}\\}_{i=1}^{N}$, each tile
centred on one of the pixels in the image. Padding is not necessary, simply
ignoring the borders is preferable given the size of the HSIs. A tile with
dimensions $p$ and $q$, $\textbf{T}_{i}\in\mathbb{R}^{p\times q}$, can be
considered it’s own sparse coding problem, each producing a matrix of sparse
coefficients $\textbf{Q}_{i}\in\mathbb{R}^{p\times q}$. Some pursuit
algorithms have been altered to perform sparse coding on matrices, one such
algorithm is simultaneous orthogonal matching pursuit (SOMP), which is
outlined in [17]. The dictionary update step is also quite similar to (1),
with $\textbf{T}_{i}$ replacing $\textbf{x}_{i}$ and $\textbf{Q}_{i}$
replacing $\boldsymbol{\alpha}_{i}$, and can be seen in [18].

Figure 1: Standard ODL top, JSR bottom

###  II-D Related Work

Much research has been conducted into clustering the pixels of HSIs, although
typically that research is focused on advancing clustering algorithms
themselves, as opposed to using feature extraction techniques to improve
existing approaches. This can be seen in [7] and [8]. The former presents an
algorithm known as k Stochastic Expectation Maximisation, improving on
traditional EM by selecting the number of clusters, k, using a measure of
entropy instead of specifying k as a hyperparameter. The latter work has two
significant contributions, DLSS (spectral-spatial diffusion learning) and
active learning. DLSS is a form of clustering relying on an advanced notion of
distance called diffusion distance. Active learning is essentially semi-
supervised learning, in which points that are difficult for DLSS to identify
can be specified before clustering begins, which results in improved accuracy.
Much of [8] is devoted to expounding on these two concepts and they are out of
the scope of this work.

Noting the above, some works standout as similar to ours, although none of
them deal with unsupervised learning in HSIs. Reference [12] deals
specifically with the use of sparse representation as a preprocessing step to
broader computer vision problems, like image classification and facial
recognition. Furthermore, a noteworthy conclusion drawn in [12] states that
sparsity could provide a powerful prior for inference with high dimensional
visual data that have intricate low-dimensional structure. This statement is
exactly in line with the initial hypothesis that the work in this paper is
based on. Reference [11] is also related to this research area, although that
work differs from ours in some key areas. Reference [11] shows the benefits of
constructing affinity matrices for SC from sparse coefficients. However, [11]
does not consider learning an overcomplete dictionary as a separate problem
prior to this step. One of the key findings of this work deals with the
dictionary learning process specifically, and how adjustments in the
hyperparameters of ODL affect clustering accuracy. This extension allows for
the use of dictionary learning techniques that could further improve the
construction of the affinity matrices in SC, e.g. JSR. The approach to
simultaneous sparse approximation used in [11] is a technique called non-
negative matrix factorisation (NMF). NMF is a powerful dimensionality
reduction technique similar to dictionary learning and has been extended into
a parallel dictionary learning algorithm by [16]. We will compare our method
to NMF in section IV below.

Reference [9] also more explicitly considers dictionary learning. However, in
their clustering framework, a different dictionary is trained for every
cluster in the data. The similarity of a given sample to each cluster is
defined by the quality of the reconstruction each of the dictionaries
provides. They also define a novel way to measure the quality of the
reconstruction, based on the more traditional reconstruction error defined in
(1). The paper also includes the use of SC, specifically for the
initialisation of the dictionaries defining each cluster. Unlike [9], we learn
a unified dictionary for the entire data set as opposed to learning a
different dictionary for each cluster in the data.

##  III Proposed Method

###  III-A Overview of the Proposed Method

Following the detailed background provided in section II, our proposed method
can be summarised in the following three steps:

  1. 1.

Dictionary learning: in this step, given a hyperspectral image,
$\mathbf{X}\in\mathbb{R}^{m\times n}$, we aim to learn a representative
dictionary $\mathbf{D}\in\mathbb{R}^{m\times k}$. This can be achieved using
any dictionary learning technique but we proceed by using ODL, shown in
Algorithm 1.

  2. 2.

Sparse coding: once the dictionary $\mathbf{D}\in\mathbb{R}^{m\times k}$ has
been found, we can proceed to find the sparse coefficient set
$\mathbf{A}\in\mathbb{R}^{k\times n}$ such that
$\mathbf{X}\approx\mathbf{D}\mathbf{A}$; which means for all
$\textbf{x}\in\mathbf{X}$ we have a sparse coefficient
$\boldsymbol{\alpha}\in\mathbf{A}$ such that
$\textbf{x}\approx\mathbf{D}\boldsymbol{\alpha}$.

  3. 3.

Clustering: We now apply a clustering technique to the sparse coefficients
$\textbf{A}\in\mathbb{R}^{k\times n}$ found in step 2. We have had the most
success using SC, we construct the affinity matrix using KNN and then cluster
the resulting eigenvectors using KMeans. We also present some results obtained
by applying KMeans directly to A.

The method can be evaluated by comparing the resultant clustering after step 3
to a basic clustering method, in which the hyperspectral image clustering is
carried out in the original pixel domain instead of on the sparse coefficient
vectors. Comparing cluster accuracy is a non-trivial task, given that the
numerical value representing each label is assigned arbitrarily, a simple
element-wise comparison of a clustering with the ground truth is meaningless.
Although visually assessing the accuracy of a particular clustering is
possible and oftentimes simple, having a more consistent measure is desirable.
Two popular metrics for this task are AMI (Adjusted Mutual Information score)
and ARI (Adjusted Rand Index). According to [5], when comparing clusters on
data with unbalanced or small clusters, one should prioritise the AMI. Given
that HSIs typically have these qualities, we have chosen to use AMI to
evaluate the method.

AMI will return $1.0$ if the clustering is perfect, the average return for
entirely random clusterings is $0$, this means AMI can also be negative. The
procedure for computing the AMI requires both the ground-truth and the
clustering being evaluated, so labelled data is required. AMI is an extension
of an earlier score for measuring the accuracy of clusterings, called mutual
information (MI) score. MI attempts to quantify the similarity between two
partitions, it does this probabilistically. We are given two equally-sized
partitions with the same number of clusters in each, we denote the number of
clusters as k and the total number of datum as $n$. Then, let
$\\{G_{i}\\}_{i=1}^{k}$ and $\\{L_{i}\\}_{i=1}^{k}$ denote the k clusters in
each of our partitions; in our case, the ground-truth and predicted labels,
respectively. We can now find the probability that a point drawn at random
from each partition is associated with a particular cluster, $i$, using
$\textbf{P}(i)={|G_{i}|\over{n}}$ and
$\textbf{P}^{\prime}(i)={|L_{i}|\over{n}}$. We can also compute the
probability that a random point is in cluster $i$ in the ground truth and
cluster $j$ in the prediction by evaluating $\textbf{P}(i,j)={{|G_{i}\cap
L_{j}|}\over{n}}$. With these expressions, we can define the mutual
information of partitions G and L as:

$MI(G,L)=\sum_{i=1}^{k}{\sum_{j=1}^{k}{\textbf{P}(i,j)\log{{\textbf{P}(i,j)\over{\textbf{P}(i)}\textbf{P}^{\prime}(j)}}}}$
(3)

Given this formulation of MI we also know that the entropy of partitions $G$
and $L$ are given by $E(G)=-\sum_{i=1}^{k}{\textbf{P}(i)\log{\textbf{P}(i)}}$
and
$E(L)=-\sum_{i=1}^{k}{\textbf{P}^{\prime}(i)\log{\textbf{P}^{\prime}(i)}}$,
respectively. As per [5], extending (3) and incorporating entropy, the AMI of
partitions G and L is defined as follows:

$AMI(G,L)={{MI(G,L)-\mathbb{E}[MI(G,L)]}\over{\max{[E(G),E(L)]}-\mathbb{E}[MI(G,L)]}}$
(4)

The most important distinction between MI and AMI is that AMI incorporates the
idea of chance. The score that a clustering receives from AMI will be reduced
in proportion to how easily its MI score could be produced by a random
clustering, this is encapsulated in $\mathbb{E}[MI(G,L)]$, the expected MI.
AMI is similar to the popular Cohen’s kappa coefficient (widely used in
supervised learning), both endeavour to incorporate the percentage agreement
and the possibility of that percentage agreement occurring by chance.

###  III-B Hyperparameters

The proposed method has four hyperparameters that require tuning. In this
subsection we will discuss the affects of these hyperparameters. In the first
step we perform ODL, any pursuit algorithm can be used; we elected to use OMP
so we could select exactly how many non-zero elements to extract in the sparse
coding step. When performing dictionary learning you must specify the number
of atoms in the dictionary, the regularisation parameter, $\lambda$, and the
number of iterations, $T$.

The number of atoms in the dictionary must be greater than the dimensionality
of the input signals, we found a reasonable heuristic for the number of atoms
is double the dimensionality of the signals. $\lambda$ specifies how heavily
to prioritise the sparsity of the sparse coefficients, higher $\lambda$ values
correspond to greater sparsity. In general, we require the sparse coefficient
vectors to be as sparse as possible while still providing for accurate data
reconstruction. Finally, $T$ specifies the number of iterations over which to
learn the dictionary; $T$ is particularly data specific but ODL guarantees
convergence on a specific reconstruction error past-which further training is
unnecessary.

##  IV Experiments and Results

###  IV-A Domains

In this section we present our experimental results. We ran experiments on
three widely used HSIs, namely: Salinas, Indian Pines and Pavia
University[10], the ground truth labels of each of these images can be seen in
Fig. 2. HSIs are images with a large number of spectral channels, typically
over one-hundred. An HSI’s information can be seen as operating in a
3-dimensional cube, where along the two spatial dimensions it is simply a
series of 2-dimensional images, and along the third (or spectral) dimension,
it consists of images taken at various wavelengths. In this structure, each
pixel in the hyperspectral image can be represented as a continuous signal;
this makes HSIs ideal candidates for signal processing techniques. In this
work, we do not take advantage of the spatial information within the HSIs, so
the method can be extended to other domains within signal processing.

Figure 2: Left to right: ground-truth of Pavia University, Indian Pines and
Salinas

Although these three data sets are all HSIs, they each present their own
challenges. Salinas, captured by the AVIRIS spectrometer, is $83\times 86$ and
contains $204$ spectral bands. It contains six classes and one miscellaneous
class, given its simplicity, high-accuracy predictions are expected. Indian
Pines, also captured by AVIRIS, contains sixteen different classes and its
dimensions are $145\times 145\times 200$, it is considered far more
challenging than Salinas[10]. Since both Salinas and Indian Pines were
captured by AVIRIS, they have similar spectral resolutions. Pavia University,
captured by ROSIS, has only $103$ spectral bands, however, at $610\times 340$,
it is much larger than the other two images. Pavia University has nine classes
and so is between Salinas and Indian Pines in terms of label-complexity. These
three data sets are used in the field for a wide-range of experiments, as can
be seen in [6, 7] and [8].

###  IV-B Methodology

Our Python implementation of the method described in subsection III-A will be
made available as of $18$ October
2019111https://github.com/JoshuaDBruton/SparseCoefficientClustering. The data
sets used in our experiments are available as .mat files on the Grupo de
Inteligencia Computacional website222http://tiny.cc/f6gnez. As is common when
working with these HSIs, we ignored the miscellaneous classes when clustering.
All of the experiments described in this paper were validated across all three
images.

We began our experiments by training dictionaries, these dictionaries were
saved as .npy files and used for clustering at a later stage. This means that
the same dictionaries could be used when comparing different clustering
techniques. When clustering, the number of clusters was set to the number of
classes present in the ground-truth. For the experiments shown here, we used
OMP as the pursuit algorithm. OMP’s stopping criteria was always the number of
non-zero elements in the coefficients, this allows for fairer comparison to
other feature extraction techniques, like PCA, in which the number of
components is specified.

When describing our experiments we require some terms that are not available
succinctly in the field, and so we will elaborate on their meaning here. By
“discriminative dictionary” we mean a dictionary whose sparse coefficients
provide features that result in adequate clustering accuracy. By “size” of the
dictionary we mean the number of atoms specified for the dictionary. We term
the number of non-zero elements targeted in the sparse coding step of ODL the
“sparsity of the dictionary”.

###  IV-C Sparsity and Number of Components

We first undertook the task of establishing the correct size and sparsity to
target when creating a discriminative dictionary. We provide a baseline
clustering for Salinas in Fig. 3, this baseline was created by applying KMeans
to the orignal pixels of Salinas. Since the intention of this experiment was
to compare the discriminative power of the dictionaries, our only concerns
were the relative performance of the final sparse coefficient predictions and
how they compared to this baseline. This also justifies our use of KMeans
instead of a more complex clustering technique.

Figure 3: Ground-truth left, prediction of KMeans on original pixels right,
colour of classes is arbitrary, note significant error in bottom-right class

We proceeded with a grid-search over some possible parametrisations of ODL. We
trained 16 different dictionaries for each HSI, each with various sizes and
sparsities. These dictionaries and some visualisations of the training process
are available on Comet.ml333https://www.comet.ml/joshuabruton/honours-
project/view/. We then performed sparse coding using these dictionaries, and
used the computed sparse coefficients as inputs to KMeans. We displayed the
resultant predictions in a grid to see if any obvious patterns emerged.

A $3\times 3$ sample of the grid created for Salinas can be seen in Fig. 4.
Nine dictionaries are compared in the grid, all were trained for $5000$
iterations. The number of atoms in each dictionary varies over $(220,306,408)$
and the sparsity varies over $(10,5,2)$. Every combination of those values is
shown in the grid, with the smallest and sparsest dictionary in the bottom-
right and the largest and densest dictionary in the top-left.

Figure 4: Grid-search over ODL hyperparameters, increasing sparsity on x-axis,
increasing size on y-axis, the best result is highlighted in the bottom-right

As can be seen in the grid, the sparsest and smallest dictionary, targeting
$2$ non-zero elements in sparse coefficients with size $220$, provided the
best prediction. This prediction is similar to the baseline prediction
provided in Fig. 3 and many of the errors are common to both predictions. This
result was not anticipated, the denser and larger dictionaries have
substantially lower reconstruction errors. We initially expected lower
reconstruction error to correlate with predictions more similar to the
baseline. In hindsight, this result confirms the idea that the sparse
coefficients intend to represent the variance in the signals while the
dictionary represents their generic structure. It also confirms that sparse
coding prioritises the extraction of more discriminative features. Making the
sparse coefficients too large, or too dense, effectively allows OMP to
represent any noise within the signals.

###  IV-D Joint Sparse Representation

Before proceeding with experiments relating to SC, we tested the performance
of the most trivial implementation of JSR, as outlined in II-C, with KMeans.
We trained five dictionaries using JSR, all of them had size $330$ and
targeted $2$ non-zero elements in the sparse coding step. The tile-sizes of
the dictionaries varied over $(3,5,7,9,11)$. When using JSR, after the sparse
coding step, each pixel is represented by a tile of sparse coefficients. We
thought of two approaches for reducing this tile to one datum so it could be
clustered. Either taking the mean of the entire tile, or simply extracting the
centre sparse coefficient from the tile. We attempted both approaches, they
have similar results and a comparison can be seen in Fig. 5.

Figure 5: Mean of tile top, middle pixel of tile bottom

The predictions provided by JSR are not as accurate as we hypothesised. JSR’s
predictions seem to be less accurate than the KMeans baseline. There are more
complex approaches to JSR that may have produced better results, some of which
are described in [19]. However, most of these alternatives are supervised,
placing them out of the scope of this work.

###  IV-E Feature extraction

We now examine the benefit of using the most discriminative dictionaries we
found in experiment IV-C to produce sparse coefficients as input to SC. The
intention is that these sparse coefficients will allow for more effective
construction of the affinity matrices. We compared the results of four
approaches to spectral clustering; the first using the original signals, the
second using principal components from PCA, the third using non-negative
matrices computed from NMF, and the fourth using the sparse coefficients
computed from dictionaries. For this experiment, the three smallest and
sparsest dictionaries found in the previous experiment were used, one for each
HSI. All three targeted a sparsity of $5$ in the sparse coding step that
precedes training. The dictionaries trained on Salinas and Indian Pines both
had $220$ atoms and the dictionary trained on Pavia had $120$ atoms. The
simplicity of ODL’s hyperparameter selection is an important characteristic of
the method. We used PCA to extract $50$ principal components and clustered on
those components. We provide results for NMF with 8 components, a higher
number tended to provide grainy predictions. It is worth noting that both NMF
and PCA use larger representations than the one used for the sparse
coefficients, PCA significantly so.

The results of this comparison are shown in Table I. SC generally performs
poorly on the original pixels, this is because constructing the affinity
matrix from high-dimensional data is ineffective. PCA and sparse coding
effectively compensate for this by creating less complex but still meaningful
representations. The effect of sparse coding is less noticeable on Pavia
University, we conclude that this is because it has the fewest spectral
channels. Notably, spectral clustering on the sparse coefficients works most
effectively on the two images with the largest number of channels. Dictionary
learning and sparse coding are designed specifically for application to high-
dimensional signals, this explains the under-performance of PCA and NMF in
these cases. Relative to PCA and the sparse coefficients, NMF performs poorly
on the two images with the highest dimensional pixels but performs best on
Pavia University. The improvement our method provides is especially noticeable
on Salinas, with a $0.0399$ increase in the AMI. SC is well-suited to Salinas,
however, only the sparse coefficients provide features powerful enough to take
advantage of this compatibility. A visual comparison of the top three Salinas
predictions is shown in Fig. 6.

TABLE I: Comparison of different features in spectral clustering Adjusted
mutual information scorea
---
HSI | Original pixels | PCA | NMF | Sparse coefficients
Salinas | 0.8981 | 0.9088 | 0.8731 | 0.9487
Indian Pines | 0.6310 | 0.6489 | 0.626 | 0.6509
Pavia U | 0.8056 | 0.83 | 0.849 | 0.8065
aAn AMI of 1.0 signifies a perfect clustering, higher is better

Figure 6: Left to right: SC on original pixels, SC on 50 principal components,
SC with 5 non-zero elements

###  IV-F Discussion

Our experiments here show that it is possible to represent HSIs in the
significantly compressed form provided by ODL while still maintaining enough
information to cluster as or more effectively than in the original
hyperspectral pixel domain. Clustering the sparse coefficients is also more
computationally efficient as far fewer comparisons are needed. This
computational efficiency will be cancelled out if clustering is not performed
frequently, but could be beneficial when clustering similar data repetitively.
The benefit of using ODL in this scenario is that the same dictionary can be
continually updated as new data is added, without reconsidering what data was
previously seen. As has been shown by [2], dictionary learning does allow for
generalisability. This means it would also be possible to generalise a
dictionary trained on a similar data set to other datum, meaning that a new
dictionary does not have to be found for every new clustering problem. We also
intend on further investigating joint sparse representation. By taking
advantage of the available spatial information, joint sparse representation
should greatly improve spectral clustering accuracy in HSIs.

##  V Conclusion

In this paper, we proposed a method for clustering hyperspectral images via
dictionary learning and sparse coding. In particular, we investigated a number
of different clustering approaches for hyperspectral images. This includes
clustering in the original pixel domain, clustering in the sparse coefficient
domain, spectral clustering in the sparse coefficient domain, spectral
clustering using non-negative matrix factorisation, and clustering in the PCA
coefficient domain. Our experimental results show that using sparse
coefficients led to competitive or superior performance when compared to the
other methods considered in this paper.

## References

  * [1] J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online dictionary learning for sparse coding,” $26^{th}$ International Conference on Machine Learning, Montreal, pp. 689–696, 2009.
  * [2] M. Aharon, M. Elad, and A. Bruckstein, “K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation,” IEEE Transactions on Signal Processing, vol. 54, pp. 4311–4322, November 2006.
  * [3] K. Engan, S. O. Aase, and J. H. Husøy, “Multi-frame compression: theory and design,” EURASIP Signal Process, vol. 80, pp. 2121–2140, 2000.
  * [4] S. Mallat and Z. Zhang, “Matching pursuits with time-frequency dictionaries,” IEEE Trans. Signal Processing, vol. 41, pp. 3397–3415, 1993.
  * [5] S. Romano, N. Vinh, J. Bailey, and K. Verspoor, “Adjusting for chance clustering comparison measures,” Journal of Machine Learning Research, vol. 17, pp. 1–32, 2016.
  * [6] H. Wang and T. Celik, “Sparse representation-based hyperspectral image classification,” Signal, Image and Video Processing, vol. 12, pp. 1009–1017, July 2018.
  * [7] C. Cariou and K. Chehdi, “Unsupervised nearest neighbors clustering with application to hyperspectral images,” IEEE Journal of Selected Topics in Signal Processing, vol. 9, no. 6, pp. 1105–1116, September 2015.
  * [8] J. Murphy and M. Maggioni, “Unsupervised clustering with active learning of hyperspectral images with nonlinear diffusion,” IEEE Transactions on Geoscience and Remote Sensing, vol. 57, no. 3, pp. 1829–1845, 2018.
  * [9] P. Sprechmann and G. Sapiro, “Dictionary learning and sparse coding for unsupervised clustering,” IEEE international conference on acoustics, speech and signal processing, pp. 2042–2045, 2010.
  * [10] R. Green, M. Eastwood, C. Sarture, T. Chrien, M. Aronsson, B. Chippendale, J. Faust, B. Pavri, C. Chovit, and M. Solis, et al., “Imaging spectroscopy and the airborne visible/infrared imaging spectrometer (AVIRIS),” Remote sensing of environment, vol. 65, no. 3, pp. 227–248, 1998.
  * [11] X. Feng, “Robust spectral clustering via sparse representation,” Recent Applications in Data Clustering, p. 155, August 2018.
  * [12] J. Wright, Y. Ma, J. Mairal, G. Sapiro, T. Huang, and S. Yan, “Sparse representation for computer vision and pattern recognition,” Proceedings of the IEEE, vol. 98, no. 6, pp. 1031–1044, 2010.
  * [13] Y. Weiss, “Segmentation using eigenvectors: a unifying view,” Proceedings of the seventh IEEE international conference on computer vision, vol. 2, pp. 975–982, 1999.
  * [14] A. Ng, M. Jordan, and Y. Weiss, “On spectral clustering: analysis and an algorithm,” Advances in neural information processing systems, pp. 849–856, 2002.
  * [15] C. Bishop, “Pattern recognition and machine learning,” Springer, 2006.
  * [16] Z. Tang, S. Ding, Z. Li, and L. Jiang, “Dictionary learning based on nonnegative matrix factorization using parallel coordinate descent,” Abstract and Applied Analysis, vol. 2013, Article ID 259863, 2013.
  * [17] J. Tropp, A. Gilbert, and M. Strauss, “Algorithms for simultaneous sparse approximation. Part I: greedy pursuit,” Signal Processing,vol. 86, no. 3, pp. 572–588, Elsevier, 2006.
  * [18] J. Mairal, F. Bach, J. Ponce, and G. Sapiro, “Online learning for matrix factorization and sparse coding,” inria-00408716v1, 2009.
  * [19] Z. Wang, J. Liu, and J. Hue, “Joint sparse model-based discriminative K-SVD for hyperspectral image classification,” Signal Processing, vol. 133, pp. 144–155, Elsevier, 2017.


* Olson (1990) T. S. Olson, Annals of Physics 199, 18 (1990).
  * Salazar and Zannias (2020) J. F. Salazar and T. Zannias, International Journal of Modern Physics D 29, 2030010 (2020), arXiv:1904.04368 [gr-qc] .
  * Noronha  _et al._ (2021) J. Noronha, M. Spaliński, and E. Speranza, arXiv e-prints , arXiv:2105.01034 (2021), arXiv:2105.01034 [nucl-th] .
  * Landau and Lifshitz (2013b) L. Landau and E. Lifshitz,  _Fluid Mechanics_ , v. 6 (Elsevier Science, 2013).
  * Gavassino  _et al._ (2022) L. Gavassino, M. Antonelli, and B. Haskell, Phys. Rev. D 105, 045011 (2022).
  * Camelio  _et al._ (2022b) G. Camelio, L. Gavassino, M. Antonelli, S. Bernuzzi, and B. Haskell, arXiv e-prints , arXiv:2204.11810 (2022b), arXiv:2204.11810 [gr-qc] .
  * Krotscheck and Kundt (1978) E. Krotscheck and W. Kundt, Communications in Mathematical Physics 60, 171 (1978).
  * Heimburg (2017) T. Heimburg, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions) 19, 17331 (2017), arXiv:1608.06093 [physics.chem-ph] .
  * Hiscock and Lindblom (1988) W. A. Hiscock and L. Lindblom, Physics Letters A 131, 509 (1988).


# Demonstration of electron-nuclear decoupling at a spin clock transition

Krishnendu Kundu  1 Jia Chen  2,3,4 Silas Hoffman  2,3,4 Jonathan Marbey
1,4,5 Dorsa Komijani  1,5 Yan Duan  6 Alejandro Gaita-Ariño  6 John Stanton
3,4,7 Xiao-Guang Zhang  2,3,4 Hai-Ping Cheng  2,3,4 Stephen Hill1,4,5

###### Abstract

The ability to design quantum systems that decouple from environmental noise
sources is highly desirable for development of quantum technologies with
optimal coherence. The chemical tunability of electronic states in magnetic
molecules combined with advanced electron spin resonance techniques provides
excellent opportunities to address this problem. Indeed, so-called clock
transitions (CTs) have been shown to protect molecular spin qubits from
magnetic noise, giving rise to significantly enhanced coherence. Here we
conduct a spectroscopic and computational investigation of this physics,
focusing on the role of the nuclear bath. Away from the CT, linear coupling to
the nuclear degrees of freedom causes a modulation and decay of electronic
coherence, as quantified via electron spin echo signals generated
experimentally and _in silico_. Meanwhile, the effective hyperfine interaction
vanishes at the CT, resulting in electron-nuclear decoupling and an absence of
quantum information leakage to the nuclear bath, providing opportunities to
characterize other decoherence sources.

National High Magnetic Field Laboratory, Florida State University,
Tallahassee, FL 32310, USA

Department of Physics, University of Florida, Gainesville, FL 32611, USA

Quantum Theory Project, University of Florida, Gainesville, FL 32611, USA

Center for Molecular Magnetic Quantum Materials

Department of Physics, Florida State University, Tallahassee, FL 32306, USA

Instituto de Ciencia Molecular (ICMol), Universidad de Valencia, Paterna,
Spain

Department of Chemistry, University of Florida, Gainesville, FL 32611, USA

## Main

The synthetic tunability of molecular nanomagnets provides a versatile
platform for exploring and potentially harnessing unique physical attributes
that are of utility for the development of next-generation quantum information
(QI) technologies.1; 2; 3; 4 In particular, the electronic spin associated
with a magnetic molecule may serve as the computational basis for a quantum
bit, or qubit. However, as with any such system, protection from environmental
noise that causes decoherence is of critical importance, representing one of
the main hurdles on the path towards practical applications. In an attempt to
suppress one of the more stubborn sources of decoherence arising from
electron-nuclear interactions, various synthetic strategies have been employed
such as nuclear spin patterning5; 6 and the use of nuclear spin free
ligands.7; 8 However, demonstration of long phase memory (coherence) times
typically still requires extreme dilution in order to minimize electron spin-
spin dephasing.

Rather than modifying the spin bath, an alternative approach involves so-
called clock transitions (CTs)9 at which the electron spin resonance (ESR)
frequency is insensitive to the local magnetic induction and, therefore, does
not couple to the fluctuating magnetic environment. Such CTs occur at avoided
level crossings associated with the Zeeman splitting of qubit basis states.
This approach is well established in solid-state materials such as donor atoms
in silicon10; 11 or defect states in various other host crystals. 12; 13; 14;
15; 16 Our interest is in molecular systems, for which enhanced coherence was
demonstrated at a CT for a [Ho(W5O18)2]9- molecule by Shiddiq et al.17
Subsequently, CTs have been studied in other molecular systems19; 20; 21; 22;
23 and the effects of structural distortions have been analyzed theoretically
for several HoIII and VIV complexes.24

Here we directly investigate electron-nuclear coupling in the vicinity of a CT
by means of pulsed electron-spin-echo (ESE) measurements and numerical
modelling. Away from the CT, dipolar hyperfine coupling to the nuclear bath
results in periodic modulations of the electronic coherence—the so-called ESE
envelope modulation (ESEEM) effect.25 This modulation vanishes at the CT.
Theoretically, we consider a minimal model that can host a CT: an $S=1$ spin
subject to a relatively strong axial magnetic anisotropy, with an avoided
Zeeman level crossing generated by a weaker transverse interaction [Fig.
1(a)]. We treat coupling to the nuclear bath explicitly to reproduce the ESEEM
effect via quantum dynamics simulations. The parameters in our simplified
$S=1$ model are chosen to mimic the low energy physics of the [Ho(W5O18)2]9-
molecule, the only system for which ESEEM has been characterized as a function
of applied magnetic field, ${B}$0, in the vicinity of a CT. The simulations
compare favorably with experiment. Crucially, we demonstrate electron-nuclear
decoupling at the CT. Although the experiments focus on [Ho(W5O18)2]9-, our
model applies quite generally for the coupling of an electronic spin to a
finite nuclear bath. The combined study provides a microscopic view of the
mechanism via which an electron spin qubit couples to nearby nuclei, in
essence mediating leakage of quantum information to the nuclear bath.

Pulsed ESR, which is central to most spin-based QI device implementations,26
is an extremely powerful technique enabling both sample characterization and
quantum control. The simplest illustration involves the two-pulse Hahn echo
sequence,25; 27 where a coherent superposition of spin “up” and “down” states
is first generated via a $\pi/2$ rotation on the Bloch sphere, and then the
magnetization is allowed to evolve freely in the $xy$-plane; this evolution is
later inverted via application of a $\pi$-pulse, ideally refocusing any
dephasing that occurs due to static disorder, resulting in emission of an ESE
at time $2\tau$ after the initial $\pi/2$ pulse ($\tau$ is the delay between
pulses). A dynamic environment causes decoherence,28 which manifests as a
decay of the ESE intensity upon increasing $\tau$. Meanwhile, coherent
interactions with nearby quantum systems, e.g., other electrons or atomic
nuclei, can give rise to a modulation of the ESE intensity.25 In particular,
ESEEM arises due to excitation of formally forbidden nuclear transitions
during the pulsed ESE sequence, through hyperfine coupling to the central
electron spin. Here, “central” refers to spins that have been prepared in a
prescribed coherent quantum state, e.g., via application of a $\pi/2$ pulse.
ESEEM may therefore be used to characterize this aspect of the environment,
providing uniquely sensitive fingerprints of electron-nuclear decoherence
mechanisms.

In order to gain microscopic insights into electron-nuclear coupling in the
vicinity of a CT, ESEEM measurements were performed on a crystal of
Na9[Ho0.001Y0.999(W5O18)2]$\cdot n$H2O (hereon abbreviated HoW10), i.e.,
$0.1\%$ HoW10 doped into an isostructural non-magnetic YW10 host. HoIII
possesses a ground state spin-orbit coupled angular momentum, $J=L+S=8$. The
pseudo-axial coordination geometry imposed on the HoIII ion results in a
crystal field (CF) interaction that lifts the degeneracy of the $2J+1$
projection ($m_{J}$) states, giving rise to a singlet and a series of
$m_{J}\approx\pm i$ ($i=1$ to 8) quasi-doublets, with the $m_{J}=\pm 4$ ground
doublet lying $\approx$ 40 cm-1 below the first excited CF states.29; 30 A
weak tetragonal CF interaction is effective in generating an avoided Zeeman
level crossing between the $m_{J}=\pm 4$ basis states,17; 18 thus giving rise
to a 9.18 GHz CT. The hyperfine interaction involving the $I=\frac{7}{2}$
165Ho nuclear spin further splits the $m_{J}=\pm 4$ states into $(2I+1)=8$
pairs of $m_{I}$ sub-levels, resulting in eight avoided-crossings, i.e., eight
CTs, four either side of zero applied field (see Fig. 1 in Ref.
[Shiddiq2016]). For reasons explained in Ref. [Shiddiq2016], we focus here on
the lowest field CT (at $B_{0z}=23.6$ mT), which also gives the strongest
ESEEM; note that, due to a small sample misalignment, this occurs at at
$B_{0}=B_{\mathrm{min}}=25.5$ mT in the present investigation (see Methods).

ESE time traces recorded at a frequency of 9.18 GHz are shown in Fig. 2(a) for
different detuning fields ($\Delta B=B_{0}-B_{\mathrm{min}}$) from the CT,
revealing strong temporal modulations (ESEEM) at most detunings. The first
thing to note is the variation in decay time ($\equiv$ phase memory time,
$T_{\mathrm{m}}$) and modulation depth as a function of the detuning. Most
notably, there is a complete absence of ESEEM at zero detuning, i.e., at the
CT. Fast Fourier transforms (FFTs) of the time traces reveal three prominent
peaks, highlighted by the red, green and blue circles in Fig. 2(b). The
associated ESEEM frequencies are plotted as a function $B_{0}$ in Fig. 2(c);
superimposed on the data are the 1st and 2nd harmonics of the bare proton
Larmor frequency, $\nu_{\mathrm{H}}=\gamma_{\mathrm{H}}B_{0}$, where
$\gamma_{\mathrm{H}}=42.577$ MHz/T is the proton gyromagnetic ratio. The fact
that the average of the red/green data points coincides with
$\nu_{\mathrm{H}}$ and the blue data points with $2\nu_{\mathrm{H}}$ is a
strong indication that the ESEEM is caused by dipolar coupling to protons.
This is not surprising given the significant amount of water in the lattice of
[HoW10]$\cdot n$H2O ($n\approx 35$ in fully solvated crystals). Indeed, a
strong proton ESEEM effect is expected in this field range where the Ho-H
dipolar coupling strength is comparable to the proton Larmor frequency (see
below). By contrast, all other nuclei are predominantly non-magnetic, either
due to low $\gamma$-values or low abundance of magnetic isotopes.

A qualitative understanding of the ESEEM spectrum is obtained by first
considering the simplest possible case of coupled $S=\frac{1}{2}$ and
$I=\frac{1}{2}$ spins in the high-field limit in which $\nu_{\mathrm{H}}\gg
A$, where $A$ (= $A_{zz}/h$, $A_{zz}$ is the $z$-component of the hyperfine
tensor) quantifies the bare dipolar coupling strength in frequency units.
ESEEM arises due to excitation of formally forbidden zero- and double-quantum
transitions that rotate coupled electron and nuclear spins.31 The modulation
results from combinations of the allowed
($\nu_{a}=\gamma_{e}B_{0}\pm\frac{1}{2}A$, $\gamma_{e}$ is the electron
gyromagnetic ratio) and formally forbidden
($\nu_{f}=\gamma_{e}B_{0}\pm\nu_{\mathrm{H}}$) transition frequencies at:
$|\nu^{\pm}_{a}-\nu^{\mp}_{a}|=A$,
$|\nu^{+}_{f}-\nu^{\pm}_{a}|=|\nu^{-}_{f}-\nu^{\pm}_{a}|=\nu_{\mathrm{H}}\pm\frac{1}{2}A$
and $|\nu^{\pm}_{f}-\nu^{\mp}_{f}|=2\nu_{\mathrm{H}}$.31 One may then
understand the red/green data points in Fig. 2(c) as being due to the
hyperfine coupled proton frequencies,
$\nu_{\mathrm{H}}\pm\frac{1}{2}A^{\mathrm{eff}}$, where $A^{\textrm{eff}}$ is
an effective coupling strength on account of the new physics that emerges at
the CT ($A^{\textrm{eff}}$ is further renormalized for HoW10 due to the fact
that $S\neq\frac{1}{2}$). Crucially, $A^{\mathrm{eff}}\rightarrow 0$ at the
CT, which may be understood as being a consequence of the effective electron
gyromagnetic ratio, $\gamma^{\mathrm{eff}}_{e}$, crossing through zero at
$B_{0}=B_{\textrm{min}}$ ($\gamma^{\mathrm{eff}}_{e}\propto df/dB_{0}$ or
$\langle\hat{S}_{z}\rangle$, the $z$-component spin expectation value), as
illustrated in Fig. 1(b), where the ESR (clock) frequency couples
quadratically to $B_{0}$ at the avoided crossing (CT), in contrast to the
usual linear coupling far from the CT. This is why the ordering of red and
green circles switches at the CT, i.e., there is a smooth evolution of
$A^{\mathrm{eff}}$ ($\propto\gamma^{\mathrm{eff}}_{e}$) such that it switches
sign at the CT. Remarkably, to first order, this implies that the effective
dipolar coupling to protons vanishes right at the CT; hence the ESEEM effect
also vanishes at the CT, as does the electron-nuclear decoherence, leading to
the steep rise in $T_{\mathrm{m}}$ as one approaches the CT
[$=8.43(6)~{}\mathrm{\mu}$s @ the CT].17 Meanwhile, the ESEEM modulation depth
grows with the detuning, $\Delta B$ (i.e., with $\gamma^{\mathrm{eff}}_{e}$),
away from the CT, as does the electron-nuclear contribution to the central
spin decoherence, i.e., $T_{\mathrm{m}}$ decreases to
$\sim$$1~{}\mathrm{\mu}$s far from the CT.17

The ESEEM effect is ultimately governed by collective coupling of the HoIII
ion to the entire nuclear bath. However, the $1/r^{3}$ dependence of the
dipolar interaction and large value of $\gamma_{\mathrm{H}}$ in comparison to
other nuclei results in a spectrum that is dominated by nearby protons,32; 33
the closest of which are $\sim$4 ${\mathrm{\AA}}$ from the central HoIII
ion.18 At this separation and in the linear Zeeman regime [$\Delta B_{0}>200$
mT from the CT in Fig. 1(b)], the maximum Ho-H dipolar coupling strength,
$A^{\mathrm{max}}\approx 3$ MHz (=
$2\mu_{\mathrm{o}}\mu_{\mathrm{Ho}}\mu_{\mathrm{H}}/4\pi hr^{3}$); this
assumes $m_{J}=\pm 4$ for the ground state of HoIII. The experimental results
displayed in Fig. 2 remain very far from this linear regime, which is why the
separation of the red and green data points ($A^{\mathrm{eff}}\approx 0.4$
MHz) is well below the maximum. Meanwhile, ESEEM measurements far from the CT
are hampered by the short phase memory time. Nevertheless, one would expect to
observe an ESEEM effect in this field range because $A^{\mathrm{eff}}$ is of
the same order as the proton Larmor frequency, $\nu_{H}=1.1$ MHz at 25.5 mT.
Indeed, ESEEM is also observed at the 2nd ($\nu_{H}=3.3$ MHz) and 3rd
($\nu_{H}=5.4$ MHz) CTs. Although the effect is less pronounced, the same
qualitative behavior is found, i.e., a vanishing of ESEEM at each CT and
harmonic content centered at $\nu_{\mathrm{H}}$ and $2\nu_{\mathrm{H}}$.
Therefore, the enhanced coherence in the vicinity of the CTs provides a window
through which to observe the ESEEM, which ultimately vanishes right at the CT
because $\gamma^{\mathrm{eff}}_{e}\rightarrow 0$. We note that no modulation
is discernible at the 4th CT ($\nu_{H}=7.6$ MHz), presumably because the
effective dipolar coupling is just too weak in comparison to
$\nu_{\mathrm{H}}$.

In order to gain microscopic understanding, we developed a simplified
Hamiltonian for a central electron spin coupled to a finite proton spin bath.
In order to preserve computational resources for the bath, we model the
electronic system as an $S=1$ spin with longitudinal and transverse anisotropy
[Fig. 1(a)]:

$\displaystyle\hat{H}_{S}=D[\hat{S}_{z}^{2}-\tfrac{1}{3}S(S+1)]+E(\hat{S}_{x}^{2}-\hat{S}_{y}^{2})+\gamma_{e}(B_{0}-B_{\textrm{min}})\hat{S}_{z}\,,$
(1)

where $\hat{S}_{j}$ are spin-1 generators of rotation about axis $j$, while
$D$ and $E$ are the 2nd order axial and rhombic zero-field splitting
(anisotropy) parameters, respectively. $B_{\textrm{min}}$ is introduced to
shift the CT away from $B_{0}=0$, mimicking the effect of the on-site
hyperfine interaction with the 165Ho nuclear spin; note that this field does
not act on the proton bath. The eigenvectors of Eq. (1) at the CT (i.e., when
$\Delta B=0$) are
$|\pm\rangle=\tfrac{1}{\sqrt{2}}(|\uparrow\rangle\pm|\downarrow\rangle)$ and
$|0\rangle$, with energies $-\tfrac{1}{3}|D|\pm E$ and $+\tfrac{2}{3}|D|$,
respectively. Here, $|\uparrow\rangle$, $|\downarrow\rangle$, and $|0\rangle$
are the states with $\langle\hat{S}_{z}\rangle=\pm 1$ and
$\langle\hat{S}_{z}\rangle=0$, respectively.

We set $D=-45$ GHz, $|E|=4.5$ GHz and $B_{\textrm{min}}=23.5$ mT in order to
mimic the actual low-energy electronic structure of HoW10. These parameters
ensure the same CT frequency, $\Delta=2E=9$ GHz, the same curvature of the two
lowest lying levels, and a sizeable separation to the $|0\rangle$ state (Fig.
1). As an aside, because $|\pm\rangle$ are energetically well-separated from
$|0\rangle$ in the vicinity of the CT, we can project onto the two-dimensional
subspace defined by the former, wherein,

$\displaystyle\hat{S}_{z}^{2}\rightarrow\mathbb{1}\,,\,\,\hat{S}_{z}\rightarrow\sigma_{x}\,,\,\,\hat{S}_{x}^{2}-\hat{S}_{y}^{2}\rightarrow\sigma_{z}\,,\,\,\\{\hat{S}_{x},\hat{S}_{y}\\}\rightarrow
2\sigma_{y}\,,\,\,$ $\displaystyle\hat{S}_{x}\rightarrow
0\,,\,\,\hat{S}_{y}\rightarrow
0\,,\,\,\\{\hat{S}_{y},\hat{S}_{z}\\}\rightarrow
0\,,\,\,\\{\hat{S}_{z},\hat{S}_{x}\\}\rightarrow 0\,.$ (2)

Using this notation, the Hamiltonian reduces to $\hat{H}_{S}\rightarrow
E\sigma_{z}+\gamma\Delta B\sigma_{x}$, which precisely maps onto a
‘fictitious’ spin-$\frac{1}{2}$ model subjected to an effective magnetic field
in the $xz$-plane.34 The eigenvectors, which are quantized along the effective
field direction, are still denoted $|\pm\rangle$, although these are no longer
equally weighted mixtures of $|\uparrow\rangle$ and $|\downarrow\rangle$ upon
detuning from the CT. Nevertheless, at the CT ($\Delta B=0$), one may
visualize qubit operations within this subspace in terms of pure rotations
around the $j$th axis of the Bloch sphere defined by $|\pm\rangle$, according
to the Pauli matrices, $\sigma_{j}$; the corresponding spin-1 operators are
then easily found from Eq. (2). This mapping is helpful in understanding the
simulated Hahn-echo sequence (see Methods), as there is no simple analogy to
the $S=\frac{1}{2}$ rotating frame for the actual $S=1$ spin dynamics.

The nuclear spin bath, which ultimately causes decoherence and the observed
ESEEM effect, is described by $N$ protons coupled via dipolar interactions to
the central $S=1$ state,

$\displaystyle\hat{H}_{SI}=\hat{S}_{z}\sum_{m=1}^{N}\left[A_{sc}^{m}\hat{I}_{z}^{m}+A_{psc}^{m}(\hat{I}^{m}_{x}+\hat{I}^{m}_{y})\right]\,.$
(3)

Here, we employ secular ($sc$) and pseudosecular ($psc$) approximations with
phenomenological couplings $A_{sc}^{m}$ and $A_{psc}^{m}$, respectively; the
$\hat{I}_{j}^{m}$ are generators that rotate the spin of the $m$th proton
around axis $j$. The pseudosecular interaction is often ignored due to
averaging brought about by the mismatch in the proton Larmor and hyperfine
frequencies. However, as previously discussed, this is not the case at the
first CT. Indeed, the pseudosecular interaction turns out to be essential to
the ESEEM effect because it is responsible for driving formally forbidden
nuclear transitions during the Hahn echo sequence.31 Meanwhile, the protons
also undergo their own dynamics, independent of the central spin, according to

$\displaystyle\hat{H}_{I}=-\sum_{m\neq
n}D_{mn}(3\cos^{2}\theta_{mn}-1)[2\hat{I}_{z}^{m}\hat{I}_{z}^{n}-\hat{I}_{x}^{m}\hat{I}_{x}^{n}-\hat{I}_{y}^{m}\hat{I}_{y}^{n}]-\gamma_{\mathrm{H}}B_{0}\sum_{m=1}^{N}\hat{I}_{z}^{m}\,.$
(4)

That is, each proton in the bath undergoes Larmor precession at a bare
frequency $\gamma_{\mathrm{H}}B_{0}$, and couples to other protons via a
dipolar interaction of strength $D_{mn}$ ($\sim$10 kHz); $\theta_{mn}$ is the
angle between ${B}$0 and the vector joining protons $m$ and $n$. Energy
conserving proton flip-flop processes, driven by the
$(\hat{I}_{x}^{m}\hat{I}_{x}^{n}+\hat{I}_{y}^{m}\hat{I}_{y}^{n})$ term, are
central to the electron spin decoherence process.10; 28; 32; 33 To simulate
the ESEEM, we numerically recreate the two-pulse Hahn echo sequence _in
silico_ by performing a time evolution according to the total Hamiltonian,
$\hat{H}_{tot}=\hat{H}_{S}+\hat{H}_{SI}+\hat{H}_{I}$ (see Methods).

As a warm up, we first consider the simple case of a single proton ($N=1$)
coupled to the central $S=1$ spin, with $A=A_{sc}=2A_{psc}=1$ MHz. Fig. 3
displays FFTs of the Hahn echo simulations for several detuning fields [inset
to (a) displays a representative time trace]. In analogy to the
$S=\frac{1}{2}$ case, we associate the lowest frequency FFT peak, and the
splitting of the peaks either side of $\nu_{\textrm{H}}$, with the effective
hyperfine interaction strength, $A^{\mathrm{eff}}$; the inset to Fig. 3(b)
plots this frequency as a function of $B_{0}-B_{\textrm{min}}$. As can clearly
be seen, and in analogy with the experiments, $A^{\mathrm{eff}}\rightarrow 0$
at the CT; indeed, the modulation (not shown) is also zero at the CT.
Moreover, far from the CT, such that
$\gamma_{e}|B_{0}-B_{\textrm{min}}|\gg|E|/h$, $A^{\mathrm{eff}}\rightarrow
2A$; the factor of two is due to renormalization because $S=1$ as opposed to
$\frac{1}{2}$. Thus, in the high-field limit, FFT peaks occur at $2A$,
$\nu_{\mathrm{H}}\pm A$ and $2\nu_{\mathrm{H}}$. Superimposed on the data in
the inset to Fig. 3(b) is a phenomenological fit that assumes
$A^{\mathrm{eff}}\propto\gamma_{e}^{\mathrm{eff}}$, deduced from $df/dB_{0}$
via Eq. (1). This confirms the idea that the variation in
$\gamma_{e}^{\mathrm{eff}}$ (or $\langle\hat{S}_{z}\rangle$) in the vicinity
the CT governs the dipolar coupling of the central spin to the nearby proton.
The final thing to note from the inset to Fig. 3(a) is the absence of
decoherence, i.e., the peak ESE intensity does not decay. This is because the
two-spin system executes perfectly coherent coupled dynamics, with no quantum
phase leakage, i.e., there is no bath associated with this model.

In order to better capture the physics associated with the spin bath, we
extend the model to $N=7$ nuclear spins with a distribution of dipolar
couplings to the central spin [(Fig. 4(a)], enabling simulations of the ESEEM
on reasonable timescales whilst also capturing the emergence of decoherence;
we set $\langle{A^{m}_{sc}}\rangle=2\langle{A^{m}_{psc}}\rangle=8$ MHz to best
reproduce the experimental results (see Methods for further details). Time
traces for several detunings either side of the CT are displayed in Fig. 4(b).
As can be seen, the simulations qualitatively reproduce the experimental
results in Fig. 2. A very clear ESEEM effect is observed that more-or-less
vanishes at the CT. Moreover, the modulation depth increases with the
detuning, $\Delta B$. The time traces also exhibit a very apparent decay in
the coherence of the central spin dynamics, with a phase memory time,
$T_{\mathrm{m}}$, that clearly diverges at the CT, i.e., the finite spin bath
model causes decoherence of the central spin. Remarkably for such a simplified
model, even the decoherence timescale is of the same order as the experiments.
The only exception is at zero detuning, where the numerical decay is
considerably flatter than the experiments. The residual decoherence observed
at the CT in experiments is attributed to spin-lattice relaxation,17 which is
not included in our model; we comment on this further below. Fourier tranforms
of the numerical time traces are displayed in Fig. 4(c). Again, agreement with
experiment is remarkably good. Indeed, a plot of the center frequencies of the
main FFT peaks as a function of detuning, $B_{0}-B_{\mathrm{min}}$ [Fig.
4(d)], reveals identical behavior to the experiments, i.e., a pair of peaks at
$\nu_{\mathrm{H}}\pm\frac{1}{2}A^{\mathrm{eff}}$ and a higher frequency peak
at $\sim 2\nu_{\mathrm{H}}$; the peaks have been color coded in the same way
as in Fig. 4(c). Once again, it can be seen that $A^{\mathrm{eff}}\rightarrow
0$ at $B_{0}=B_{\mathrm{min}}$, and increases with detuning from the CT.

The present experimental and theoretical investigation clearly demonstrates
effective decoupling of an electron spin qubit from the surrounding nuclear
bath at a CT, going beyond previous studies that simply show evidence for
enhanced coherence.17; 18 In fact, our simulations reveal a pronounced
enhancement in $T_{\mathrm{m}}$ at the CT, whereas the experiments on HoW10
indicate that coherence is limited there by other factors. The primary culprit
is spin-lattice ($T_{1}$) relaxation.17 In particular, molecular vibrations
that couple directly to the CF interactions(s) responsible for the CT (Fig. 1)
may be expected drive spin-lattice relaxation,18; 35; 36 an effect not
included in our model. However, weak decoherence is observed even at the CT in
our numerical simulations. We attribute this to $2^{\rm{nd}}$-order coupling,
$d^{2}f/dB_{0}^{2}=\gamma_{e}^{2}/\Delta$, i.e., $df/dB_{0}$ vanishes only
precisely at the CT, and the HoW10 qubit is therefore exposed to weak 1H
dipolar field fluctuations either side of $B_{\mathrm{min}}$. This suggests
that electron-nuclear decoupling should improve upon increasing the CT
frequency, since the $2^{\rm{nd}}$-order coupling scales inversely with
$\Delta$.

Electron spin-spin interactions have also been omitted from our model, since
we consider only one HoIII ion. One may expect the secular part of this
interaction (i.e., $\hat{S}^{m}_{z}\hat{S}^{n}_{z}$) to decouple at a CT in
exactly the same way that the proton bath decouples in this study, provided
that the interaction is not too strong. As noted above, perfect decoupling
occurs only to first-order ($df/dB_{0}\rightarrow 0$) at the CT. However,
$2^{\rm{nd}}$-order coupling should be weak if the spin-spin interaction
strength is substantially weaker than the CT frequency ($\Delta=2E$),10; 23 as
is the case for well-separated ($>$nm) qubits. Meanwhile, although one may
safely ignore angular momentum conserving electron-nuclear dipolar flip-flop
processes in the present work because of the vastly different CT ($\Delta$)
and proton Larmor ($\gamma_{\mathrm{H}}B_{0}$) frequencies, this is not the
case for electron spin-spin interactions. Dipolar coupling within arrays of
nominally identical qubits will cause decoherence due to flip-flop processes
between resonant electron spins ($\Delta_{1}=\Delta_{2}$) via the
$\hat{S}^{m}_{x}\hat{S}^{n}_{x}+\hat{S}^{m}_{y}\hat{S}^{n}_{y}$
interaction.10; 37 Correctly modeling this physics is more challenging,
requiring a much larger bath with resonant and non-resonant qubits, due both
to disorder (distributions in $\Delta$) and a dynamic distribution of dipolar
interactions within the ensemble. Such a model contains complex many-body
physics that lies outside of the realm of the present investigation.

One may anticipate that future QI devices based on molecular spins will
feature controllable entangling interactions between individual qubits.1
Crucially, this control would enable mitigation of resonant electron-electron
spin flip-flop processes. Likewise, quantum sensing applications involving
single qubits are immune to this mode of decoherence. However, it is virtually
impossible to remove all sources of magnetic noise (particularly due to the
nuclear bath), whilst maintaining the flexibility that molecular design
principles allow. The present investigation therefore highlights the
importance of CTs for suppressing electron-nuclear spin-spin decoherence.
Moreover, one may expect these principles to apply quite generally to any type
of CT. In this regard, hyperfine CTs show the most promise due to weaker
coupling to molecular vibrations.23

## Methods

###  0.1 Experimental details.

Since extensive discussions of sample preparation and handling, experimental
setup and conditions, as well as the electronic properties that give rise to
CTs in HoW10 have been presented previously,29; 17 only brief descriptions of
essential details are given here. Single crystals of
Na9[Ho0.001Y0.999(W5O18)2]$\cdot n$H2O were prepared according to the method
described in Ref. [AlDamen09]. ESEEM measurements were performed using a
commercial Bruker E680 X-band spectrometer equipped with a cylindrical TE011
dielectric resonator (model ER 4118 X-MD5, with an unloaded center frequency
of 9.75 GHz), which was overcoupled to increase bandwidth and, thus, allow
measurements at frequencies down to 9.1 GHz.17; 23 The sample temperature was
controlled using an Oxford Instruments CF935 helium flow cryostat and ITC503
temperature controller.

All of the data presented in this study (Fig. 2) were obtained for a single
crystal. However, similar ESEEM behavior has been observed in experiments
performed on many other crystals of varying HoIII concentration.17 Although
_in situ_ rotation of the crystal about a single-axis is possible, the low
symmetry HoW10 structure29 and the need for rapid sample loading to avoid
degradation due to solvent loss resulted in an $\sim$22.5∘ misalignment
between the applied field and the molecular Ising axis. This simply leads to a
re-scaling of the CT fields: in this study, the lowest field CT occurs at
$B_{0}=B_{\mathrm{min}}=25.5$ mT, which is equivalent to a longitudinal field,
$B_{0z}=23.6$ mT, where the $z$-direction defines the approximate HoW10 four-
fold symmetry axis. ESEEM results were derived from ESE decay curves generated
using a standard two-pulse Hahn-echo sequence ($\pi/2-\tau-\pi-\tau-$ echo) as
a function of detuning from the CT field, $B_{\mathrm{min}}$. The frequency
domain plots in Fig. 2(b) were obtained by performing FFTs of the time traces,
which were zero padded by twice the number of data points and further smoothed
using a 5 point average.

The spin Hamiltonian of the HoW10 molecule may be described in terms of set of
axial CF parameters, $B_{k}^{q}$ ($k=2,4,6$, representing the rank of the
associated CF operator, $\hat{O}_{k}^{q}$, and $q=0$ the rotational order).
Distortions away from the approximate $D_{4d}$ point symmetry of the HoW10
molecule engage the tetragonal CF interaction,
$B_{4}^{4}\hat{O}_{4}^{4}\propto(\hat{S}_{+}^{4}+\hat{S}_{-}^{4})$,29 which is
effective in generating avoided crossings between the eight hyperfine sub-
level pairs associated with the $m_{J}=\pm 4$ ground doublet, resulting in CTs
at magnetic fields, $B_{\mathrm{min}}=\pm 23.6$, $\pm 70.9$, $\pm 118.1$ and
$\pm 165.4$ mT (for an applied field, $B_{0}$, parallel to the molecular
$z$-axis).17 The W and O nuclei in the HoW10 molecular core are predominantly
non-magnetic, with the exception of 17O ($I=\frac{5}{2}$, $\gamma=5.77$ MHz/T)
and 183W ($I=\frac{1}{2}$, $\gamma=1.77$ MHz/T) with $0.04\%$ and $14.3\%$
natural abundance, respectively. Moreover, their associated $\gamma$-values,
along with those of the more distant 23Na and 89Y nuclei (both $100\%$
abundance) are considerably smaller than those of the proton. Consequently,
one would not expect to see strong ESEEM from coupling to these other nuclei,
i.e., assignment of the observed ESEEM to protons is unambiguous.

###  0.2 Theoretical details.

The two-pulse Hahn echo sequence was recreated _in silico_ by performing a
time evolution according to the total Hamiltonian,
$\hat{H}_{tot}=\hat{H}_{S}+\hat{H}_{SI}+\hat{H}_{I}$.33 The initial density
matrix at thermal equilibrium was defined in the lab frame as

$\displaystyle\rho_{eq}=\frac{\exp(-\beta\hat{H}_{tot})}{\mathrm{Tr}(\exp(-\beta\hat{H}_{tot}))},$
(5)

where $\beta=h/k_{B}T$ and $T=5$ K. Instantaneous $\pi/2$ and $\pi$ pulses
were performed according to the procedure described in the following
paragraph. The density matrix was then allowed to evolve according to
$\hat{H}_{tot}$ for a time interval $\tau$ after each pulse. Finally, the echo
intensity was evaluated by computing the expectation value of the
$z$-component of the HoIII magnetization in the lab frame,
$\textrm{Tr}_{I}(\rho\hat{S}_{z})$, with the trace taken over the nuclear
states. Exact matrix diagonalization demands considerable computational
resources. Therefore, in order to carry out these calculations on reasonable
time scales, a number of compromises were necessary. Foremost among these was
the limitation on the size of the spin bath to $N=7$ protons. Meanwhile, based
on _a priori_ knowledge of the spin dynamics, we could also optimize the time
step and duration of the simulations, i.e., the time step (100 ns) results in
a frequency cut-off, which we set to well above the $2\nu_{\mathrm{H}}$
frequency seen in the experimental spectra (Fig. 2), and the duration (100
$\mathrm{\mu s}$) was chosen to ensure a FFT resolution comparable to the
experiments.

As discussed in the main text, the low energy $|\pm\rangle$ eigenvectors at
the CT are not the usual $|\uparrow\rangle$ and $|\downarrow\rangle$ states
relevant to the $S=\frac{1}{2}$ case; indeed, there is no simple rotating
frame analogy that can easily be visualized in the case of the ‘real’ $S=1$
system. One must therefore take care applying appropriate $\pi/2$\- and
$\pi$-pulses in order to generate the echo. In fact, one may reduce the
problem to the simple Bloch sphere picture via projection onto the two-
dimensional $|\pm\rangle$ subspace according to Eq. (2), i.e., a ‘fictitious’
spin-$\frac{1}{2}$ subjected to an effective magnetic field in the $xz$-plane
($\hat{H}_{S}\rightarrow E\sigma_{z}+\gamma\Delta B\sigma_{x}$).34 The
appropriate pulses can then be implemented via rotations about any axis that
is perpendicular to the effective field, $\vec{B}^{eff}$
($=\sqrt{E^{2}+(\gamma\Delta B)^{2}}$). Exactly at the CT ($\Delta B=0$),
where $\vec{B}^{eff}\parallel z$, this is easily achieved using either the
pure $\sigma_{x}$ or $\sigma_{y}$ Pauli matrix, corresponding to the spin-1
operators $\hat{S}_{z}$ and $\\{\hat{S}_{x},\hat{S}_{y}\\}$, respectively.
Away from the CT, $\vec{B}^{eff}$ tilts towards $x$ within the $|\pm\rangle$
subspace. We therefore employ a pure $\sigma_{y}$ rotation, which does not
depend on the orientation of $\vec{B}^{eff}$, i.e., we implement pulses of the
form $\exp[i\phi\\{\hat{S}_{x},\hat{S}_{y}\\}/2]$, where $\phi$ denotes the
rotation angle in radians. Although the $\\{\hat{S}_{x},\hat{S}_{y}\\}$
operator has no direct correspondence with the microwave $B_{1}$ field
employed in the experiments, it conveniently achieves the desired result.
Moreover, it is formally equivalent to operating with $\hat{S}_{z}$ at the CT,
which does correspond directly to the experimental parallel mode $B_{1}$
field. However, upon moving away from the CT, the ideal magnetic
$\hat{S}_{z}$-pulse evolves with the applied field, $B_{0}$, as the
eigenvectors acquire unequal $|\uparrow\rangle$ and $|\downarrow\rangle$
weights. Indeed, the durations of the $\pi/2$ and $\pi$ pulses employed in the
real experiments had to be optimized for each field step, something that could
be avoided in the simulations by implementing pure $\sigma_{y}$ rotations.

Additional subtleties of the calculations concern the precise details of the
microscopic interactions. For example, in order to recreate a realistic proton
bath, a distribution of electron-nuclear hyperfine coupling strengths,
$A^{m}=A_{sc}^{m}=2A_{psc}^{m}$ ($m=1$ to $N$), was generated with random
values in the range from 7 to 9 MHz such that $\langle A\rangle=8$ MHz.
Likewise, a distribution of proton-proton dipolar interactions was implemented
by fixing the coupling strength in Eq. (4),
$D_{mn}=\mu_{\mathrm{o}}\mu_{\mathrm{H}}^{2}/8\pi hr^{3}\approx 10$ kHz
($\equiv 1.8$ $\mathrm{\AA}$ distance), and randomizing the angle
$\theta_{mn}$. To compensate for the small size of the nuclear bath, the
simulations were repeated ten times for different $A^{m}$ and $\theta_{mn}$
randomizations, then averaged; this approach is obviously vastly more
efficient computationally compared to increasing the bath size. Not only do
these measures better mimic the real [HoW10]$\cdot n$H2O system, they avoid
the highly unphysical situation in which the seven protons are
indistinguishable, with identical couplings to the central spin and to each
other. The 8 MHz value for $\langle A\rangle$ was chosen so as to reproduce
the effective hyperfine interaction seen in the experiments, i.e., splitting
of the red/blue circles in Figs. 2c and 4d. This corresponds to a Ho-1H
separation of $\sim 2.9$ $\mathrm{\AA}$, which is a little under the known
closest distance ($\sim 4$ $\mathrm{\AA}$), which we attribute to the fact
that the model under-counts the number of nearby protons by about an order of
magnitude, i.e., $n=35$ H2O molecules, or 70 protons per HoIII ion.
Consequently, the smaller distance employed in the simulations effectively
renormalizes the collective hyperfine coupling strength. Finally, the
frequency-domain plots in Fig. 4(c,d) were obtained by first subtracting a
stretched exponential background from the time traces [Fig. 4(b)], then
performing FFTs of the residual ESEEM modulations; Fourier-transform filtering
was used to smooth the resulting frequency-domain plots.

In spite of the aforementioned simplifying assumptions, the employed model
captures the essential physics associated with the ESEEM effect in the
vicinity of a realistic CT. Moreover, the simulations reproduce the
experimentally observed electron-nuclear decoherence. Approximate cluster
correlation expansion (CCE) methods are able to consider a much larger and
more realistic bath consisting of thousands of protons. Indeed, such studies
applied to simple spin-$\frac{1}{2}$ qubits (with no CT) obtain essentially
perfect quantitative agreement with experimental phase memory times.32
However, they also reveal that decoherence is dominated by stochastic flip-
flop processes associated with proton pairs that are relatively close (a few
$\mathrm{\AA}$) to the central electron spin. As such, the exact quantum
calculations considered here contain the same ingredients. It is therefore
unsurprising that the obtained phase memory times agree with experiment to
within approximately a factor of two. Indeed, the $N=7$ model enables
exploration of many other microscopic aspects of the bath that influence
decoherence. We hope to explore this further in the future. We wish to
emphasize, however, that it was not our original intent to quantitatively
reproduce the decoherence, but rather to qualitatively reproduce the ESEEM
effect in the vicinity of a CT, something that this investigation has most
definitely achieved.

## Code availability

The computer code to run our simulation is available from K. K. upon request.

## References

  * (1) Gaita-Ariño, A., Luis, F., Hill, S. & Coronado, E.  Molecular spins for quantum computation.  _Nature Chemistry_ 11, 301–309 (2019).  URL https://doi.org/10.1038/s41557-019-0232-y.
  * (2) Atzori, M. & Sessoli, R.  The second quantum revolution: Role and challenges of molecular chemistry.  _Journal of the American Chemical Society_ 141, 11339–11352 (2019).  URL https://doi.org/10.1021/jacs.9b00984.
  * (3) Godfrin, C. _et al._ Operating quantum states in single magnetic molecules: Implementation of Grover’s quantum algorithm.  _Phys. Rev. Lett._ 119, 187702 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevLett.119.187702.
  * (4) Bayliss, S. L. _et al._ Optically addressable molecular spins for quantum information processing.  _Science_ 370, 1309–1312 (2020).  URL https://science.sciencemag.org/content/370/6522/1309.
  * (5) Graham, M. J., Yu, C.-J., Krzyaniak, M. D., Wasielewski, M. R. & Freedman, D. E.  Synthetic approach to determine the effect of nuclear spin distance on electronic spin decoherence.  _Journal of the American Chemical Society_ 139, 3196–3201 (2017).  URL https://doi.org/10.1021/jacs.6b13030.
  * (6) Jackson, C. E., Lin, C.-Y., Johnson, S. H., van Tol, J. & Zadrozny, J. M.  Nuclear-spin-pattern control of electron-spin dynamics in a series of V(IV) complexes.  _Chem. Sci._ 10, 8447–8454 (2019).  URL http://dx.doi.org/10.1039/C9SC02899D.
  * (7) Zadrozny, J. M., Niklas, J., Poluektov, O. G. & Freedman, D. E.  Millisecond coherence time in a tunable molecular electronic spin qubit.  _ACS Central Science_ 1, 488–492 (2015).  URL https://doi.org/10.1021/acscentsci.5b00338.
  * (8) Yu, C.-J. _et al._ Long coherence times in nuclear spin-free vanadyl qubits.  _Journal of the American Chemical Society_ 138, 14678–14685 (2016).  URL https://doi.org/10.1021/jacs.6b08467.
  * (9) Bollinger, J. J., Prestage, J. D., Itano, W. M. & Wineland, D. J.  Laser-cooled-atomic frequency standard.  _Phys. Rev. Lett._ 54, 1000–1003 (1985).  URL https://link.aps.org/doi/10.1103/PhysRevLett.54.1000.
  * (10) Wolfowicz, G. _et al._ Atomic clock transitions in silicon-based spin qubits.  _Nature Nanotechnology_ 8, 561–564 (2013).  URL https://doi.org/10.1038/nnano.2013.117.
  * (11) Morello, A., Pla, J. J., Bertet, P. & Jamieson, D. N.  Donor spins in silicon for quantum technologies.  _Advanced Quantum Technologies_ 3, 2000005.  URL https://onlinelibrary.wiley.com/doi/abs/10.1002/qute.202000005.
  * (12) Stark, A. _et al._ Clock transition by continuous dynamical decoupling of a three-level system.  _Scientific reports_ 8, 1–8 (2018).  URL https://www.nature.com/articles/s41598-018-31984-4.
  * (13) Miao, K. C. _et al._ Universal coherence protection in a solid-state spin qubit.  _Science_ 369, 1493–1497 (2020).  URL https://science.sciencemag.org/content/369/6510/1493.
  * (14) Hemmer, P.  Multiplicative suppression of decoherence.  _Science_ 369, 1432–1433 (2020).  URL https://science.sciencemag.org/content/369/6510/1432.
  * (15) Probst, S. _et al._ Hyperfine spectroscopy in a quantum-limited spectrometer.  _Magn. Reson._ 1, 315–330 (2020).  URL https://doi.org/10.5194/mr-1-315-2020.
  * (16) Onizhuk, M. _et al._ Probing the coherence of solid-state qubits at avoided crossings.  _PRX Quantum_ 2, 010311 (2021).  URL https://link.aps.org/doi/10.1103/PRXQuantum.2.010311.
  * (17) Shiddiq, M. _et al._ Enhancing coherence in molecular spin qubits via atomic clock transitions.  _Nature_ 531, 348–351 (2016).  URL https://doi.org/10.1038/nature16984.
  * (18) Liu, J. _et al._ Quantum coherent spin–electric control in a molecular nanomagnet at clock transitions.  _Nature Physics_ 17, 1205–1209 (2021).  URL https://doi.org/10.1038/s41567-021-01355-4.
  * (19) Harding, R. T. _et al._ Spin resonance clock transition of the endohedral fullerene ${}^{15}\mathrm{N}@{\mathrm{c}}_{60}$.  _Phys. Rev. Lett._ 119, 140801 (2017).  URL https://link.aps.org/doi/10.1103/PhysRevLett.119.140801.
  * (20) Zadrozny, J. M., Gallagher, A. T., Harris, T. D. & Freedman, D. E.  A porous array of clock qubits.  _Journal of the American Chemical Society_ 139, 7089–7094 (2017).  URL https://doi.org/10.1021/jacs.7b03123.
  * (21) Collett, C. A. _et al._ A clock transition in the Cr7Mn molecular nanomagnet.  _Magnetochemistry_ 5, 4 (2019).  URL https://www.mdpi.com/2312-7481/5/1/4.
  * (22) Collett, C. A., Santini, P., Carretta, S. & Friedman, J. R.  Constructing clock-transition-based two-qubit gates from dimers of molecular nanomagnets.  _Phys. Rev. Research_ 2, 032037 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.2.032037.
  * (23) Kundu, K. _et al._ 9.2 GHz clock transition in a Lu(II) molecular spin qubit arising from a 3467 MHz hyperfine interaction.  _Nat. Chem._ 14, 392–397 (2020).  URL https://doi.org/10.1038/s41557-022-00894-4.
  * (24) Giménez-Santamarina, S., Cardona-Serra, S., Clemente-Juan, J. M., Gaita-Ariño, A. & Coronado, E.  Exploiting clock transitions for the chemical design of resilient molecular spin qubits.  _Chem. Sci._ 11, 10718–10728 (2020).  URL http://dx.doi.org/10.1039/D0SC01187H.
  * (25) Schweiger, A. & Jeschke, G.  _Principles of pulse electron paramagnetic resonance_ (Oxford University Press on Demand, 2001).
  * (26) Wolfowicz, G. & Morton, J. J. L.  Pulse Techniques for Quantum Information Processing.  _eMagRes_ 5, 1515–1528 (2016).  URL https://doi.org/10.1002/9780470034590.emrstm1521.
  * (27) Hahn, E. L.  Spin echoes.  _Phys. Rev._ 80, 580–594 (1950).  URL https://link.aps.org/doi/10.1103/PhysRev.80.580.
  * (28) Prokof’ev, N. V. & Stamp, P. C. E.  Theory of the spin bath.  _Reports on Progress in Physics_ 63, 669–726 (2000).  URL https://doi.org/10.1088/0034-4885/63/4/204.
  * (29) Ghosh, S. _et al._ Multi-frequency EPR studies of a mononuclear holmium single-molecule magnet based on the polyoxometalate [Ho(W5O18)2]9-.  _Dalton Trans._ 41, 13697–13704 (2012).  URL http://dx.doi.org/10.1039/C2DT31674A.
  * (30) Vonci, M. _et al._ Magnetic excitations in polyoxotungstate-supported lanthanoid single-molecule magnets: An inelastic neutron scattering and ab initio study.  _Inorganic Chemistry_ 56, 378–394 (2017).  URL https://doi.org/10.1021/acs.inorgchem.6b02312.
  * (31) Goldfarb, D. & Stoll, S. (eds.) _EPR Spectroscopy: Fundamentals and Methods, in eMagRes_ (John Wiley and Sons, 2018).
  * (32) Canarie, E. R. & Stoll, S.  Quantitative Structure-Based Prediction of Electron Spin Decoherence in Organic Radicals.  _J. Phys. Chem. Lett_ 11, 3396–3400 (2020).  URL https://dx.doi.org/10.1021/acs.jpclett.0c00768.
  * (33) Chen, J. _et al._ Decoherence in Molecular Electron Spin Qubits: Insights from Quantum Many-Body Simulations.  _J. Phys. Chem. Lett_ 11, 2074–2078 (2020).  URL https://doi.org/10.1021/acs.jpclett.0c00193.
  * (34) Vega, S. & Pines, A.  Operator formalism for double quantum NMR.  _J. Chem. Phys._ 66, 5624-5644 (1977).  URL https://doi.org/10.1063/1.433884.
  * (35) Kragskow, J. G. C. _et al._ Analysis of vibronic coupling in a 4f molecular magnet with FIRMS.  _Nature Communications_ 13, 825 (2022).  URL https://doi.org/10.1038/s41467-022-28352-2.
  * (36) Ullah, A. _et al._ Electrical two-qubit gates within a pair of clock-qubit magnetic molecules.  _arXiv_ (2022, unpublished).  URL https://arxiv.org/abs/2204.09592.
  * (37) Escalera-Moreno, L., Gaita-Ariño, A. & Coronado, E.  Decoherence from dipolar interspin interactions in molecular spin qubits.  _Phys. Rev. B_ 100, 064405 (2019).  URL https://link.aps.org/doi/10.1103/PhysRevB.100.064405.
  * (38) AlDamen, M. A. _et al._ Mononuclear lanthanide single molecule magnets based on the polyoxometalates [Ln(W5O18)2]9– and [Ln($\mathrm{\beta}_{2}$-SiW11O39)2]13– (LnIII = Tb, Dy, Ho, Er, Tm, and Yb).  _Inorganic Chemistry_ 48, 3467–3479 (2009).  URL https://doi.org/10.1021/ic801630z.

## Acknowledgements

We are grateful to Eugenio Coronado and Shimon Vega for insightful discussion,
and we thank Haechan Park and James Fry for assistance in estimating the
closest Ho-1H distances. The spectroscopic and theoretical work reported in
this paper was supported by the Center for Molecular Magnetic Quantum
Materials (M2QM), an Energy Frontier Research Center funded by the US
Department of Energy, Office of Science, Basic Energy Sciences under Award DE-
SC0019330. Experimental work performed at the National High Magnetic Field
Laboratory is supported in part by the National Science Foundation (under
DMR-1644779) and the State of Florida. Synthesis of the HoW10 sample was
supported by: the EU (ERC-2018-AdG-788222 MOL-2D, the QUANTERA project SUMO,
and FET-OPEN grant 862893 FATMOLS); the Spanish MCIU (grant CTQ2017-89993 and
PGC2018-099568-B-I00 co-financed by FEDER, grant MAT2017-89528; the Unit of
excellence ‘Maríade Maeztu’ CEX2019-000919-M); and the Generalitat Valenciana
(Prometeo Program of Excellence).

## Author information

###  0.3 Contributions

K.K. and J.C. performed the simulations. S. Hoffman developed the theory. The
manuscript was written by J.C., J.M. and S. Hill with contributions from S.
Hoffman, K.K., A.G. and X.-G.Z. The experiments were conceived and designed by
S. Hill and D. K., while D.K. performed the measurements. Y.D. prepared the
samples. J.S., A.G., X.-G.Z., S. Hill and H.-P.C. supervised the research.

###  0.4 Corresponding authors

Correspondence and requests for materials should be addressed to S. H., S. H.,
and H.-P. C. (email<EMAIL_ADDRESS><EMAIL_ADDRESS>and
hping@ufl.edu).

## Ethics declaration

The authors declare that they have no competing financial interests.

Figure 1: (a) Zeeman diagram according to the Hamiltonian of Eq. (1), with the
parameters given in the main text. An avoided crossing (a CT) between the two
lowest lying states (blue curves - see labeling) is seen at
$B_{0z}=B_{\mathrm{min}}=23.5$ mT. (b) ESR frequency, $f$, corresponding to
the (clock) transition between the $m_{S}=\pm 1$ states in (a), and the
associated effective gyromagnetic ratio,
$\gamma_{e}^{\textrm{eff}}=df/dB_{0}$. Note that the ESR frequency couples
linearly to $B_{0}$ far from the CT and quadratically at the CT, such that
$\gamma_{e}^{\textrm{eff}}$ crosses through zero. Figure 2: (a) HoW10 ESE
decay curves recorded at 9.18 GHz and 5 K as a function of detuning,
$B_{0}-B_{\mathrm{min}}$; the white dash curve is a fit to a mono-exponential
decay, from which the optimum $T_{\mathrm{m}}=8.43(6)~{}\mu$s is deduced. (b)
FFTs of the decay curves in (a); prominent peaks in the ESEEM spectra are
marked with red, green and blue circles. (c) Plot of ESEEM frequencies in (b)
versus $B_{0}$; the dashed lines correspond to harmonics of the proton Larmor
frequency, the open circles are colored according to the same scheme as those
in (b), and the vertical red line marks the CT.  Figure 3: FFTs of Hahn echo
simulations for the simple case of a single proton ($N=1$) coupled to the
central electron spin (see text for employed parameters) for different
detunings, $B_{0}-B_{\mathrm{min}}$ = +5 mT (a), +20 mT (b) and +50 mT (c);
the inset to (b) shows a representative ESE intensity time trace. Several
relevant frequencies are labeled in the FFT spectra. The inset to (a) plots
$A^{\mathrm{eff}}$ deduced from the first FFT peak versus
$B_{0}-B_{\mathrm{min}}$; the red curve is a simple fit that assumes
$A^{\mathrm{eff}}\propto df/dB_{0}$ from Fig. 1(b). Figure 4: (a) Schematic of
a central spin coupled to $N=7$ nuclear spins, with $\langle A\rangle=8$ MHz.
(b) Simulated ESEEM time traces as a function of detuning,
$B_{0}-B_{\mathrm{min}}$. (c) FFTs of the time traces in (b). (d) Centers of
the main FFT peaks in (c) as a function of $B_{0}$; the circles in (c) and (d)
are colored according to the same scheme as Fig. 2.