# TEAL: New Selection Strategy for Small Buffers
in Experience Replay Class Incremenal Learning

Shahar Shaul-Ariel$\dagger$, Daphna Weinshall$\dagger$
School of Computer Science & Engineering$\dagger$
The Hebrew University of Jerusalem
Jerusalem 91904, Israel
<EMAIL_ADDRESS>

###### Abstract

Continual Learning is an unresolved challenge, whose relevance increases when
considering modern applications. Unlike the human brain, trained deep neural
networks suffer from a phenomenon called Catastrophic Forgetting, where they
progressively lose previously acquired knowledge upon learning new tasks. To
mitigate this problem, numerous methods have been developed, many relying on
replaying past exemplars during new task training. However, as the memory
allocated for replay decreases, the effectiveness of these approaches
diminishes. On the other hand, maintaining a large memory for the purpose of
replay is inefficient and often impractical. Here we introduce _TEAL_ , a
novel approach to populate the memory with exemplars, that can be integrated
with various experience-replay methods and significantly enhance their
performance on small memory buffers. We show that _TEAL_ improves the average
accuracy of the SOTA method XDER as well as ER and ER-ACE on several image
recognition benchmarks, with a small memory buffer of 1-3 exemplars per class
in the final task. This confirms the hypothesis that when memory is scarce, it
is best to prioritize the most typical data.

##  1 Introduction

With the recent advances in deep neural networks, research interest in
incremental learning has increased. Here, the problem of _catastrophic
forgetting_ , as defined by McCloskey and Cohen [20], can be severe. The
necessity to incorporate new task knowledge into an already trained network
has become more relevant given that training over extensive datasets is time-
consuming. It is often impractical to retrain the network from scratch on the
combined original and new task data. Furthermore, access to the already-
trained data might be limited or non-existent.

Figure 1: Catastrophic forgetting when training CIFAR-10 [13] sequentially
with only finetuning.

Various methods have been developed to tackle this challenge across different
frameworks, categorized by Van de Ven et al. [25] into three types: task-
incremental, domain-incremental, and class-incremental learning. The core
concept of incremental learning is that a model must learn tasks sequentially,
one after another. Among these, class-incremental learning is considered the
most challenging. In this scenario, each task introduces new classes, and the
model must identify the specific class of each input without access to the
task ID on which it was trained.

There are several approaches to incremental learning, each with its own
distinct assumptions and configurations. In this paper, we adhere to the
constrained setup described by De Lange et al. [7], which does not rely on
task boundaries during training or testing. This setup maintains a constant
memory size throughout incremental training, ensuring that it does not exceed
a predetermined limit established at the outset. At each incremental step, we
can add new exemplars to the memory buffer only after space has been allocated
for them.

We focus our attention on a prevalent and rather successful framework called
Experience Replay (ER), which involves storing a set of exemplars in memory
and reusing them for rehearsal purposes while training on new tasks. Within
this framework, several strategies exist for selecting which exemplars to
retain in memory. Masana et al. [19] demonstrate that the most successful
strategies are either random-sampling or Herding as defined by Welling [27].
Not surprisingly, the smaller the memory buffer is, the less effective the
strategy is at mitigating catastrophic forgetting. This leads us to the
following question: Are these strategies necessarily optimal for all memory
sizes? In other words, could different strategies prove to be more effective
depending on the size of the memory buffer?

Our proposed method, _TEAL_ , operates on similar principles, aiming to
identify a set of representative exemplars that also exhibit diversity. An
exemplar is deemed representative if its likelihood, when considered within
the distribution of all points in an appropriate latent space, is high.
Leveraging clustering ensures the selection of a diverse set of exemplars.
This approach enhances the model’s capacity to retain previous knowledge while
effectively learning new tasks, even when constrained by a small memory
buffer. Central to the success of this approach is our ability to derive an
appropriate latent space. The subsequent sections elaborate on our
methodology, experimental setup, comparative performance against alternative
selection strategies, and the enhancements TEAL offers to various existing
class-incremental methods.

Related work Several approaches exist for incremental learning [see 7].
Alongside ER, alternatives include Generative Replay (GR) [23], Parameter
Isolation, and Regularization. Similar to ER, GR replays data from previous
tasks during new task training, but uses a generative model to create new
samples instead of retaining exemplars seen by the model [6, 8, 9]. Parameter
Isolation assigns distinct parameters to each task to reduce forgetting [18]
by fixing parameters assigned to previous tasks, while Regularization-based
methods [15] incorporate additional terms into the loss function to retain
prior knowledge while learning from new data.

These methods differ in how they utilize memory: ER stores exemplars, GR
stores a generative model, and Parameter Isolation stores task-specific
parameters. Regularization-based methods do not rely on memory. This diversity
makes it challenging to compare methods directly. In particular, we note that
since generative models tend to be very large, this approach is hardly
suitable to the domain of small memory buffers addressed below. Hence, our
focus in this paper is on ER methods only.

Some CIL methods use an online learning setting, where each training example
is seen only once. In this paper, we primarily use offline methods. Note that
while few-shot incremental learning [24] may appear similar to our work on
class-incremental learning (CIL) with a small memory buffer, there is a
significant difference. We assume that the data for new tasks is initially
sufficient for effective learning, whereas few-shot incremental learning
inherently deals with a scarcity of labeled samples.

Summary of contribution We present _TEAL_ , a novel strategy designed to
select effective exemplars for small memory buffers. _TEAL_ can be seamlessly
integrated with replay-based class-incremental learning methods, significantly
enhancing their performance.

##  2 Problem formulation

CIL involves learning from a sequence of disjoint tasks, where each task
requires the classification of a set of classes, and every class belongs
exclusively to a single task. At each incremental step, the learner receives
the data of the new task and, in some cases, has access to a memory buffer
containing exemplars from previous tasks. After each step, the learner is
tested on the successful classification of all the classes it has seen,
including those seen in previous steps.

More formally [see 19], a class-incremental learning problem $\mathcal{T}$
consists of a sequence of $T$ tasks. Each task $t$ contains labeled samples
from $n_{t}$ different classes $X^{1},\dots,X^{n_{t}}$, where
$X^{i}=\\{(x_{1},i),\dots,(x_{m^{i}},i)\\}$. The tasks do not share classes,
i.e., $X^{i}\cap X^{j}=\emptyset~{}\forall i\neq j$. We denote by $N^{t}$ the
total number of classes in all tasks up to task $t$:
$N^{t}=\sum_{i=1}^{t}{n_{i}}$.

We consider the scenario where there is access to a memory buffer
${\mathcal{M}}$, which has a fixed size throughout the learning of
$\mathcal{T}$. An incremental learner is a learning model (we consider only
deep neural networks) that is trained sequentially on the tasks. Accordingly,
the training on task $t$ is performed on data
$X^{1},\dots,X^{n_{t}},{\mathcal{M}}$, where ${\mathcal{M}}$ contains
exemplars from tasks $1,\dots,t-1$.

Different CIL methods may employ different loss functions, but at test time,
they are all required to classify a given example into its predicted class
while considering all previously seen classes.

##  3 Our method: _TEAL_

In active learning, it has been shown (empirically and formally) that when the
number of labeled examples is small, it is best to choose the most typical
examples for training. In the context of CIL and when the memory buffer is too
small to truly represent the distribution of each class, we propose to adopt
this strategy for the selection of points that are intended to populate the
replay memory buffer. In other words, when the buffer is considerably small,
it should contain representative exemplars and thus retain a more significant
fraction of the previously acquired knowledge.

Our proposed method _TEAL_ , Typicality Election Approach to continual
Learning, aims to address only one component of the general Incremental
Learning (IL) problem, namely, how to populate the memory buffer at the end of
each IL iteration. Accordingly, after task $t$, _TEAL_ should select for each
seen class $n=\frac{|{\mathcal{M}}|}{N_{t}}$ exemplars to populate memory
buffer ${\mathcal{M}}$.

Ideally, any new mechanism for memory population can be combined with any
competitive Experience-Replay-Incremental-Learning (ER-IL) method in which
this mechanism is a separate module. In accordance, we evaluate _TEAL_ by
considering alternative ER-IL methods, where we replace their native mechanism
for buffer population by our proposed method. In almost all cases, this
mechanism involves random sampling from the known labeled set, in order to
faithfully represent the class distribution. In our evaluation study, we also
consider an alternative population mechanism called Herding [22, 12].

###  3.1 Description of method

We begin with a definition of point _typicality_.

######  Definition 3.1 (Typicality).

For each exemplar $x$, let

$Typicality(x)=\left(\frac{1}{K}\sum_{x_{i}\in
K\textrm{-}NN(x)}{\|x-x_{i}\|_{2}}\right)^{-1}$

where $K$ is a fixed number of nearest neighbors111We use $K=20$, but other
options yield similar results..

The essence of _TEAL_ is to populate the memory buffer with points that are
both _typical_ (as captured by the definition above) and _diverse_.

After training on task $t$, the incremental learner has to update the memory
buffer ${\mathcal{M}}$ with data from the new classes
$\\{X^{N^{t-1}+1},\dots,X^{N^{t}}\\}$. Given the fixed size of
${\mathcal{M}}$, the learner must first reduce the number of exemplars already
in the buffer to accommodate new ones from the new classes. For each class,
this involves repeatedly removing examples from the buffer after each IL step.

To manage this, _TEAL_ maintains a list of the selected exemplars from each
class ordered by typicality, reflecting the order in which they are to be
deleted from ${\mathcal{M}}$. This allows the learner to remove the least
typical points in the buffer, whenever new classes appear. To achieve this
goal, whenever a new class is introduced, _TEAL_ selects the initial set of
representative exemplars through an iterative process. The selected list is
ordered such that the most typical exemplars, which are to be retained for as
long as possible, are positioned first, while those that will be removed first
are positioned last.

More specifically, in each repeat of this iterative selection process and for
each class separately, _TEAL_ selects a small number of exemplars (a fraction
of $n$, where $n$ denotes the number of exemplars assigned to the respective
class). Let $S_{i}$ denote the set of exemplars after the $i^{th}$ iteration,
for a total of $k$ iterations. Note that $S_{1}\subseteq
S_{2}\subseteq\dots\subseteq S_{k}$ and that $|S_{k}|=n$. The sizes
$|S_{1}|,\dots,|S_{k}|$ define the _selection pace_ , indicating the pace of
selecting exemplars from new class data.

To initialize the selection process, we need a suitable embedding space for
the class exemplars. To this end we train a deep model on all the available
training data, and use activations in its penultimate layer as a
representation for the new classes [cf. 22]. Subsequently, the selection
process in the $i^{th}$ iteration involves 2 steps (see pseudocode in Alg. 1):

Step 1: Clustering

Since typical examples tend to cluster in bunches [10], and in order to ensure
diversity, we seek typical exemplars from different regions of the embedding
space. We achieve this by dividing the set of labeled points into $|S_{i}|$
clusters, where $|S_{i}|$ represents the number of exemplars that should be
selected from this class after the $i^{th}$ iteration.

Step 2: Typicality

We then select the most typical point from each of the $|S_{i}|-|S_{i-1}|$
largest uncovered clusters, where an uncovered cluster is a cluster from which
no point has already been selected. We order the selected points by their
typicality score.

Input: $X^{c}$ COMMENT# a set of exemplars from class $c$

Input:  $n$ COMMENT# number of examples to choose

Require: $\phi:{\mathcal{X}}\to{\mathbb{R}}^{d}$ COMMENT# current feature
function

Require: $s_{1},\dots,s_{k}$ COMMENT# iterations pace, $s_{1}\leq\dots\leq
s_{k}=n$

Output: ${\mathcal{M}}_{c}$ COMMENT# a set of n exemplars chosen from $X^{c}$

$X^{c}_{emb}\leftarrow\phi(X^{c})$

${\mathcal{M}}_{c}\leftarrow\emptyset$

for all $i=1,\dots,k$ do

$C_{1},\dots,C_{s_{i}}\leftarrow$ clustering_algorithm$(X^{c}_{emb},\ s_{i})$
COMMENT# $|C_{1}|\geq\dots\geq|C_{s_{i}}|$

for all $j=1,\dots,s_{i}-s_{i-1}$ do

add $\underset{x\in C_{j}}{\mathrm{argmax}}\,\\{Typicality(x)\\}$ to
${\mathcal{M}}_{c}$

end for

end for

return ${\mathcal{M}}_{c}$

Algorithm 1 _TEAL_ ConstructExemplarSet

###  3.2 Using _TEAL_ for incremental learning

The integration of _TEAL_ into a general class-incremental algorithm is
described in Algs. 2-3.

Input: ${\mathcal{M}}_{c}=(m_{1},\dots,m_{|{\mathcal{M}}_{c}|})$ COMMENT#
current exemplars set

Input:  $n$ COMMENT# number of examples to keep

return $(m_{1},\dots,m_{n})$

Algorithm 2 ReduceExemplarSet

Input: $X^{t}=\\{X^{N^{t-1}+1},\dots,X^{N^{t}}\\}$ COMMENT# training examples
in per-class sets

Require: $\theta^{t-1}$ COMMENT# current model parameters

Require: ${\mathcal{M}}=({\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N^{t-1}})$
COMMENT# current exemplar sets

$\theta^{t}\leftarrow\textrm{training\\_model}(X^{t},{\mathcal{M}};\theta^{t-1})$

$n\leftarrow\frac{|{\mathcal{M}}|}{N^{t}}$ COMMENT# number of exemplars per
class

for $c\in\\{1,\dots,N^{t-1}\\}$ do

${\mathcal{M}}_{c}\leftarrow ReduceExemplarSet({\mathcal{M}}_{c},n)$

end for

for $c\in\\{N^{t-1}+1,\dots,N^{t}\\}$ do

${\mathcal{M}}_{c}\leftarrow ConstructExemplarSet(X^{c},n,\theta^{t})$

end for

${\mathcal{M}}\leftarrow({\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N^{t}})$

Algorithm 3 IncrementalTrain

##  4 Empirical evaluation

We report two settings:  (i) _Integrated:_ we evaluate the beneficial
contribution of _TEAL_ to competitive ER-IL methods, by replacing their native
selection strategy with _TEAL_ (Section 4.2.1);  (ii) _Alternative selection
strategies:_ we evaluate the beneficial contribution of _TEAL_ as compared to
Herding (Section 4.2.2).  In both cases, we incrementally train a deep-
learning model with a fixed-size replay buffer and monitor the average
accuracy defined below upon completion of each task $t$.

###  4.1 Methodology

The majority of the experiments are conducted using the open-source CL library
Avalanche [16], with the exception of those involving XDER [1], which is not
integrated into Avalanche. For XDER, we utilize the code provided by the
authors. In order to guarantee a fair comparison and in all experiments, we
employ identical network architectures and maintain consistent experimental
conditions.

We make certain in both settings that the buffer remains class-balanced,
namely, each update maintains an equal representation of exemplars across all
classes. When the buffer size is not evenly divisible by the number of
classes, there may be instances where certain classes possess one additional
exemplar. To maintain this balance, we follow this procedure: Prior to adding
new exemplars we calculate $n=\frac{|{\mathcal{M}}|}{N^{t}}$, where $N^{t}$ is
the total number of existing and new classes. Subsequently, we adjust the
number of exemplars from existing classes to accommodate $n$ by removing the
least likely points, followed by the addition of $n$ exemplars from each new
class.

In the second setting, we employ a simple baseline ER-IL model. This model
updates a fixed-size buffer of exemplars after training each task $t$ and
replays it while training on task $t+1$. Additionally, we utilize a weighted
data loader to ensure a balanced mix of data from both the new classes and the
exemplars stored in the buffer in each batch. Across experiments, the only
variation lies in the selection strategy employed: some experiments utilize
random-sampling, others utilize Herding, and the remaining experiments
leverage _TEAL_. In some experiments, we replace the native buffer selection
strategy of an existing ER-IL method with _TEAL_ or Herding to assess their
impact. Other than this change in the method’s buffer population mechanism,
everything else remains the same.

Datasets We use several well known Continual Learning benchmarks.  (i) Split
CIFAR-100[22, 5], a dataset created by splitting CIFAR-100 [13] into 10 tasks,
each containing 10 different classes of 32X32 images, with 500 images per
class for training and 100 for testing.  (ii) Split tinyImageNet, a dataset
created by splitting tinyImageNet [14] into 10 tasks, each containing 20
different classes of 64X64 images, with 500 images per class for training and
50 for testing.  (iii) Split CUB-200[3], a dataset created by splitting the
CUB-200 [26] high-resolution image classification dataset, consisting of 200
categories of birds, with around 30 images per class for training and 30 for
testing.

Architectures In the first setting, except for one experiment, we employ a
ResNet-18 model in all our experiments. The exception is the experiment
involving the training of XDER on the Split CUB-200 dataset. Due to
computational limitations, we used instead a pre-trained ResNet-50 model [11],
and this was maintained under all the relevant conditions. In the second
setting, we use a smaller version of ResNet-18 [11] as a simple baseline ER-IL
model [see 17].

Metrics As customary, we employ a score that reflects the accuracy during each
stage of the incremental training. Let $T$ denote the total number of tasks,
and $a_{t,i}$ denote the accuracy of task $i$ at the end of task $t$ ( $i\leq
t\leq T$). For every task $t={1,\dots,T}$, define it average accuracy
$A_{t}=\frac{1}{t}\sum_{i=1}^{t}a_{t,i}$. This metric provides a single value
for each incremental step, enabling us to directly compare different methods
at each step. Since all of the experiments we conduct are task-balanced in the
sense that each task has an equal number of classes, there is no need to add a
weighting term to the metric.

Baselines The following methods are suitable for our investigation as they are
designed to address ER-IL: XDER [1] updates the memory buffer by adding
current information to past memories and preparing the model for new tasks.
ER-ACE [2], or ER with Asymmetric Cross-Entropy, uses a separate loss for new
and past tasks. BiC [28] introduces a bias correction layer to address class
imbalance in incremental learning. iCaRL [22] combines the incremental
classifier and representation learning with a nearest-mean-of-exemplars
strategy, and Herding for exemplar selection. Since the classifier is based on
the nearest-mean-of-exemplars, Herding is crucial for the selection strategy,
as it guarantees that the features of the exemplars in the buffer approximate
the original classes means. GEM [17] uses Gradient Episodic Memory to prevent
gradient leakage between new and past tasks through constrained optimization.
GDumb [21] employs a greedy sampling strategy for memory storage and retrains
the neural network from scratch during inference. ER [4] explores different
selection strategies in a simple ER framework.

In order to evaluate the competitiveness of these methods, we conducted
multiple runs on Split CIFAR-100, utilizing fixed buffer sizes of 300, 500,
and 2000. Fig. 2 shows the results.

(a) $|{\mathcal{M}}|=300$

(b) $|{\mathcal{M}}|=500$

(c) $|{\mathcal{M}}|=2000$

Figure 2: Baseline CIL methods comparison on Split CIFAR-100 with different
fixed buffer sizes. For each method, each point represents the average
accuracy achieved on the corresponding task $t$, $A_{t}$.

Probably because the memory buffer is so small, and as can be seen in Fig. 2,
GEM and GDumb perform poorly in a relative comparison. BiC performs well on
buffer size 2000, but significantly degrades for smaller buffers for which it
is not suited. iCarl demonstrates high performance, but it is highly non-
modular, being heavily dependent on the Herding selection strategy while using
a K-Means classifier. Hence we cannot integrate it with _TEAL_.

We are therefore left with methods XDER, ER-ACE and ER, which are both
suitable and competitive, to be examined with and without _TEAL_ as the
selection mechanism. For ER-ACE, since it depends on populating the buffer
with exemplars at the beginning of the task before training on the new
classes, we cannot rely on the representation provided by a trained model.
Therefore, we adjust the integration of _TEAL_ by filling the buffer in two
stages: First, at the beginning of task $t$, we update the buffer with
exemplars from the new classes using the random-sampling selection strategy.
Second, after training task $t$, we replace the exemplars from the new classes
in the buffer with exemplars chosen using _TEAL_. It’s important to note that
throughout this process, the buffer maintains the fixed size of exemplars it
allowed, ensuring compliance with the class-incremental framework.

###  4.2 Main results

(b) $|{\mathcal{M}}|=500$ (c) $|{\mathcal{M}}|=2000$

Split CIFAR-100

(a) $|{\mathcal{M}}|=300$ (e) $|{\mathcal{M}}|=1000$ (f)
$|{\mathcal{M}}|=2000$

Split tinyImageNet

(d) $|{\mathcal{M}}|=400$ (h) ER (i) XDER

Split CUB-200

(g) ER-ACE

Figure 3: The performance of ER-IL methods with and without _TEAL_. Each plot
displays the average accuracy after training incrementally on a different
number of classes. Each color corresponds to a different ER-IL method, where
the continuous line represents the vanilla method, while the dashed line
represents the method with _TEAL_ as its selection strategy. The error bars
correspond to standard error based on 4-10 repetitions.

####  4.2.1 _TEAL_ integrated with SOTA ER-IL methods

Split CIFAR-100

(a) $|{\mathcal{M}}|=300$

(b) $|{\mathcal{M}}|=500$

(c) $|{\mathcal{M}}|=2000$ (e) $|{\mathcal{M}}|=1000$ (f)
$|{\mathcal{M}}|=2000$

Split tinyImageNet

(d) $|{\mathcal{M}}|=400$ (h) ER (i) XDER

Split CUB-200

(g) ER-ACE

Figure 4: Performance improvement of _TEAL_ when integrated with different ER-IL methods. Each graph depicts the difference in accuracy between _TEAL_ and another method (XDER, ER-ACE, and ER) across all tasks. The first two rows present the three methods improvement on Split CIFAR-100 and Split tinyImageNet with different fixed buffer sizes denoted as $|{\mathcal{M}}|$. The last row pertains to Split CUB-200, where each graph showcases different buffer sizes for a different method. Table 1: Final average accuracy for method with and without _TEAL_ Dataset | Split CIFAR-100 | Split tinyImageNet | Split CUB-200
---|---|---|---
$|{\mathcal{M}}|$ | 300 | 500 | 2000 | 400 | 1000 | 2000 | 200 | 400
XDER |  41.97±0.25 |  47.97±0.22 |  57.69±0.21 |  20.49±0.19 |  32.68±0.24 |  39.76±0.31 | – |  43.63±1.78
XDER+_TEAL_ |  45.05±0.24 |  50.29±0.2 |  58.79±0.17 |  23.61±0.32 |  34.28±0.24 |  40.76±0.23 | – |  49.55±0.47
Improvement | 3.08 | 2.32 | 1.1 | 3.13 | 1.59 | 1.0 | – | 5.96
ER-ACE |  35.63±0.36 |  40.6±0.34 |  53.6±0.17 |  13.02±0.2 |  17.37±0.34 |  21.71(±0.17) |  8.56±0.15 |  11.25±0.99
ER-ACE+_TEAL_ |  40.59±0.2 |  44.25±0.35 |  56.11±0.19 |  14.43±0.29 |  18.88±0.21 |  24.12±0.24 |  10.41±0.16 |  12.33±0.26
Improvement | 4.97 | 3.65 | 2.51 | 1.41 | 1.51 | 2.4 | 1.85 | 1.08
ER |  14.82±0.21 |  18.21±0.29 |  36.54±0.31 |  7.75±0.06 |  7.7±0.06 |  8.27±0.08 |  8.98±0.18 |  12.01±0.52
ER+_TEAL_ |  17.97±0.25 |  22.39±0.31 |  38.17±0.37 |  7.86±0.05 |  7.99±0.06 |  8.82±0.09 |  10.37±0.23 |  14.13.76±0.3
Improvement | 3.15 | 4.18 | 1.62 | 0.11 | 0.29 | 0.56 | 1.39 | 2.12

As mentioned above, we investigate the baseline methods XDER, ER-ACE and ER,
comparing their performance when using either their native selection
strategies or _TEAL_. We examined Split CIFAR-100 with buffer sizes of 300,
500 and 2000, Split tinyImageNet with buffer sizes of 400, 1000 and 2000, and
Split CUB-200 with buffer sizes of 200222Due to computational limitation of
the relevant package, we did not run XDER with buffer size 200. and 400. Fig.
3 shows the average accuracy $A_{t}$ after each task $t$. Fig. 4 shows the
corresponding improvement of _TEAL_. Finally, the results of the final average
accuracy $A_{t}$ and the differences ’improved-vanilla’ are reported in Table
1.

As demonstrated in Fig. 4, in most cases, the improvement achieved by _TEAL_
increases as incremental training progresses, resulting in a more significant
improvement by the final task. Indeed, we do not expect significant
differences immediately after the first task, as catastrophic forgetting
hasn’t occurred yet. Reassuringly, for the majority of the experiments, the
improvement becomes apparent starting from the second task.

Table 1 shows that the difference (’improvement’ row) is always positive,
namely, _TEAL_ always improves the outcome. Looking at the improvement row for
the XDER method, it is evident that the improvement provided by _TEAL_
increases as the buffer size decreases. This trend is also demonstrated for
ER-ACE and partially for ER on Split CIFAR-100.

####  4.2.2 Comparison of selection strategy

We begin with a simple ER-IL model on Split CIFAR-100333The order of the
classes and the partition to tasks is randomly selected and fixed for all
comparisons. with various buffer sizes, and three different selection
strategies: random, Herding and _TEAL_. As depicted in Fig. 5, our method
enhances performance as compared to random selection by almost 3% for smaller
buffers (300, 400, 500) and about 2.5% for larger buffers (1000, 2000). When
compared to Herding, our method is more effective for smaller buffer sizes,
enhancing performance by up to 1.2% for a buffer size of 300, but less
effective with large buffers. Similar results, using different class ordering,
are shown in App A.1.

Figure 5: Split CIFAR-100: performance difference between _TEAL_ and either
random-sampling or Herding selection. Each bar corresponds to a fixed buffer
size.

Dataset | Split CIFAR-100
---|---
$|{\mathcal{M}}|$ | 300 | 500 | 2000
XDER+_TEAL_ __OneTime_ |  44.13±0.26 |  49.87±0.22 |  59.05±0.14
XDER+_TEAL_ | 45.05±0.24 | 50.29±0.2 |  58.79±0.17
ER-ACE+_TEAL_ __OneTime_ |  35.81±0.42 |  40.16±0.3 |  51.52±0.24
ER-ACE+_TEAL_ | 36.4±0.12 |  40.26±0.23 |  51.68±0.35
ER+_TEAL_ __OneTime_ |  17.64±0.12 |  22.01±0.18 | 40.65±0.32
ER+_TEAL_ | 18.1±0.1 |  22.07±0.13 |  39.96±0.23

Table 2: The improvement of 2 variants of _TEAL_ when integrated with
different CIL methods: XDER, ER-ACE and ER. In almost all cases, the original
_TEAL_ matches or outperforms the OneTime variant.

Given the results shown in Fig. 5, wherein _TEAL_ seems less effective than
herding with large buffers when training the simple baseline model, we repeat
this comparison while integrating both methods - _TEAL_ and Herding - into
competitive ER-IL methods. Fig. 6 shows the difference between the improvement
obtained by _TEAL_ and the one obtained by Herding. Since in almost all cases
this difference is positive, we may conclude that _TEAL_ is significantly more
beneficial than Herding.

(a) $|{\mathcal{M}}|=300$

(b) $|{\mathcal{M}}|=500$

(c) $|{\mathcal{M}}|=2000$

Figure 6: The difference between the improvement obtained by _TEAL_ and the
one obtained by Herding, showing 2 ER-IL methods and 3 buffer sizes while
training on the Split CIFAR-100 dataset.

###  4.3 Ablation study

We explore an alternative variant of _TEAL_ , which selects the set of
exemplars in a single pass. We call this variant _TEAL_ _OneTime. While still
selecting exemplars for each class separately, this variant first partitions
the exemplars of each class into $n=\frac{|{\mathcal{M}}|}{N^{t}}$ clusters,
and then selects the most typical exemplar from each cluster in descending
order of cluster size, preserving the order of selections.

The experiments on the two variants are conducted using Split CIFAR-100 with
the smaller version of ResNet-18 mentioned in Section 4.1 and buffer sizes of
300, 500 and 2000. We investigate the performance enhancement of each variant
when integrated with XDER, ER-ACE, and ER. The results are shown in Table 2.
Whenever there is a significant difference in performance, the iterative
variant outperforms the other one in 4 out of 5 cases. However, the
differences are not substantial.

##  5 Summary and discussion

We propose a new mechanism to select exemplars for the memory buffer in
replay-based CIL methods. This method is based on the principles of diversity
and representativeness. When the memory buffer is relatively small, our method
_TEAL_ is shown to outperform both the native mechanism of each ER-IL method
(usually random selection) and an alternative selection mechanism called
Herding. Even when the buffer is large, our method is beneficial in almost all
studied cases.

Our future work will investigate ways to determine at which point, if any, the
added value of our method diminishes or possibly becomes harmful. We will also
integrate _TEAL_ with other replay-based CIL methods that utilize different
selection strategies, and analyze the impact on their performance.

##### Acknowledgement

This work was supported by grants from the Israeli Council of Higher Education
and the Gatsby Charitable Foundations.

## References

  * Boschini et al. [2022] Matteo Boschini, Lorenzo Bonicelli, Pietro Buzzega, Angelo Porrello, and Simone Calderara.  Class-incremental continual learning into the extended der-verse.  _IEEE transactions on pattern analysis and machine intelligence_ , 45(5):5497–5512, 2022.
  * Caccia et al. [2021] Lucas Caccia, Rahaf Aljundi, Nader Asadi, Tinne Tuytelaars, Joelle Pineau, and Eugene Belilovsky.  New insights on reducing abrupt representation change in online continual learning.  _arXiv preprint arXiv:2104.05025_ , 2021.
  * Chaudhry et al. [2018] Arslan Chaudhry, Marc’Aurelio Ranzato, Marcus Rohrbach, and Mohamed Elhoseiny.  Efficient lifelong learning with A-GEM.  _CoRR_ , abs/1812.00420, 2018.  URL http://arxiv.org/abs/1812.00420.
  * Chaudhry et al. [2019a] Arslan Chaudhry, Marcus Rohrbach, Mohamed Elhoseiny, Thalaiyasingam Ajanthan, P Dokania, P Torr, and M Ranzato.  Continual learning with tiny episodic memories.  In _Workshop on Multi-Task and Lifelong Reinforcement Learning_ , 2019a.
  * Chaudhry et al. [2019b] Arslan Chaudhry, Marcus Rohrbach, Mohamed Elhoseiny, Thalaiyasingam Ajanthan, Puneet Kumar Dokania, Philip H. S. Torr, and Marc’Aurelio Ranzato.  Continual learning with tiny episodic memories.  _CoRR_ , abs/1902.10486, 2019b.  URL http://arxiv.org/abs/1902.10486.
  * Choi et al. [2021] Yoojin Choi, Mostafa El-Khamy, and Jungwon Lee.  Dual-teacher class-incremental learning with data-free generative replay.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 3543–3552, 2021.
  * De Lange et al. [2021] Matthias De Lange, Rahaf Aljundi, Marc Masana, Sarah Parisot, Xu Jia, Aleš Leonardis, Gregory Slabaugh, and Tinne Tuytelaars.  A continual learning survey: Defying forgetting in classification tasks.  _IEEE transactions on pattern analysis and machine intelligence_ , 44(7):3366–3385, 2021.
  * Gao and Liu [2023] Rui Gao and Weiwei Liu.  Ddgr: Continual learning with deep diffusion-based generative replay.  In _International Conference on Machine Learning_ , pages 10744–10763. PMLR, 2023.
  * Gautam et al. [2024] Chandan Gautam, Sethupathy Parameswaran, Ashish Mishra, and Suresh Sundaram.  Generative replay-based continual zero-shot learning.  In _Towards Human Brain Inspired Lifelong Learning_ , pages 73–100. World Scientific, 2024.
  * Hacohen et al. [2022] Guy Hacohen, Avihu Dekel, and Daphna Weinshall.  Active learning on a budget: Opposite strategies suit high and low budgets.  In _International Conference on Machine Learning_. PMLR, 2022.
  * He et al. [2016] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _2016 IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2016, Las Vegas, NV, USA, June 27-30, 2016_ , pages 770–778. IEEE Computer Society, 2016.  doi: 10.1109/CVPR.2016.90.  URL https://doi.org/10.1109/CVPR.2016.90.
  * Hou et al. [2019] Saihui Hou, Xinyu Pan, Chen Change Loy, Zilei Wang, and Dahua Lin.  Learning a unified classifier incrementally via rebalancing.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 831–839, 2019.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  _Online_ , 2009.
  * Le and Yang [2015] Ya Le and Xuan S. Yang.  Tiny imagenet visual recognition challenge.  2015\.  URL https://api.semanticscholar.org/CorpusID:16664790.
  * Li and Hoiem [2017] Zhizhong Li and Derek Hoiem.  Learning without forgetting.  _IEEE transactions on pattern analysis and machine intelligence_ , 40(12):2935–2947, 2017.
  * Lomonaco et al. [2021] Vincenzo Lomonaco, Lorenzo Pellegrini, Andrea Cossu, Antonio Carta, Gabriele Graffieti, Tyler L Hayes, Matthias De Lange, Marc Masana, Jary Pomponi, Gido M Van de Ven, et al.  Avalanche: an end-to-end library for continual learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 3600–3610, 2021.
  * Lopez-Paz and Ranzato [2017] David Lopez-Paz and Marc’Aurelio Ranzato.  Gradient episodic memory for continual learning.  _Advances in neural information processing systems_ , 30, 2017.
  * Mallya and Lazebnik [2018] Arun Mallya and Svetlana Lazebnik.  Packnet: Adding multiple tasks to a single network by iterative pruning.  In _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , pages 7765–7773, 2018.
  * Masana et al. [2022] Marc Masana, Xialei Liu, Bartłomiej Twardowski, Mikel Menta, Andrew D Bagdanov, and Joost Van De Weijer.  Class-incremental learning: survey and performance evaluation on image classification.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 45(5):5513–5533, 2022.
  * McCloskey and Cohen [1989] Michael McCloskey and Neal J Cohen.  Catastrophic interference in connectionist networks: The sequential learning problem.  In _Psychology of learning and motivation_ , volume 24, pages 109–165. Elsevier, 1989.
  * Prabhu et al. [2020] Ameya Prabhu, Philip HS Torr, and Puneet K Dokania.  Gdumb: A simple approach that questions our progress in continual learning.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part II 16_ , pages 524–540. Springer, 2020.
  * Rebuffi et al. [2017] Sylvestre-Alvise Rebuffi, Alexander Kolesnikov, Georg Sperl, and Christoph H Lampert.  icarl: Incremental classifier and representation learning.  In _Proceedings of the IEEE conference on Computer Vision and Pattern Recognition_ , pages 2001–2010, 2017.
  * Shin et al. [2017] Hanul Shin, Jung Kwon Lee, Jaehong Kim, and Jiwon Kim.  Continual learning with deep generative replay.  _Advances in neural information processing systems_ , 30, 2017.
  * Tian et al. [2024] Songsong Tian, Lusi Li, Weijun Li, Hang Ran, Xin Ning, and Prayag Tiwari.  A survey on few-shot class-incremental learning.  _Neural Networks_ , 169:307–324, 2024.
  * Van de Ven et al. [2022] Gido M Van de Ven, Tinne Tuytelaars, and Andreas S Tolias.  Three types of incremental learning.  _Nature Machine Intelligence_ , 4(12):1185–1197, 2022.
  * Wah et al. [2011] Catherine Wah, Steve Branson, Peter Welinder, Pietro Perona, and Serge Belongie.  The caltech-ucsd birds-200-2011 dataset.  2011\.
  * Welling [2009] Max Welling.  Herding dynamical weights to learn.  In _Proceedings of the 26th annual international conference on machine learning_ , pages 1121–1128, 2009.
  * Wu et al. [2019] Yue Wu, Yinpeng Chen, Lijuan Wang, Yuancheng Ye, Zicheng Liu, Yandong Guo, and Yun Fu.  Large scale incremental learning.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 374–382, 2019.

## Appendix

##  Appendix A Implementation details

The source code used in this study would be uploaded to our GitHub repository
and will be made public upon acceptance.

_TEAL_ , as described in Alg. 1, requires an iterations pace
$s_{1}\leq\dots\leq s_{k}=n$ which indicates the pace of selecting exemplars
from new class data. In both settings we use a logarithmic pace. We set a base
$b=1.4$, and define $s_{1}=\lfloor b^{4}\rfloor,s_{2}=\lfloor
b^{5}\rfloor,\dots,s_{k}=\lfloor b^{log_{b}n}\rfloor$.

###  A.1 Stand-alone setting

For all selection strategies, we use a smaller ResNet-18, as mentioned above,
trained for 200 epochs. Optimization is performed with an SGD optimizer using
Nesterov momentum of 0.9, a weight decay of 0.0002, and a learning rate that
starts at 0.1 and decays by a factor of 0.3 every 66 epochs. Training is
conducted with a batch size of 128 examples, and data augmentation is applied
through random cropping and horizontal flips.

We run this setting on Split CIFAR-100 with three different random orders of
classes:

  1. 1.

The Order we display on Fig. 5 is: [44, 19, 93, 90, 71, 69, 37, 95, 53, 91,
81, 42, 80, 85, 74, 56, 76, 63, 82, 40, 26, 92, 57, 10, 16, 66, 89, 41, 97, 8,
31, 24, 35, 30, 65, 7, 98, 23, 20, 29, 78, 61, 94, 15, 4, 52, 59, 5, 54, 46,
3, 28, 2, 70, 6, 60, 49, 68, 55, 72, 79, 77, 45, 1, 32, 34, 11, 0, 22, 12, 87,
50, 25, 47, 36, 96, 9, 83, 62, 84, 18, 17, 75, 67, 13, 48, 39, 21, 64, 88, 38,
27, 14, 73, 33, 58, 86, 43, 99, 51]

  2. 2.

The order in 7(a) is: [45, 15, 90, 32, 35, 63, 17, 72, 79, 96, 48, 36, 16, 11,
23, 80, 22, 58, 3, 62, 50, 33, 66, 99, 43, 76, 7, 57, 81, 82, 6, 10, 24, 52,
95, 73, 91, 21, 38, 31, 85, 59, 13, 69, 75, 70, 64, 8, 77, 34, 46, 39, 92, 0,
44, 98, 49, 9, 4, 61, 12, 83, 28, 78, 40, 88, 54, 5, 26, 41, 89, 20, 84, 2, 1,
55, 19, 74, 25, 37, 42, 14, 30, 18, 67, 71, 68, 27, 60, 51, 29, 56, 93, 47,
97, 94, 86, 87, 65, 53]

  3. 3.

The order in 7(b) is: [48, 97, 1, 81, 90, 49, 10, 8, 7, 20, 70, 73, 75, 14,
91, 38, 47, 21, 74, 52, 80, 98, 59, 12, 71, 85, 6, 34, 55, 82, 95, 63, 78, 15,
94, 60, 99, 76, 25, 40, 88, 0, 62, 96, 87, 51, 16, 18, 9, 19, 29, 45, 86, 53,
56, 31, 28, 61, 30, 33, 4, 67, 64, 58, 50, 54, 3, 13, 37, 27, 66, 77, 84, 69,
2, 41, 22, 92, 42, 44, 11, 36, 46, 79, 65, 72, 23, 17, 39, 5, 89, 35, 24, 83,
43, 57, 93, 32, 68, 26]

(a) Order 2

(b) Order 3

Figure 7: Enhanced accuracy of _TEAL_ compared to random-sampling and Herding
selection strategies, additional two different classes order.

###  A.2 Integrated setting

Here we train a ResNet-18 model for 100 epochs using the same optimizer as in
A.1, with the same batch size and data augmentations, with some exceptions.
ER-ACE starts with a learning rate of 0.01, and all experiments on Split
CUB-200 are conducted for 30 epochs with a batch size of 16 due to the dataset
size and the resolution of its images.

###  A.3 Baselines setting

The setting for the experiments of the baseline methods is the same as in A.2
with some exceptions. For BiC, the number of training epochs is 250, and the
learning rate scheduler decays the learning rate by a factor of 0.1 on epochs
100, 150 and 200. For GEM, the batch size is 32 and the learning rate starts
from 0.03.For GDumb, the batch size is 32. For iCaRL the learning rate starts
from 2, the weight decay is 0.00001 and the learning rate scheduler decays the
learning rate by a factor of 0.2 on epochs 49 and 63.

##  Appendix B Compute resources

All experiments involved training deep learning models, necessitating the use
of GPUs. For the Split CIFAR-100 experiments, 10 GB of GPU memory was used.
The Split tinyImageNet experiments required 22 GB of GPU memory, while the
Split CUB-200 experiments utilized 45 GB of GPU memory. Any other experiments
conducted for the full research and not reported in the paper required the
same compute resources.


# PULSEE: A software for the quantum simulation of an extensive set of nuclear
magnetic and quadrupole resonance observables

Davide Candoli  Samuele Sanna  Vesna F. Mitrović  Department of Physics and
Astronomy “A. Righi”, University of Bologna, 40127 Bologna, Italy  Department
of Physics, Brown University, Providence, 02912 Rhode Island, USA

###### Abstract

We present an open-source software for the simulation of observables in
nuclear magnetic/quadrupole resonance experiments (NMR/NQR) on solid-state
samples, developed to assist experimental research in the design of new
strategies for the investigation of quantum materials inspired by the early
NMR/NQR quantum computation protocols. The software is based on a quantum
mechanical description of nuclear spin dynamics in NMR/NQR experiments and has
been widely tested on both theoretical and experimental available results.
Moreover, the structure of the software allows an easy generalization of basic
experiments to more sophisticated ones, as it includes all the libraries
required for the numerical simulation of generic spin systems. In order to
make the program easily accessible to a large user base, we developed a user-
friendly graphical interface and fully-detailed documentation. Lastly, we
portray several examples of the execution of the code that demonstrate the
potential of NMR/NQR for the scopes of quantum control and quantum information
processing.

###### keywords:

Nuclear magnetic resonance; Nuclear quadrupole resonance; Quadrupolar
interaction; Nuclear spin dynamics; NMR/NQR quantum computing; Python 3

††journal: Computer Physics Communications

PROGRAM SUMMARY/NEW VERSION PROGRAM SUMMARY

Program Title: PULSEE (Program for the simULation of nuclear Spin Ensemble
Evolution)
CPC Library link to program files: (to be added by Technical Editor)
Developer’s repository link: https://github.com/DavideCandoli/PULSEE
Code Ocean capsule: (to be added by Technical Editor)
Licensing provisions: GPLv3
Programming language: Python 3
Nature of problem: Applications of nuclear magnetic/quadrupole resonance
techniques to study the properties of materials often requires extensive
spectral simulations. On the other hand, applications of magnetic resonance
techniques to quantum information science (QIS) involves different set of
observables. Available simulation software address only one of these
applications, that is, either detailed spectral simulations [1] and/or QIS
relevant observables [2]. Therefore, there is a need for an up-to-date and
easily accessible software for the simulation of an extensive set of nuclear
magnetic/quadrupole resonance experimental observables that allow to reproduce
the behaviour/responce of nuclear systems with varying degree of complexity.
Solution method: The open-source Python code provides an extensive set of
libraries for the simulation of the time evolution of nuclear spins in the
presence of specific interactions and the reproduction of spectra and other
observables measured in nuclear magnetic/quadrupole resonance experiments. The
ready-to-use software features a user-friendly graphical interface.

## References

  * [1] F. A. Perras, C. M. Widdifield, and D. L. Bryce, “QUEST - Quadrupolar Exact Software: A fast graphical program for the exact simulation of NMR and NQR spectra for quadrupolar nuclei,” Solid State Nuclear Magnetic Resonance, vol. 45-46, pp. 36-44, (2012).
  * [2] D. Possa, A. C. Gaudio, and J. C. C. Freitas, “Numerical simulation of NQR/NMR: Applications in quantum computing,” Journal of Magnetic Resonance, vol. 209, pp. 250-260, (2011).

##  1 Introduction

Nuclear magnetic and quadrupole resonance (NMR/NQR) have a long-standing
reputation as accurate methods for the microscopic investigation of materials
based on remarkably simple working principles. In addition to being a dominant
tool in chemistry, materials science, structural biology and medicine, NMR
represents an essential tool in QIS [1]. NMR also allows for fundamental tests
of quantum mechanics [2] and fundamental condensed matter physics, as well as
for probing microscopic spin and charge properties of materials [3, 4, 5, 6].
These features are the reason for the success of magnetic resonance techniques
in implementing one of the first quantum information processors: the high
degree of control of nuclear spins that they provide lends naturally to the
purposes of basic quantum computing, and has made it possible to witness for
the first time the experimental realization of several quantum algorithms [7,
8, 9, 10, 11, 12, 13, 14, 15]. Specifically, the handling of quantum systems
to perform data processing tasks in NMR is accomplished through the
application of specific radio frequency (RF) pulses (logic gates) on ensemble
states adequately prepared, called pseudopure states (PPS) [13, 14, 16, 17,
18]. The logic gates can be executed with high fidelity due to the high level
of control of the quantum evolution of nuclear spins. The true power of a
quantum computing lays in its ability to coherently manipulate entangled
qubits. A 12-qubit NMR based quantum computer holds the record for the largest
quantum computer, _i.e._ high fidelity implementation of a quantum algorithm
with coherent manipulation of 12 qubits [19].

Nonetheless, the long term interest in the applications of NMR and NQR in
quantum computing has faded since they present some major limitations when it
comes to implementing a large scale quantum computer. However, the protocols
developed for the manipulation of NMR/NQR qubits promise to be valuable
resources for the exploration of unknown properties of materials. These
methods may be used to highlight various features of NMR/NQR spectra, and in
this way reveal important details of the underlying nuclear interactions.

Realizing protocols that enhance the power of NMR and NQR through the careful
manipulation of nuclear spin degrees of freedom requires the development of a
software that simulates such experimental techniques that feature the
representation of nuclear spin states. To our knowledge, already existing
programs that simulate NMR/NQR experiments are coded with outdated programming
languages [20, 21], are based on closed-source software systems [22] or lack
the density matrix visualization of spin states [23]: therefore, our purpose
was to write an up-to-date software which could integrate all the features
relevant for the current needs of research, making them fully accessible and
extensively documented.

Our software PULSEE (Program for the simULation of nuclear Spin Ensemble
Evolution) [24] is grounded on the quantum mechanical description of NMR and
NQR, and is able to simulate the time evolution of nuclear spins in solid-
state samples in a wide variety of configurations observed experimentally.
After simulating the acquisition of the characteristic observable measured in
NMR/NQR laboratories, such as for example the free induction decay signal
(FID), one generates the NMR/NQR spectrum in a form which can be put in direct
comparison with results coming from real experiments.

In section 2 we give an overview of the theory of NMR and NQR, including both
the description of nuclear spin dynamics and the generation of the spectra
from the analysis of the FID. In section 3 we present the simulation software,
providing practical information about its installation, structure and usage.
In section 4 we report several examples of simulations carried out with our
software, which have been chosen due to their relevance for the purposes of
quantum control and quantum information processing.

##  2 Theoretical background

Nuclear magnetic and quadrupole resonance (NMR/NQR) involve the time evolution
of resonantly perturbed nuclear spins in matter. Experimentally, the
distinction between the two methods lies in the different nuclear interactions
being probed: NMR pertains to nuclei coupled to a local magnetic field (that
is, an externally applied magnetic field), while NQR deals with the
quadrupolar interaction between each nucleus and the surrounding electronic
charges. From a theoretical point of view, though, it is convenient to treat
the problem where both interactions are simultaneously present, since it
includes all the possible intermediate configurations between pure NMR and
pure NQR [5].

Our work focuses on the NMR study of single crystals, where assuming in first
approximation that the nuclei are identical and non-interacting, each of them
will be described by the same stationary Hamiltonian at thermal equilibrium:

$\mathcal{H}_{0}=\mathcal{H}_{Z}+\mathcal{H}_{Q}.$ (1)

Here, $\mathcal{H}_{Z}$ and $\mathcal{H}_{Q}$ stand for the Zeeman and
quadrupolar interaction terms, respectively. The former represents the direct
coupling between the nuclear intrinsic magnetic moment
$\gamma\hslash\mathbf{I}$ and the magnetic field $\mathbf{B}_{0}$ externally
applied in the laboratory:

$\mathcal{H}_{Z}=-\gamma\hslash\mathbf{I}\cdot\mathbf{B}_{0},$ (2)

where $\gamma$ is the gyromagnetic ratio of the spin, while
$\hslash\mathbf{I}$ is the spin operator of the nucleus. The term
$\mathcal{H}_{Q}$, on the other hand, represents the interaction between the
electric quadrupole moment of the nucleus and the electric field gradient
(EFG) generated by the surrounding electrons, and in the coordinate system of
the principal axes of the EFG it reads:

$\mathcal{H}_{Q}=\frac{e^{2}qQ}{4I(2I-1)}\left(3I_{Z}^{2}-I(I+1)+\frac{1}{2}\eta(I_{+}^{2}+I_{-}^{2})\right)$
(3)

where $e$ is the elementary charge, $eq=V_{ZZ}$ is the largest eigenvalue of
the EFG tensor, $eQ$ is the electric quadrupole moment, and $\eta$ is the
asymmetry parameter of the EFG.

In NMR/NQR, one probes these interactions by sending a pulse of radiation onto
the system, which accounts for a perturbing term to be included in the full
Hamiltonian:

$\mathcal{H}_{1}(t)=\left(2\mathbf{B}_{1}\cos(2\pi\nu_{P}t-\varphi_{P})\right)\cdot\mathbf{I}$
(4)

where $2\mathbf{B}_{1}$ is the magnetic component of the radiation pulse,
$\nu_{P}$ and $\varphi_{P}$ its frequency and phase, respectively. This
$B_{1}$ radiation field is in a plane perpendicular to externally applied
magnetic field, $B_{0}$, that defines Zeeman quantization axis.

Before the application of any pulse, the system is found in its thermal
equilibrium state, which at room temperature can be approximated as
$\rho(0)=\exp(-\mathcal{H}_{0}/k_{B}T)/\mathcal{Z}=(1-\mathcal{H}_{0}/k_{B}T)/\mathcal{Z}$.
Computing the evolution of this state under the action of a pulse is
equivalent to finding the corresponding evolution operator $U(t_{P},0)$, where
$t_{P}$ is the time duration of the pulse. In the most general case, this
operator cannot be computed directly, since the full Hamiltonian
$\mathcal{H}=\mathcal{H}_{0}+\mathcal{H}_{1}(t)$ may depend on time. The
procedure followed by our software consists of two steps:

  1. 1.

The problem is cast to the interaction picture, where the only relevant term
of the Hamiltonian is the perturbing one:

$\mathcal{H}(t)\rightarrow\mathcal{H}_{\mathcal{H}_{0}}(t)=\exp(i\mathcal{H}_{0}t/\hslash)\mathcal{H}_{1}(t)\exp(-i\mathcal{H}_{0}t/\hslash).$
(5)

  2. 2.

The evolution operator is approximated retaining the first few terms of the
Magnus expansion [25, 26], computed through the following formulas:

$\Omega_{1}(t_{P},0)=\int_{0}^{t_{P}}dt_{1}\widetilde{\mathcal{H}}(t_{1})$ (6)
$\Omega_{2}(t_{P},0)=\frac{1}{2}\int_{0}^{t_{P}}dt_{1}\int_{0}^{t_{1}}dt_{2}[\widetilde{\mathcal{H}}(t_{1}),\widetilde{\mathcal{H}}(t_{2})]$
(7)

from which the evolution operator is readily computed as

$U(t_{P},0)=\exp(\Omega_{1}(t_{P},0)+\Omega_{2}(t_{P},0)+\Omega_{3}(t_{P},0)).$
(8)

The software is also able to generate and analyze the typical observable
measured in real NMR/NQR experiments, i.e. the free induction decay signal
(FID). In laboratories, the FID is the electrical signal induced in a coil
wound around the sample after the electromagnetic RF pulse which generates
$B_{1}$ field is switch-off. This signal is simply related to the component of
the sample’s magnetization along the axis of the coil, related to the FID
itself [27]. If the coil is oriented along $\mathbf{\hat{n}}$, then the FID
signal will be given by:

$S(t)=\mathrm{Tr}\left[\rho(t)\mathbf{\hat{n}}\cdot\mathbf{I}\exp(-t/T_{2})\right]\qquad
t>t_{P}$ (9)

where we replaced the magnetization with the spin operator $\mathbf{I}$ of a
single nucleus in the ensemble, since they are equal up to a scaling factor.
Moreover, we introduced an empirical exponential decay $\exp(-t/T_{2})$ to
account for the decoherence of spin ensemble. In effect, our software
generates a complex FID whose imaginary part represents the signal produced by
an additional coil orthogonal to $\mathbf{\hat{n}}$, so that in the most usual
situation the FID reads
$S(t)=\mathrm{Tr}\left[\rho(t)I_{+}\exp(-t/T_{2})\right]$.

Once the FID is acquired, one typically computes its Fourier transform to
obtain what is called the NMR/NQR spectrum. This is the main outcome of the
experiment and provides information about the interactions experienced by the
system and the that of transitions occurred in its evolution. This is shown by
the expansion of the FID in Fourier components:

$S(t)=\displaystyle\sum_{\varepsilon,\eta}\matrixelement{\varepsilon}{I_{+}}{\eta}\matrixelement{\eta}{\rho(t_{P})}{\varepsilon}\exp(i\omega_{\varepsilon,\eta}t)$
(10)

where $\varepsilon$, $\eta$ run over the energy eigenvalues of the system,
$\ket{\varepsilon}$, $\ket{\eta}$ are the corresponding eigenstates, and
$\omega_{\varepsilon,\eta}$ is the frequency of transition between these two.
This formula proves that the peaks of the NMR/NQR spectrum are located at the
resonance frequencies $\omega_{\varepsilon,\eta}$, and that some of these
frequencies may not show up in the spectrum if the associated transition has
not occurred or the detection setup is not oriented appropriately.

##  3 Structure and usage of the software

PULSEE is not simply a simulator of the time evolution of nuclear spin states,
but also reproduces all the main features of NMR/NQR experiments. In this way,
the program is a valuable tool in experimental research, as the outcomes of a
simulation are generated in a form that can be directly compared with the
results measured in laboratories.

###  3.1 Download, dependencies and launching

The software can be downloaded from the following GitHub repository:
https://github.com/DavideCandoli/PULSEE

PULSEE has been written entirely in Python 3.7. The program makes wide use of
many of the standard Python modules (namely numpy, scipy, pandas, matplotlib)
for its general purposes. Tests have been carried out using the pytest
framework and the hypothesis module. The software includes a GUI which has
been implemented with the tools provided by the Python library kivy. In order
to run the GUI, the installation of the additional modules kivy.garden and
garden.matplotlib.backend_kivy is required.

The GUI is launched in the directory Code entering the following command from
terminal

      $ python PULSEE_GUI.py

Otherwise, one may use the functions defined in the module Simulation to write
a custom simulation, as outlined in subsection 3.3.

###  3.2 Modules of the software

The program consists of 6 modules. Below, the content and role of each module
is reported briefly:

  1. 1.

Operators

This module, together with Many_Body, is to be considered a sort of toolkit
for the simulation of generic quantum systems. It contains the definition of
Python classes and functions related to the basic mathematical objects which
enter the treatment of a quantum system. Operators concerns with the
simulation of a single system, while Many_Body extends it to systems made up
of several particles.

  2. 2.

Many_Body

It contains the definitions of the functions tensor_product and partial _trace
which allow to pass from a single particle Hilbert space to a many particle
space, and viceversa.

  3. 3.

Nuclear_Spin

Here is where the classes representing the spin of an atomic nucleus or a
system of nuclei are defined.

  4. 4.

Hamiltonians

This file includes the definition of the relevant terms of the Hamiltonian of
a nuclear spin system in a typical NMR/NQR experiment: namely, these are the
Zeeman interaction, the quadrupolar interaction, the $J$-coupling between
nuclei, and the interaction with an RF pulse of radiation.

  5. 5.

Simulation

This is the module the user should refer to in order to implement a custom
simulation. The functions defined here allow to set up the nuclear system,
evolve it under the action of a sequence of pulses, generate the FID signal,
and compute the NMR spectrum from it.

  6. 6.

NMR_NQR_GUI

This is the graphical user interface of the program. The GUI provides a simple
and intuitive way to perform a simulation, but it has limited features with
respect to a custom simulation code: in order to reproduce complicated
experiments involving long pulse sequences or many spin systems, one should
make use of the latter.

###  3.3 Building up a simulation

The starting point of any simulation is the set up of the system under study,
which is done by calling the function nuclear_system_setup:

    nuclear_system_setup(spin_par, quad_par=None, zeem_par=None, \
                         j_matrix=None, initial_state=’canonical’,\
                         temperature=1e-4)

This function returns three objects representing the spin system, the
unperturbed Hamiltonian, and the initial state, respectively.

The next step is to evolve the state of the system under the action of a pulse
of radiation, a task carried out by the function evolve:

    evolve(spin, h_unperturbed, dm_initial, \
           mode=None, pulse_time=0, \
           picture=’RRF’, RRF_par={’nu_RRF’: 0,
                                   ’theta_RRF’: 0,
                                   ’phi_RRF’: 0}):

The function not only allows to specify the features of the pulse applied to
the system, but also the reference frame where the evolution is computed.

Once the evolved state is obtained as the outcome of evolve, one can generate
the FID signal associated with this state by calling the function FID_signal

    FID_signal(spin, h_unperturbed, dm, acquisition_time, T2=100, \
               theta=0, phi=0, reference_frequency=0)

The arguments of this function allow the user to set the time window of
acquisition of the FID, the decoherence time $T_{2}$, the orientation, and the
frequency of rotation of the detecting coils.

Eventually, one computes the NMR/NQR spectrum from the FID signal by passing
the FID through the function fourier_transform_signal

    fourier_transform_signal(times, signal, frequency_start, \
                             frequency_stop, opposite_frequency=False)

The module Simulation includes also the functions for plotting the density
matrix of the evolved state as well as the FID signal and NMR/NQR spectrum.

##  4 Examples of execution

In this section we illustrate some noteworthy simulations performed with
PULSEE. In addition to being valid examples of the execution of the code,
these simulations have been chosen because they clearly demonstrate the
precision of NMR and NQR in the control of nuclear spin degrees of freedom,
which reflects their accuracy in the determination of unknown nuclear
interactions in a sample under study.

###  4.1 Selective transitions between quadrupolar states by means of properly
polarized pulses in NQR experiments

The structure of the energy spectrum of quadrupolar nuclei allows for the
selective excitation of its states by applying a pulse of radiation with the
proper polarization.

A first notable example is represented by the pure NQR of a spin 3/2 nuclei
whose energy levels and available transitions are depicted in Figure 1. This
system may undergo two single photon transitions at the same frequency, namely
$\ket{1/2}\leftrightarrow\ket{3/2}$ and $\ket{-1/2}\leftrightarrow\ket{-3/2}$,
in contrast with a pure NMR experiment, where all the transitions are
characterised by the same variation of the magnetic quantum number $\Delta m$.
These two transitions imply an opposite change in the angular momentum of the
system, so that each of them can occur only under exchange of a photon with
circular polarization (c.p.) $\sigma^{+}$ and $\sigma^{-}$ respectively.
Therefore, when one irradiates the system by a linearly polarized (l.p.)
resonant pulse, both transitions will be induced; in contrast, by choosing the
proper polarization of the pulse one is able to select only one of the two.
The potential of circularly and in general elliptically polarized RF pulses in
NQR has been widely explored [28, 29, 30].

Figure 1: Scheme of the energy spectrum and the available transitions of a
quadrupolar nucleus of spin 3/2. The transitions labeled with $\lvert\Delta
m\rvert=1$ involve the exchange of a single photon, while those labeled
$\lvert\Delta m\rvert=2$ involve two photons.

These theoretical expectations are correctly reproduced by our software. We
simulated the pure NQR of a spin 3/2 35Cl nuclei in a potassium chlorate
crystal (KClO3), whose gyromagnetic ratio is $\gamma/2\pi=4.17$ MHz/T and the
quadrupolar resonance frequency is $\nu_{Q}=28.1$ MHz [31]. We prepared the
system in the initial state depicted in Figure 2. Then we performed two
distinct simulations evolving the system under the action of a $\pi$ pulse
with polarization $\sigma^{+}$ or $\sigma^{-}$ respectively (in a classical
picture such pulses rotate the initial nuclear magnetization by 180∘ clockwise
and anticlockwise respectively), and the results obtained are shown in Figure
3. We note that the $\pi$ pulse is defined as that whose amplitude $B_{1}$ and
time duration $t_{P}$ satisfy the equation

$\gamma\alpha B_{1}t_{P}=\pi$ (11)

where $\gamma$ is the gyromagnetic ratio of the nucleus and
$\alpha=\sqrt{I(I+1)-m(m+1)}$ is a factor depending on the transition being
induced.

Figure 2: Real part of the density matrix representing the initial state of
the spin 3/2 nucleus in the simulation of a pure NQR experiment. The system is
prepared in a classical distribution where the states at the ground level are
equally populated. Figure 3: Real part of the simulated density matrix
representing the evolved state of the spin 3/2 quadrupolar nucleus after: (a)
a $\pi$ pulse with polarization $\sigma^{+}$; (b) a $\pi$ pulse with
polarization $\sigma^{-}$.

The evolved density matrices clearly show that a pulse with circular
polarization $\sigma^{+}$ ($\sigma^{-}$) couples only with the transition
between states $\ket{1/2}\leftrightarrow\ket{3/2}$
($\ket{-1/2}\leftrightarrow\ket{-3/2}$) by acting selectively on the two
relevant energy eigenstates to induce a full inversion of their respective
populations, in such a way that the total angular momentum is conserved.

Another experiment where the pulse can be set up to selectively induce
transitions is the NQR of a spin 1 nucleus in the presence of an asymmetric
EFG [5]. Due to the non-vanishing asymmetry parameter, the energy eigenstates
of this system are no longer the spin eigenstates, but they read:

$\ket{0}\qquad\&\qquad\ket{\xi_{\pm}}=\left(\ket{1}\pm\ket{-1}\right)/\sqrt{2}$
(12)

The energies of these states and the frequencies of transition between them
are displayed in Figure 4.

Figure 4: Scheme of the energy spectrum and the available transitions for a
quadrupolar nucleus of spin 1 with an asymmetric EFG. The subscripts of the
transition frequencies $\nu_{x/y/z}$ refer to the direction of linear
polarization of the pulse required to induce each transition.

What is peculiar with this system is that in order for each of the three
transitions to occur, the pulse should have a distinct linear polarization.
Indeed, developing calculations one finds that

$I_{x}=\blockarray{cccc}\ket{\xi_{+}}&\ket{0}\ket{\xi_{-}}\\\
\block{(ccc)c}010\ket{\xi_{+}}\\\ 100\ket{0}\\\ 000\ket{\xi_{-}}\\\ \quad
I_{y}=\left(\begin{array}[]{ccc}0&0&0\\\ 0&0&1\\\ 0&1&0\\\
\end{array}\right)\quad I_{z}=\left(\begin{array}[]{ccc}0&0&1\\\ 0&0&0\\\
1&0&0\\\ \end{array}\right)$ (13)

from which it is easy to prove that an $\mathbf{\hat{x}}$-,
$\mathbf{\hat{y}}$\- or $\mathbf{\hat{z}}$-polarized pulse will only affect
the transition $\ket{\xi_{+}}\leftrightarrow\ket{0}$,
$\ket{\xi_{-}}\leftrightarrow\ket{0}$ or
$\ket{\xi_{+}}\leftrightarrow\ket{\xi_{-}}$, respectively. This behaviour can
be assessed by observing the spectrum generated by each of these pulses:
recalling equation (10), once the proper orientation of the detection coils is
set, one is able to visualize if a certain transition has occurred, depending
on the whether the term $\matrixelement{\eta}{\rho(t_{P})}{\varepsilon}$
vanishes or not.

These results have been simulated in a fictitious spin 1 nucleus with
$e^{2}qQ/h=1\,\rm{MHz}$ and asymmetry $\eta=0.6$, for which the transition
frequencies are $\nu_{x}\equiv\nu_{\xi_{+}\leftrightarrow 0}=0.9\,\rm{MHz}$,
$\nu_{y}\equiv\nu_{\xi_{-}\leftrightarrow 0}=0.6\,\rm{MHz}$ and
$\nu_{z}\equiv\nu_{\xi_{+}\leftrightarrow\xi_{-}}=0.3\,\rm{MHz}$. The
resulting NQR spectra are displayed in Figure 5.

Figure 5: Spectra resulting from three distinct simulations of the NQR of a
spin 1 nucleus in an asymmetric EFG, where different pulses have been applied
with polarization along $\mathbf{\hat{x}}$ (a), $\mathbf{\hat{y}}$ (b) and
$\mathbf{\hat{z}}$ (c), respectively.

###  4.2 Generation of quantum coherences in a spin 3/2 quadrupolar nucleus

Let us consider the spin 3/2 35Cl nucleus in the same KClO3 crystal introduced
above. In the previous example, we showed how to induce a full inversion of
the populations of two of its energy eigenstates, say $\ket{m=1/2}$ and
$\ket{m+1=3/2}$, by means of a $\sigma^{+}$ c.p. $\pi$ pulse of radiation. In
general, when the angle on the right-hand side of Equ. (11) is set to a value
different from $n\pi$, where $n$ is an integer, the final density matrix
exhibits non-zero off-diagonal elements, meaning that the evolved state
includes a quantum superposition of $\ket{m}$ and $\ket{m+1}$. In NMR, such
elements are typically called “single quantum coherences”, where “single”
specifies the fact that $\Delta m=1$.

In Figure 6 we display the results of three simulated experiments where a
pulse resonant with the $\ket{1/2}\leftrightarrow\ket{3/2}$ transition is
applied and the angle in Equ. (11) is set to the values $\pi/3$, $\pi/2$ and
$2\pi/3$, respectively. These simulations demonstrate it is possible to fine-
tune the amplitudes of two states linked by a single-photon transition through
the careful manipulation of the parameters of the pulse.

Figure 6: Results of the simulation of the NQR of 35Cl nuclei in a KClO3
crystal where a pulse resonant with the $\ket{1/2}\leftrightarrow\ket{3/2}$
transition is applied. The three histograms show the real part of the density
matrix of the system evolved after a: (a) $\pi/3$ pulse; (b) $\pi/2$ pulse;
(c) $2\pi/3$ pulse. As the angle of rotation approaches $\pi$, the populations
of the states involved in the transition undergo a continuous exchange, and at
the same time non-zero off-diagonal elements emerge between them, meaning that
the two states are in a quantum superposition.

###  4.3 Preparation of an ensemble of spin 3/2 nuclei in a pseudopure state
by means of NQR

NMR and NQR are naturally suited to the implementation of simple quantum
information processors, as they are an efficient and high-precision method for
manipulation of nuclear spins [32, 33, 34, 35, 36]. Nonetheless, in typical
NMR/NQR experiments the system under study is a macroscopic sample made up of
a huge number of nuclei, which makes it impossible to prepare it in a pure
state as would be required by an ordinary quantum computation protocol. This
problem has been addressed following a different strategy [7]: by means of a
properly designed pulse sequence, the ensemble of nuclear spins can be
prepared in a pseudopure state, _i.e._ a state
$\rho=a\mathbbm{1}+b\outerproduct{\psi}{\psi}$ which differs from a pure state
by a term proportional to the identity. A state of this kind is called
pseudopure because, under evolution, it behaves like a full-fledged pure
state. This property makes it the ideal starting point for any NMR/NQR quantum
computation protocol [37, 38, 39].

In what follows, we describe a simulation of the NQR protocol aimed at
realizing a 2-qubit pseudopure state in the ensemble of spin 3/2 quadrupolar
35Cl nuclei of a KClO3 crystal [22], whose energy levels and available
transitions have been already shown in Figure 1. The states of the 2-qubit
computational basis correspond to the spin ones as follows:

$\begin{array}[]{ccc}\ket{00}&\equiv&\ket{3/2}\\\ \ket{01}&\equiv&\ket{1/2}\\\
\ket{10}&\equiv&\ket{-1/2}\\\ \ket{11}&\equiv&\ket{-3/2}.\\\ \end{array}$ (14)

Here, we remark that the simulated protocol is not aimed at implementing the
pseudopure state in the physical system itself: such a state is obtained as
the average of the results of three distinct experiments (as depicted in
Figure 7), following a common practice employed in NMR/NQR called temporal
averaging [40]. In each of the three experiments, the system is handled in a
distinct way:

  1. 1.

In the first, the system is left in its original thermal equilibrium state.

  2. 2.

In the second, the system is irradiated by a c.p. pulse with resonant
frequency $\nu_{Q}=(E_{\pm 3/2}-E_{\pm 1/2})/h$, inducing one of the single
photon transitions ($\lvert\Delta m=1\rvert$). The time duration of the pulse
is set to a value such that the populations of the states linked by the
transition are exchanged.

  3. 3.

In the third, a c.p. pulse at half the resonance frequency is applied,
yielding one of the two-photons transitions ($\lvert\Delta m=2\rvert$). Again,
the time duration of the pulse accounts for the exchange of the populations of
the states linked by the transition.

Figure 7: Population diagrams of the states to be combined through temporal
averaging in order to realise the $\ket{11}=\ket{-3/2}$ pseudopure state in a
spin 3/2 quadrupolar nucleus.

If the polarization of the pulses applied in the steps 2. and 3. are
appropriately chosen the average of the density matrices resulting from the
three experiments will have the properties of a pseudopure state belonging to
the computational basis in Equ. (14). The outcomes of the simulation,
illustrating the real part of the density matrices representing the four
pseudopure states of the computational basis of 2 qubits, are shown in Figure
8.

Figure 8: Outcomes of the protocol for the realization of pseudopure states in
an ensemble of spin 3/2 quadrupolar nuclei, as simulated by our software. The
histograms display the real part of the density matrices representing the four
pseudopure states of the computational basis of 2 qubits.

###  4.4 NQR and NMR implementation of a CNOT gate on a couple of qubits

Implementing a CNOT gate in the system we have already discussed in subsection
4.3 is a straightforward task. That is, in a 2-qubit system with $\ket{0}$ and
$\ket{1}$ as only allowed input values for both qubits, the CNOT gate flips
the second (target) qubit from $\ket{0}$ to $\ket{1}$ if and only if the first
(control) qubit is in $\ket{1}$ initial state. Indeed, one can easily check
that the action performed by a CNOT1 gate on the 2-qubit system is equivalent
to that of a pulse which yields an exchange of the populations of the states
$\ket{-1/2}$ and $\ket{-3/2}$, as illustrated in Figure 9. The effect of this
gate as simulated by our software is depicted in Figure 10.

Figure 9: Symbolic notation of a CNOT gate operating on the state $\ket{10}$
(on top) and the action of the pulse which carries out the equivalent
operation on the NQR version of the 2-qubit system (at the bottom). Figure 10:
Input (on the left) and output (on the right) states of the simulated NQR
CNOT1 gate when the initial state is $\ket{10}$.

It is possible to implement an analogous operation by means of NMR as well,
but in a different nuclear system. As explained in [35], this time the 2
qubits are encoded in 2 distinct spin 1/2 nuclei (following the convention
$\ket{0}\equiv\ket{1/2}$, $\ket{1}\equiv\ket{-1/2}$), and in order for them to
work as a control-target qubit couple they must interact with each other.
Thus, we assume that they are linked by the typical $J$-coupling, whose
contribution to the Hamiltonian is:

$\mathcal{H}_{J}=hJI_{z}^{(1)}I_{z}^{(2)}$ (15)

where $J$ is the coupling constant and $I_{z}^{(i)}$ is the $z$ component of
the spin of the $i$-th nucleus. The experimental protocol for the
implementation of an NMR CNOT gate employs both selective rotations of each
spin as well as the free evolution of the whole system under the action of
$J$-coupling, according to the sequence:

$\text{CNOT}_{1}=\left(-\frac{\pi}{2}\right)^{\mathbf{I}_{1}}_{z}\left(\frac{\pi}{2}\right)^{\mathbf{I}_{2}}_{z}\left(-\frac{\pi}{2}\right)^{\mathbf{I}_{2}}_{x}U\left(\frac{1}{2J}\right)\left(-\frac{\pi}{2}\right)^{\mathbf{I}_{2}}_{y}$
(16)

Here, factors of the type $\left(\alpha\right)^{\mathbf{I}_{i}}_{x/y/z}$
represent pulses resonant with the $i$-th spin which make it rotate an angle
$\alpha$ around the axis specified in the subscript; $U\left(1/2J\right)$, on
the other hand, stands for the free evolution of the system for a time
duration of $1/2J$. We point out that in order to be able to perform selective
rotations of one of the two spins, the nuclei’s gyromagnetic ratios must be
appreciably different, leading to well separated gyromagnetic frequencies
$\nu_{0}^{(i)}=-(\gamma^{(i)}/2\pi)B_{0}$.

We have carried out a simulation of this protocol starting from ideal pure
input states. The outcomes match closely our expectations, _i.e._ initial ket
is flipped since control (second) qubit was in $\ket{1}$ initial state, as is
shown in Figure 11.

Figure 11: Input (on the left) and output (on the right) of the simulated NMR
CNOT1 gate when the initial state is set to the ideal pure state $\ket{10}$.
The output state presents a slight discrepancy with respect to the expected
$\ket{11}$ state, which is thought to be a consequence of the discrete
approximations taken in the simulation.

##  5 Conclusions

PULSEE is an open-source software for the simulation of nuclear
magnetic/quadrupole resonance experiments (NMR/NQR) on solid-state samples.
The main purpose of this program is to provide a numerical tool for the
development of new methods of investigation of emergent properties in complex
materials inspired by the NMR/NQR protocols established in the context of
quantum information processing [41].

The software follows the principles of wide accessibility and intuitive
utilization, being available for download from a public GitHub repository [24]
and providing a GUI, as well as a complete and detailed documentation.

The examples of execution we have discussed were meant to illustrate the
features of the software, which include the ability to simulate both the
evolution of spin states and the corresponding experimental observables and
also highlight the possibilities to manipulate nuclear spin states through
NMR/NQR, which could be relevant for the design of highly sensitive protocols
for the study of emergent quantum properties of materials.

##  6 CRediT author statement

Davide Candoli: Investigation, Software Programing, Development and
Validation, Formal analysis, Visualization, Data curation, Writing-Original
draft preparation. Samuele Sanna: Methodology, Investigation, Supervision,
Formal analysis, Resources, Reviewing and Editing. Vesna F. Mitrović:
Conceptualization, Methodology, Supervision, Formal analysis, Visualization,
Resources, Writing-Reviewing and Editing.

##  7 Aknowledgments

We thank Prof. Enrico Giampieri and Dr. Stephen Carr for their helpful advice
during the development of the program. We are grateful to Jonathan Frassineti
and Ilija Nikolov for their feedback as the very first users of the program.
We also thank Prof. Paolo Santini and Prof. Alessandro Chiesa for reading the
manuscript and providing helpful feedback. V. F. M. acknowledges support form
the U.S. National Science Foundation grants OIA-1921199 and DMR-1905532.

## References

  * [1] L. M. K. Vandersypen, I. L. Chuang, NMR techniques for quantum control and computation, Rev. Mod. Phys. 76 (2005) 1037–1069.  doi:10.1103/RevModPhys.76.1037.
URL https://link.aps.org/doi/10.1103/RevModPhys.76.1037

  * [2] K. Modi, A. Brodutch, H. Cable, T. Paterek, V. Vedral, The classical-quantum boundary for correlations: Discord and related measures, Rev. Mod. Phys. 84 (2012) 1655–1707.  doi:10.1103/RevModPhys.84.1655.
URL https://link.aps.org/doi/10.1103/RevModPhys.84.1655

  * [3] E. N. Kaufmann, R. J. Vianden, The electric field gradient in noncubic metals, Rev. Mod. Phys. 51 (1979) 161–214.  doi:10.1103/RevModPhys.51.161.
URL https://link.aps.org/doi/10.1103/RevModPhys.51.161

  * [4] W. P. Halperin, Quantum size effects in metal particles, Rev. Mod. Phys. 58 (1986) 533–606.  doi:10.1103/RevModPhys.58.533.
URL https://link.aps.org/doi/10.1103/RevModPhys.58.533

  * [5] A. Abragam, Principles of Nuclear Magnetism, Oxford University Press, 1961.
  * [6] R. Blinc, Magnetic resonance and relaxation in structurally incommensurate systems, Phys. Rep. 79 (1981) 331.
  * [7] D. G. Cory, A. F. Fahmy, T. F. Havel, Ensemble quantum computing by NMR spectroscopy, Proceedings of the National Academy of Sciences 94 (5) (1997) 1634–1639.  arXiv:http://www.pnas.org/content/94/5/1634.full.pdf, doi:10.1073/pnas.94.5.1634.
URL http://www.pnas.org/content/94/5/1634

  * [8] I. L. Chuang, L. M. K. Vandersypen et al., Experimental realization of a quantum algorithm, Nature 393 (1998) 143–146.  doi:https://doi.org/10.1038/30181.
  * [9] J. A. Jones, M. Mosca, Implementation of a quantum algorithm on a nuclear magnetic resonance quantum computer, The Journal of Chemical Physics 109 (5) (1998) 1648–1653.  doi:10.1063/1.476739.
URL http://dx.doi.org/10.1063/1.476739

  * [10] J. A. Jones, M. Mosca, R. H. Hansen, Implementation of a quantum search algorithm on a quantum computer, Nature 393 (6683) (1998) 344–346.  doi:10.1038/30687.
URL http://dx.doi.org/10.1038/30687

  * [11] L. M. K. Vandersypen, M. Steffen, G. Breyta, C. S. Yannoni, R. Cleve, I. L. Chuang, Experimental realization of an order-finding algorithm with an nmr quantum computer, Physical Review Letters 85 (25) (2000) 5452–5455.  doi:10.1103/physrevlett.85.5452.
URL http://dx.doi.org/10.1103/PhysRevLett.85.5452

  * [12] G. Long, H. Yan, Y. Li, C. Tu, J. Tao, H. Chen, M. Liu, X. Zhang, J. Luo, L. Xiao, X. Zeng, Experimental NMR realization of a generalized quantum search algorithm, Physics Letters A 286 (2) (2001) 121–126.
  * [13] N. Sinha, T. S. Mahesh, K. V. Ramanathan, A. Kumar, Toward quantum information processing by nuclear magnetic resonance: Pseudopure states and logical operations using selective pulses on an oriented spin 3/2 nucleus, The Journal of Chemical Physics 114 (10) (2001) 4415–4420.  arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.1346645, doi:10.1063/1.1346645.
URL https://aip.scitation.org/doi/abs/10.1063/1.1346645

  * [14] T. F. Havel, D. G. Cory, S. Lloyd, N. Boulant, E. M. Fortunato, M. A. Pravia, G. Teklemariam, Y. S. Weinstein, A. Bhattacharyya, J. Hou, Quantum information processing by nuclear magnetic resonance spectroscopy, American Journal of Physics 70 (3) (2002) 345–362.  arXiv:https://doi.org/10.1119/1.1446857, doi:10.1119/1.1446857.
URL https://doi.org/10.1119/1.1446857

  * [15] T. Xin, B.-X. Wang, K.-R. Li, X.-Y. Kong, S.-J. Wei, T. Wang, D. Ruan, G.-L. Long, Nuclear magnetic resonance for quantum computing: Techniques and recent achievements, Chinese Physics B 27 (2) (2018) 020308.
URL http://stacks.iop.org/1674-1056/27/i=2/a=020308

  * [16] I. Oliveira, T. Bonagamba, R. Sarthour, J. Freitas, E. de Azevedo, NMR Quantum Information Processing, Elsevier, 2007.
  * [17] K. R. K. Rao, T. S. Mahesh, A. Kumar, Efficient simulation of unitary operators by combining two numerical algorithms: An NMR simulation of the mirror-inversion propagator of an $XY$ spin chain, Phys. Rev. A 90 (2014) 012306.  doi:10.1103/PhysRevA.90.012306.
URL https://link.aps.org/doi/10.1103/PhysRevA.90.012306

  * [18] J. Teles, R. Auccaise, C. Rivera-Ascona, A. G. Araujo-Ferreira, J. P. Andreeta, T. J. Bonagamba, Spin coherent states phenomena probed by quantum state tomography in Zeeman perturbed nuclear quadrupole resonance, Quantum Information Processing 17 (7) (2018) 177.  doi:10.1007/s11128-018-1947-1.
URL https://doi.org/10.1007/s11128-018-1947-1

  * [19] D. Lu, K. Li, J. Li, H. Katiyar, A. J. Park, G. Feng, T. Xin, H. Li, G. Long, A. Brodutch, J. Baugh, B. Zeng, R. Laflamme, Enhancing quantum control by bootstrapping a quantum processor of 12 qubits, npj Quantum Information 3 (1) (2017) 45.  doi:10.1038/s41534-017-0045-z.
URL https://doi.org/10.1038/s41534-017-0045-z

  * [20] A. Allouche, G. Pouzard, Computer simulation of ft-nmr multiple pulse experiment, Computer Physics Communications 54 (1) (1989) 171–176.  doi:https://doi.org/10.1016/0010-4655(89)90042-8.
URL https://www.sciencedirect.com/science/article/pii/0010465589900428

  * [21] M. Bak, J. T. Rasmussen, N. C. Nielsen, Simpson: A general simulation program for solid-state nmr spectroscopy, Journal of Magnetic Resonance 147 (2) (2000) 296–330.  doi:https://doi.org/10.1006/jmre.2000.2179.
URL https://www.sciencedirect.com/science/article/pii/S1090780700921797

  * [22] D. Possa, A. C. Gaudio, J. C. C. Freitas, Numerical simulation of NQR/NMR: Applications in quantum computing, Journal of Magnetic Resonance 209 (2) (2011) 250–260.
  * [23] F. A. Perras, C. M. Widdifield, D. L. Bryce, QUEST – Quadrupolar Exact Software: A fast graphical program for the exact simulation of NMR and NQR spectra for quadrupolar nuclei, Solid State Nuclear Magnetic Resonance 45-46 (2012) 36–44.  doi:https://doi.org/10.1016/j.ssnmr.2012.05.002.
URL https://www.sciencedirect.com/science/article/pii/S0926204012000586

  * [24] D. Candoli, PULSEE (Program for the simULation of nuclear Spin Ensemble Evolution (2021).
URL https://github.com/DavideCandoli/PULSEE

  * [25] S. Blanes, F. Casas, J. A. Oteo, J. Ros, A pedagogical approach to the Magnus expansion, Eur. J. Phys. 31 (2010) 907–918.
  * [26] C. P. Slichter, Principles of Magnetic Resonance, Springer-Verlag, 1990.
  * [27] P. J. Hore, J. A. Jones, S. Wimperis, NMR: The Toolkit. How Pulse Sequences Work, Oxford, 2015.
  * [28] M. J. Weber, E. L. Hahn, Selective spin excitation and relaxation in nuclear quadrupole resonance, Phys. Rev. 120 (1960) 365–375.  doi:10.1103/PhysRev.120.365.
URL https://link.aps.org/doi/10.1103/PhysRev.120.365

  * [29] Y. Lee, H. Robert, D. Lathrop, Circular Polarization Excitation and Detection in 14N NQR, Journal of Magnetic Resonance 148 (2) (2001) 355–362.
  * [30] J. Miller, B. Suits, A. Garroway, Circularly polarized rf magnetic fields for spin-1 nqr, Journal of Magnetic Resonance 151 (2) (2001) 228–234.  doi:https://doi.org/10.1006/jmre.2001.2366.
URL https://www.sciencedirect.com/science/article/pii/S1090780701923663

  * [31] T. P. Das, E. L. Hahn, Nuclear Quadrupole Resonance Spectroscopy.
  * [32] D. Cory, R. Laflamme, E. Knill, L. Viola, T. Havel, N. Boulant, G. Boutis, E. Fortunato, S. Lloyd, R. Martinez, C. Negrevergne, M. Pravia, Y. Sharf, G. Teklemariam, Y. Weinstein, W. Zurek, Nmr based quantum information processing: Achievements and prospects, Fortschritte der Physik 48 (9‐11) (2000) 875–907.  doi:https://doi.org/10.1002/1521-3978(200009)48:9/11<875::AID-PROP875>3.0.CO;2-V.
  * [33] K. V. R. M. Murali, N. Sinha, T. S. Mahesh, M. H. Levitt, K. V. Ramanathan, A. Kumar, Quantum-information processing by nuclear magnetic resonance: Experimental implementation of half-adder and subtractor operations using an oriented spin-7/2 system, Phys. Rev. A 66 (2002) 022313.  doi:10.1103/PhysRevA.66.022313.
URL https://link.aps.org/doi/10.1103/PhysRevA.66.022313

  * [34] H. Kampermann, W. S. Veeman, Characterization of quantum algorithms by quantum process tomography using quadrupolar spins in solid-state nuclear magnetic resonance, The Journal of Chemical Physics 122 (21) (2005) 214108.  arXiv:https://doi.org/10.1063/1.1904595, doi:10.1063/1.1904595.
URL https://doi.org/10.1063/1.1904595

  * [35] I. S. Oliveira, T. J. Bonagamba, et al., NMR Quantum Information Processing, Elsevier Science, 2007.
  * [36] R. M. Serra, I. S. Oliveira, Nuclear magnetic resonance quantum information processing, Phil. Trans. R. Soc. A 370 (2012) 4615–4619.  doi:https://doi.org/10.1098/rsta.2012.0332.
  * [37] K. Dorai, Arvind, A. Kumar, Implementing quantum-logic operations, pseudopure states, and the deutsch-jozsa algorithm using noncommuting selective pulses in nmr, Phys. Rev. A 61 (2000) 042306.  doi:10.1103/PhysRevA.61.042306.
URL https://link.aps.org/doi/10.1103/PhysRevA.61.042306

  * [38] N. Sinha, T. S. Mahesh, K. V. Ramanathan, A. Kumar, Toward quantum information processing by nuclear magnetic resonance: Pseudopure states and logical operations using selective pulses on an oriented spin 3/2 nucleus, The Journal of Chemical Physics 114 (10) (2001) 4415–4420.  arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.1346645, doi:10.1063/1.1346645.
URL https://aip.scitation.org/doi/abs/10.1063/1.1346645

  * [39] Y.-P. Tan, X.-F. Nie, J. Li, H.-W. Chen, X.-Y. Zhou, X.-H. Peng, J.-F. Du, Preparing pseudo-pure states in a quadrupolar spin system using optimal control, Chinese Physics Letters 29 (12) (2012) 127601.  doi:10.1088/0256-307x/29/12/127601.
URL https://doi.org/10.1088/0256-307x/29/12/127601

  * [40] E. Knill, I. Chuang, R. Laflamme, Effective pure states for bulk quantum computation, Phys. Rev. A 57 (1998) 3348–3363.  doi:10.1103/PhysRevA.57.3348.
URL https://link.aps.org/doi/10.1103/PhysRevA.57.3348

  * [41] I. Nikolov, A. D. Maestro, C. Ramanathan, V. F. Mitrović, Spin Squeezing as a Probe of Quantum Materials, APS - March Meeting Bulletin - V44.00014.