# Deep learning for improved keV-scale recoil identification in high
resolution gas time projection chambers

J. Schueler,11footnotetext: Corresponding author. M. Ghrear  S. E. Vahsen  P.
Sadowski  C. Deaconu

###### Abstract

Recoil-imaging gaseous time projection chambers (TPCs) with directional
sensitivity are attractive for dark matter (DM) searches. Detectors capable of
reconstructing 3D nuclear recoil directions would be uniquely sensitive to the
predicted dipole angular distribution of DM recoils. Observation of this
directional distribution would unambiguously establish the galactic origin of
a claimed DM signal. Recoil directionality also would provide powerful
discrimination against background recoils from solar neutrino scattering.
These advantages of directionality can only be exploited however, if electron
recoil backgrounds from gamma rays can be sufficiently suppressed. We
introduce a deep learning-based recoil event classifier that uses a 3D
convolutional neural network (3DCNN) to identify event species based on their
recoil images. We compare electron background rejection performance of
discriminants determined by the 3DCNN both to the traditional discriminant of
track length, as well as discriminants obtained from state-of-the-art shallow
learning methods. We train the 3DCNN classifier using recoil charge
distributions with ionization energies ranging from
0.5-$10.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$, for
$25\text{\,}\mathrm{c}\mathrm{m}$ of drift in an 80:10:10 mixture of $\rm
He$:$\rm CF_{4}$:$\rm CHF_{3}$. The charges are initially segmented into
$(100\times 100\times 100)$$\text{\,}\um^{3}$ bins when determining track
length and the shallow learning discriminants, but are rebinned with a reduced
segmentation of about $(850\times 850\times 850)$$\text{\,}\um^{3}$ for the
3DCNN. Despite the courser binning, compared to using track length, we find
that classifying events with the 3DCNN reduces electron backgrounds by a
factor of up to 1,000 and effectively reduces the energy threshold of our
simulated TPC by $30\%$ for fluorine recoils and $50\%$ for helium recoils. We
also find that the 3DCNN reduces electron backgrounds by up to a factor of 20
compared to the shallow machine learning approaches, corresponding to a
$2\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ reduction in the energy
threshold.

##  1 Introduction

Directly detecting the constituents of dark matter (DM) remains one of the key
goals of contemporary physics. Directional detection would provide unique,
robust, and unambiguous confirmation of the galactic origin of a signal in the
form of a dipole distribution in galactic coordinates [1, 2].

As more regions of DM parameter space are being ruled out, and cross-section
limits approach the solar neutrino fog, directional detection is also
increasingly of interest for maximizing DM sensitivity within the fog, and for
studying the solar neutrinos themselves. Detectors capable of reconstructing
the full 3D vector direction, energy, and time of individual recoil events are
preferable in this context. Such “recoil-imaging” detectors are widely
applicable and hence garnering increasing interest [3, 4].

Gaseous time projection chambers (TPCs) are the most mature recoil imaging
technology with event time measurement. Unlike liquid noble gas detector,
which benefit from strong self-shielding, one key issue in gas detectors is
the rejection of electron recoil backgrounds from gamma rays. TPCs with high-
definition charge readout (HD TPCs) are, however, capable of reconstructing
the topology of both nuclear recoils and electron recoils in great detail,
which we expect should maximize particle ID capabilities and again help to
reject such backgrounds. A key question then, is whether the electron
rejection in gas TPCs is sufficient to achieve a background free DM search for
a given experiment and exposure.

Studies of the proposed $1000\text{\,}\mathrm{m}^{3}$ gas TPC experiment
CYGNUS [5], running for six years, suggest the electron background will be
flat in energy, and electron backgrounds must be rejected offline by factors
exceeding $6\times 10^{4}$ per $\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$.
Due to the diffusion in gas-based detectors, electron identification
performance falls exponentially with recoil energy. As a result, the effective
analysis energy threshold above which a large gas TPC will remain background
free is determined by the particular energy at which the electron rejection
meets the performance requirement. A relatively small improvement in electron
rejection, and consequent reduction of the analysis energy threshold can have
large benefits on DM reach, because the expected DM recoil spectrum is
generally steeply falling with energy.

Because recoil-imaging HD TPCs offer extremely rich three-dimensional charge
density measurements, one obvious question is how to best exploit this
information to maximize recoil identification. This is a problem ideally
suited for machine learning techniques.

There have been many efforts to improve electron background rejection for keV-
scale recoil events, which often use multivariate combinations of discriminant
observables for event classification. The MIMAC group, for instance, used
boosted decision trees (BDTs) [6] to improve electron rejection factors by a
factor of 20 over more traditional methods [7]. More recently, we introduced a
set of nine observables that are based on the shape of the recoil charge cloud
measured in a 3D recoil-imaging TPC [8]. A joint combination of these nine
event-shape variables using a hard-cuts-based approach led to up to a two
order of magnitude improvement in electron rejection over using the length
along a track’s principal axis, which is a common observable for event
classification.

Given the electron background rejection improvement observed both by
multivariate combinations of discriminant variables and the usage of event-
shape variables, we attempt to combine the best aspects of both of these
approaches by introducing a deep learning-based classifier for event
identification. In particular, we introduce a 3D convolutional neural network
(3DCNN, Section 3) that is directly fed the 3D ionization density distribution
of events, binned into a $32\times 32\times 32$ voxel grid, and outputs class
probabilities of the recoil species of the event. We expect this end-to-end
approach to enable us to capture more information than using predefined
observables [9, 10, 11, 12], thus leading to better background rejection
performance. In Section 4, we separately combine the nine event-shape
observables with two shallow learning techniques, (1) a BDT and (2) a fully
connected neural network (FCNN), both of which also output recoil-species
class probabilities that can be used as multivariate classification
discriminants. Comparing the electron background rejection performance between
these shallow learning-produced discriminants and the hard-cuts-based combined
observable from Ref. [8] allows us to assess the relative effectiveness of
different techniques of combining the nine event-shape observables. Comparing
the electron rejection performance of the 3DCNN with the best-performing
combination of these observables provides insight on if any crucial event-
shape information is missing from the set of nine observables.

##  2 Overview of simulation

Before introducing our new classifiers, we first describe the simulated
detector and recoil characteristics used in our study. We build off of
previous work from our group and generate a large simulation sample using
identical parameters to those in Ref. [8]. Doing this allows us to directly
compare the electron rejection performance of our new classifiers with already
established electron rejection improvements.

We use SRIM [13] and retrim [14] to simulate recoiling He and F nuclei and
DEGRAD [15] to simulate electron recoils in an 80:10:10 mixture of
He:$\text{CF}_{4}$:$\text{CHF}_{3}$ at a total pressure of
$60\text{\,}\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{r}$ and temperature of
$25^{\circ}$C. The SRIM computation requires a compound correction for every
nucleus in the gas mixture (He, H, C, and F). The calculation of these
compound corrections and their associated values are detailed in Ref. [8]. The
retrim step requires the average energy per ion pair, $W$, and Fano factor,
$\mathcal{F}$, which were calculated using Garfield++/Heed [16] as
$W=$35\text{\,}\mathrm{e}\mathrm{V}$$ and $\mathcal{F}=0.19$, respectively.
DEGRAD can simulate recoils isotropically; however, for SRIM and retrim we
isotropize the recoil simulations after generating the initial ionization
distributions. We assume a drift field of
$40.6\text{\,}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ parallel to the drift
direction and use Magboltz [17] to determine the corresponding transverse and
longitudinal diffusion coefficients as
$(\sigma_{T},\sigma_{L})=(398,425)$$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}\mathrm{/}\sqrt{cm}$.
The diffusion is applied to our simulations assuming a uniform drift length of
$25\text{\,}\mathrm{c}\mathrm{m}$.

After applying diffusion, we assume each individual electron is detected. We
do not simulated any charge amplification or digitization, but instead
simulate a highly efficient, pixelated readout by binning each primary track
into $($100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$)$ voxels. We note that
our choice of a $40.6\text{\,}\mathrm{V}\mathrm{/}\mathrm{c}\mathrm{m}$ drift
field corresponds to $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$ per clock
cycle on a $40\text{\,}\mathrm{M}\mathrm{H}\mathrm{z}$ clock which is
equivalent to the readout refresh rate used in existing TPCs with pixel
readout [18].

##  3 Convolutional neural network classifier

Recoil-imaging TPCs with high readout segmentation reconstruct detailed 3D
images of ionization charge distributions of recoil tracks. These tracks have
many identifiable characteristics that can be used for particle
identification. For sufficiently high energy recoil tracks, the length of the
track along its principal axis is often enough to reliably distinguish between
electron recoils and nuclear recoils, as electron recoils tend to create
longer tracks following a more meandering path than nuclear recoils of
equivalent energy. At lower energies, however, diffusion during drift has a
proportionately larger effect on the overall event topology, leading to more
spherical ionization distributions. Figure 1 shows
$\sim$$6\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ examples of each of the
three recoil species investigated here. Comparing the background electron
recoil (left) with the signal fluorine and helium recoils (middle and right,
respectively), we see that both of the nuclear recoil species have a dense
cluster of charge near the center of the event, whereas the electron recoil
appears to be a more diffuse charge cloud with larger gaps. It’s difficult to
unambiguously identify the principal axis direction of low energy tracks,
leading to poor angular resolution. Furthermore, a few stray charges can
significantly bias the length along the principal axis, which may lead to
event misclassification with this traditional approach. Since there are so
many features present in these tracks, using a deep learning classifier where
the classifier learns the best patterns for event selection on its own,
independent of an identified principal axis, is attractive. To this end, we
construct and train a 3DCNN classifier that is directly fed 3D voxel grids of
the ionization distribution of individual recoil events with the recoil
species assigned as the class label of the event.

Figure 1:  3D charge distributions of
$6\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ electron, fluorine, and helium
recoils after diffusion. The top row shows the original
$($100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro
m}$\times$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$)$ binning of these
events, which is fine enough that there is only a single charge in each filled
bin. The bottom row shows the events re-binned into the $(32\times 32\times
32)$ voxel grid that are input into the 3DCNN.

###  3.1 Building blocks of feature extraction in 3DCNNs

Figure 2:  Schematic representation of the dimensionality reduction of a 3D
convolution operation with a convolutional stride of $S=(2,2,2)$. Here
$N\in\mathbb{R}^{16\times 16\times 16}$, and $K\in\mathbb{R}^{2\times 2\times
2}$. The black lines connecting the two red shaded regions illustrate the
action of $K$ on $N[14,0,14]$, a $2\times 2\times 2$ sub-block of $N$, leading
to the $1\times 1\times 1$ element $C_{14,0,14}$ of $C\in\mathbb{R}^{8\times
8\times 8}$. Each index $\alpha$, $\beta$, and $\gamma$ of $C$ runs along
$\\{0,2,4,6,8,10,12,14\\}$ due to the convolutional stride of $S=(2,2,2)$. The
pink regions show the next element of the convolution operation after striding
2 units along the $y$ direction on $N$.

3DCNNs are not yet widely used for natural image classification due to the
large amounts of GPU memory required to store and perform operations on large
3D images with high voxel density. The lower energy recoil events that are of
interest in directional DM searches, however, have their entire event
topologies contained within a relatively small region, enabling the use of
3DCNNs with current-generation hardware, even in detectors with high spatial
resolution. Indeed, since convolutional neural networks are translationally
invariant, given a set of 3D recoil tracks, we can shift the coordinates of
each recoil track to a fixed origin without any change in performance of a
3DCNN up to uncertainties in bin assignment. Translating events to a common
origin allows for the fiducial binning of track charges into a much smaller 3D
grid of voxels than binning the entire fiducial volume of a detector. Here we
briefly introduce 3D convolutions and 3D pooling; two core building blocks of
our 3DCNN.

####  3.1.1 3D Convolutions

In 3DCNNs, the 3D convolution operation is a 3D cross correlation of an input
grid of data, $N$, with a kernel (also called a filter), $K$, that slides
along the input grid. To be more specific, let $N\in\mathbb{R}^{n_{x}\times
n_{y}\times n_{z}}$ (often called an input feature map), and let
$K\in\mathbb{R}^{k_{x}\times k_{y}\times k_{z}}$, where $k_{j}\leq n_{j};\quad
j=x,y,z$. We denote the shapes of $N$ and $K$ as $(n_{x}\times n_{y}\times
n_{z})$ and $(k_{x}\times k_{y}\times k_{z})$, respectively. The action of $K$
on $N[\alpha,\beta,\gamma]$, a $k_{x}\times k_{y}\times k_{z}$ sub-block of
$N$, is an element, $C_{\alpha\beta\gamma},$ of the output “convolved” feature
map, $C$, and is a mapping from $\mathbb{R}^{k_{x}\times k_{y}\times
k_{z}}\rightarrow\mathbb{R}$ defined by

$\displaystyle(K\star
N[\alpha,\beta,\gamma])_{\alpha\beta\gamma}=\sum_{l=0}^{k_{x}-1}\sum_{m=0}^{k_{y}-1}\sum_{n=0}^{k_{z}-1}N_{\alpha+l,\beta+m,\gamma+n}K_{lmn}\equiv
C_{\alpha\beta\gamma},$ (3.1)

where $\alpha+k_{x}\leq n_{x}$, $\beta+k_{y}\leq n_{y}$, $\gamma+k_{z}\leq
n_{z}$, and $(K\star N[\alpha,\beta,\gamma])_{\alpha\beta\gamma}$ is a single
element of the 3D cross correlation operation.

$C$ is formed by sliding $K$ in integer steps of $S_{x}$, $S_{y}$ and $S_{z}$
along the $x$, $y$, and $z$ extents of $N$, respectively. In particular, $C$
is an ordered grid composed of all
$\lfloor\frac{n_{x}-k_{x}}{S_{x}}+1\rfloor\times\lfloor\frac{n_{y}-k_{y}}{S_{y}}+1\rfloor\times\lfloor\frac{n_{z}-k_{z}}{S_{z}}+1\rfloor$
elements of $(R\star N[r_{x},r_{y},r_{z}])_{r_{x}r_{y}r_{z}}$, where
$r_{j}\in\\{0,S_{j},2S_{j},\ldots,\lfloor\frac{n_{j}-k_{j}}{S_{j}}+1\rfloor\\};\quad
j=x,y,z$, and $\lfloor\rfloor$ represents the floor operator, telling us that
we only compute elements of $C$ when the shape of $K$ is entirely contained
within $N$. We call $S_{x}$, $S_{y}$, and $S_{z}$, the convolutional stride in
$x$, $y$, and $z$, and will more compactly use $S\equiv(S_{x},S_{y},S_{z})$ to
denote the convolutional stride of a 3D convolution operation. Figure 2 shows
an example of the shapes of input and output feature maps, where a $2\times
2\times 2$ filter acts on a $16\times 16\times 16$ input feature map with a
convolutional stride of $S=(2,2,2)$, leading to an $8\times 8\times 8$ output
grid. We can also use zero padding where we pad the outer perimeter of our
input grid, $N$ with $P=(P_{x},P_{y},P_{z})$ layers of zeros. Here, $P_{j}$ is
a positive integer or zero and $j=x,y,z$. More generally, a $k_{x}\times
k_{y}\times k_{z}$ kernel acting with a stride of $(S_{x},S_{y},S_{z})$ on an
$n_{x}\times n_{y}\times n_{z}$ input feature map with additional zero padding
of $(P_{x},P_{y},P_{z})$, will lead to an output feature map of shape
$(C_{x}\times C_{y}\times C_{z})$, where

$\displaystyle
C_{j}=\left\lfloor\frac{n_{j}-k_{j}+2P_{j}}{S_{j}}+1\right\rfloor;\quad
j=x,y,z.$ (3.2)

In a 3DCNN, the elements of the convolutional filters, $K_{lmn}$ are learnable
parameters, so the network attempts to learn the features of interest to
extract from a feature map.

####  3.1.2 3D Pooling

Pooling is an operation commonly used in convolutional neural networks to
downsample a feature map while still retaining important information for
classification. In our 3DCNN we include an average pooling layer (called
AvgPool; see Section 3.3). Similar to a 3D convolution, the average pooling
operation involves a 3D filter of a given size sliding over an input feature-
map, $N$, but the average pooling operation simply computes the mean of all
elements contained in the sub-block of $N$ that the filter is sliding over.
The average pooling operation will thus lead to an output feature map $C$ of
shape given by Equation (3.2). Since average pooling simply computes the mean
of all elements in a sub-block of $N$, there are no learnable parameters
associated with average pooling.

###  3.2 Data processing for the 3DCNN

We simulate 12,380,422 electron recoil events ranging in ionization energies
between $0.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ and
$10.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$, 367,984 F recoil events,
and 338,909 He recoils. The nuclear recoils are simulated with a uniform
energy spectrum between $5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{r}$ and
$50\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{r}$. While our simulated detector
is sensitive to individual electrons, nuclear recoils still lose energy due to
quenching, so in the absence of reconstructing the true recoil energy, we
compute the ionization energy of all nuclear recoil species in terms of
electron equivalent energy to allow for a direct comparison of electron
rejection at a given ionization energy in our detector. Expressed in terms of
electron equivalent energy, the ionization energy of a nuclear recoil
$E_{\text{ionization}}$, is computed as

$\displaystyle E_{\text{ionization}}=N_{e}W_{e},$ (3.3)

where $N_{e}$ is the number of electrons produced in the recoil event, and
$W_{e}$ is the average energy per electron-ion pair. Following the approach in
Ref. [8] we set $W_{e}$ to $32.4\text{\,}\mathrm{e}\mathrm{V}$. After
computing $E_{\text{ionization}}$, we restrict our nuclear recoil sample to
$$0.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$\leq
E_{\text{ionization}}\leq$10.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$$,
leaving us with 224,299 F recoils and 93,774 He recoils.

| Recoil
---
species
| Original energy
---
range [$\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{r}$]
| Energy
---
range [$\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$]
|  $\\#$ Events
---
training set
|  $\\#$ Events
---
validation set
|  $\\#$ Events
---
test set
e | —- | 0.5-10.5 | 2,947,702 | 491,429 | 8,941,291
F | 5-50 | 0.5-10.5 | 53,760 | 8,812 | 161,727
He | 5-50 | 0.5-10.5 | 22,443 | 3,743 | 67,588
Table 1:  Event samples used with the 3DCNN classifier. Events in each sample
are binned in $1\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ wide bins
centered at integer ionization energy steps between
$1\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ and
$10\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$. We train and evaluate the
3DCNN classifier separately for each energy bin.

The charges in our simulated events are already binned with $(100\times
100\times 100)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$
segmentation, however because diffusion is larger than our bin size, we use a
low density gas, and we do not simulate any charge amplification, we find that
the charge clouds are relatively sparse at this resolution. While it is
feasible to train a 3DCNN on current-generation hardware using the native
$(100\times 100\times
100)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$ resolution of our
simulation, we opt to reduce the spatial segmentation of our events to a
$32\times 32\times 32$ grid of bins evenly spaced within a cube of width
$2.72\text{\,}\mathrm{c}\mathrm{m}$, leading to bin sizes of about $(850\times
850\times 850)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$. We do
this for two reasons: First, all fluorine recoil events are entirely contained
within this cube, so no recoil information is lost due to cropping, and
second, the reduced spatial segmentation leads to more than a factor of 600
reduction in grid volume compared to the original $(100\times 100\times
100)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$ segmentation,
leading to a substantial reduction in computational cost. Comparing the top
and bottom rows of Figure 1 it appears that the structure of recoil tracks
isn’t significantly altered with this reduced resolution, suggesting that
event classification performance for the simulated detector configuration will
not be significantly hampered by our reduction of spatial segmentation. We
note that uncertainties in bin placement from re-binning may slightly alter
performance compared to binning the unbinned initial coordinates of 3D charge
into the $32\times 32\times 32$ grid.

We use the PyTorch [19] software library for all neural network computations
and store each event as a tuple containing two entries: (1) a $(32\times
32\times 32)$ voxel grid filled with binned charge that is stored as a PyTorch
tensor data structure, and (2) an integer representing the class label of the
event. We label electron recoils as 0, F recoils as 1, and He recoils as 2.
Charge in each bin of the voxel grid is stored as unsigned 8 bit integers
leading to an effective dynamic range of 0 to 255 electrons per voxel. No bins
are saturated in any event with this dynamic range.

###  3.3 Network architecture

Figure 3:  Architecture of our 3DCNN classifier. Within the feature extraction
portion, each cube represents a unique feature map (represented by the
different shade of each map), with the size of the cubes shown approximately
to scale in reference to the dimensions of the input image. We show the bin
content of one such-feature map in each layer to represent the steps of a
convolutional chain. Solid black lines connect larger pale red cubes, which
illustrate the convolutional or pooling filters acting on a portion of the
input feature map, to the smaller magenta cubes which are the outputs of the
convolutional filter. The final 16 feature maps are then flattened and passed
into a fully connected neural network for classification. The FCNN diagram was
produced using [20].

Figure 3 and Table 2 together outline the network architecture of our 3DCNN.
As is shown in Figure 3, our network consists of a series of convolutional and
pooling layers for feature extraction, followed by a dense fully connected
neural network (FCNN) for event classification. Within the feature extraction
portion of the neural network, we employ five convolutional blocks (ConvBlocks
1-5) and a pooling layer (AvgPool) to downsample the feature maps while still
maintaining important features. Each of the five convolutional blocks contain
the following components: (i) A 3D convolution with a convolutional filter
size listed in Table 2, (ii) a 3D batch normalization [21], (iii) a scaled
exponential linear unit (SELU) activation layer [22], and (iv) a randomized
dropout [23] of 0.03 to reduce overfitting.

Layer |  | $\\#$
---
Filters
Filter size | $S$ |  |  $\\#$ Learnable
---
parameters
Dropout | Output shape
ConvBlock 1 | 4 | $(2\times 2\times 2)$ | (2,2,2) | 44 | 0.03 | $(4\times 16\times 16\times 16)$
ConvBlock 2 | 8 | $(3\times 3\times 3)$ | (1,1,1) | 888 | 0.03 | $(8\times 14\times 14\times 14)$
ConvBlock 3 | 16 | $(3\times 3\times 3)$ | (1,1,1) | 3,504 | 0.03 | $(16\times 12\times 12\times 12)$
ConvBlock 4 | 32 | $(3\times 3\times 3)$ | (1,1,1) | 13,920 | 0.03 | $(32\times 10\times 10\times 10)$
AvgPool | 32 | $(2\times 2\times 2)$ | (2,2,2) | —- | —- | $(32\times 5\times 5\times 5)$
ConvBlock 5 | 16 | $(3\times 3\times 3)$ | (1,1,1) | 13,872 | 0.03 | $(16\times 3\times 3\times 3)$
FC1 | —- | —- | —- | 27,712 | 0.05 | $(64\times 1)$
FC2 | —- | —- | —- | 1,040 | 0.05 | $(16\times 1)$
FC3 | —- | —- | —- | 51 | —- | $(3\times 1)$
Table 2:  More specific details of each layer shown in Figure 3. We assume a
single ($1\times 32\times 32\times 32$) image is fed into the network. The
output shape column gives the shape of the output after each layer. The output
of ConvBlocks 1-5 and AvgPool is a tensor of shape
$(\text{D}\times\text{L}\times\text{W}\times\text{H})$, where D is the layer
depth (number of feature maps), and L, W, and H are the length, width and
height of each feature map, respectively. $S$ is the convolutional stride of
the layer. The output shapes of the FC layers are 1 dimensional vectors for
each node in a given layer. In addition to the weights and biases associated
with each node in the FCNN, the entries within each convolutional filter are
also learnable parameters, so we list the total number of learnable parameters
associated with each layer in the network.

Walking through the network architecture, we start with a ($32\times 32\times
32$) input image. The vast majority of voxels in all input images are filled
with 0, so we implement a stride of 2 in ConvBlock 1 to immediately downsample
this block’s four output feature maps to ($16\times 16\times 16$) to reduce
the computational overhead of training and evaluating our network. Following
the chain in Figure 3, we next perform three successive convolutional blocks
(ConvBlocks 2-4) where we gradually increase the number of convolutional
filters employed in each layer to produce more feature maps for
classification. ConvBlocks 2-4 each use ($3\times 3\times 3$) convolutional
filters with a stride of 1 and no padding, so each of these convolutions
decreases the feature map dimension by two, leaving us with 32 feature maps of
size ($10\times 10\times 10)$ at the end of ConvBlock4. After this, we perform
average pooling. We use a ($2\times 2\times 2$) pooling filter with a stride
length of 2, leaving us with 32 ($5\times 5\times 5$) feature maps that we
feed into our final convolutional block (ConvBlock5) which leaves us with 16
($3\times 3\times 3$) feature maps. Since the values composing each
convolutional filter are learnable parameters, we expect that when the network
is trained, there will be useful features encoded in some of these feature
maps. We finally flatten these feature maps into a $((16\times 3\times 3\times
3)\times 1)=(432\times 1)$ vector, that contains each extracted feature at the
end of the convolutional chain. This flattened $(432\times 1)$ feature vector
is then fed into a fully connected dense neural network with two hidden layers
and output class assignments corresponding to e, F, and He recoils. Each of
the three fully connected layers (FC1, FC2, and FC3) use a SELU activation
function and FC1 and FC2 include a random dropout of 0.05. The raw model
output of each event is a $(3\times 1)$ vector, $\mathbf{z}$, with entries
corresponding to each of the three class outputs (e recoil, F recoil, and He
recoil). The softmax function

$\displaystyle\sigma(\mathbf{z})_{i}=\frac{e^{z_{i}}}{\sum_{j=1}^{3}e^{z_{j}}}.$
(3.4)

is applied to $\mathbf{z}$ to map the class outputs $z_{i}\in\mathbf{z}$ to
class probabilities. We henceforth label $\sigma(\mathbf{z})_{1}$,
$\sigma(\mathbf{z})_{2}$, and $\sigma(\mathbf{z})_{3}$, as $\rm p_{e}$, $\rm
p_{F}$, and $\rm p_{He}$, which represent the model-predicted class
probabilities of e, F, and He recoils, respectively.

###  3.4 Training the network

We first shuffle the order of all events and then split the data into distinct
training, validation, and test sets with 3,023,905 events in our training
sample, 503,984 events in our validation sample, and 9,172,513 events in our
test sample. We set the test sample aside and implement the following
procedure to train our model:

  1. 1.

Form a PyTorch tensor of shape $(256\times 32\times 32\times 32)$, which is a
minibatch consisting of 256 randomly selected voxel images from the training
sample.

  2. 2.

Feed the minibatch and the corresponding truth label of each image of the
minibatch into the 3DCNN.

  3. 3.

Use PyTorch’s built in CrossEntropy loss function to compute the loss of the
batch. We wish to minimize this loss. We use an Adam [24] optimizer with a
learning rate of 0.0002.

  4. 4.

Update model weights using backpropagation [25].

  5. 5.

Repeat steps 1-4 until we’ve run through all events in the training set. This
is called a training epoch.

  6. 6.

At the end of each training epoch, repeat steps 1-3 for the validation sample.
We do not implement step 4 as we don’t wish to train the 3DCNN on the
validation set. Compute the sum of the losses of each minibatch of the
validation sample.

  7. 7.

If the summed losses over the validation set minibatches are less than in the
previous epoch, we treat this as the model learning and save all model
weights.

  8. 8.

Implement early stopping [26] where steps 1-7 are repeated until the total
validation loss doesn’t decrease at all over 10 successive epochs.

  9. 9.

The model state corresponding to the epoch with the lowest validation loss is
our trained model.

We were able to train the 3DCNN with an Nvidia GeForce RTX 2070 consumer-grade
laptop GPU in less than 6 hours.

Figure 4:  Comparison of traditional (track length) and 3DCNN-based ($\rm
p_{e}$, $\rm p_{F}$ and $\rm p_{He}$) discriminants for events in the test
sample satisfying
$$6.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$<E_{\text{ionization}}<$7.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$$.
The shaded regions in the panels indicate events that are rejected after
selections corresponding to $50\%$ F recoil efficiency (dashed black vertical
line) and/or $50\%$ He recoil efficiency (solid black vertical line).

###  3.5 Interpreting model output

We use the 3DCNN output variables $\rm p_{e}$, $\rm p_{F}$, and $\rm p_{He}$
as event classification discriminants. Figure 4 shows a comparison of the
distributions of track length, $\rm p_{e}$, $\rm p_{F}$, and $\rm p_{He}$ for
all true F recoils (blue), all true He recoils (green), and all true electrons
(orange) in the test sample of data. Since $\rm p_{e}$, $\rm p_{F}$, and $\rm
p_{He}$ represent class probabilities, for each event, $\rm p_{e}+\rm
p_{F}+\rm p_{He}=1$, meaning if our 3DCNN is a good classifier, we should then
expect that each class probability peaks strongly toward 1 for its
corresponding class, and strongly toward 0 for the other two classes. This
expected behavior is observed in Figure 4, suggesting that $\rm p_{e}$, $\rm
p_{F}$, and $\rm p_{He}$ are strong event classification discriminants.

In addition to rejecting electron backgrounds, these discriminants can also be
used to distinguish nuclear recoil signal candidates. To quantify
classification performance, we define the signal efficiency as

$\displaystyle\varepsilon_{s}=\frac{N_{s}^{\prime}}{N_{s}},$ (3.5)

where $N_{s}$ is the total number of signal events in the sample and
$N_{s}^{\prime}$ is the number of remaining signal events after making a
selection. Now, suppose we want to quantify the signal purity of F recoils at
an F recoil efficiency of $\varepsilon_{F}=0.5$. We can do this by determining
$\rm p_{F}^{\prime}$, the median $\rm p_{F}$ of true F recoils. Then, we can
determine the number of F, He, and e recoils where $\rm p_{F}>p_{F}^{\prime}$,
and compute the F recoil purity as
$N_{\text{F}}^{\prime}/(N_{\text{e}}^{\prime}+N_{\text{F}}^{\prime}+N_{\text{He}}^{\prime})$.
The bottom left panel of Figure 4 shows an example of these selections for the
set of recoil events between $6.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$
and $7.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$, where the black
vertical dashed line corresponds to $\rm p_{F}^{\prime}$. We can perform an
analogous procedure with $\rm p_{He}$ to compute the He recoil at a fixed He
recoil efficiency.

In certain instances it may be advantageous to treat any nuclear recoil
species as signal and electrons as background. In these cases we could treat
$\rm p_{e}$ as a binary classification variable between electrons versus not
electrons (i.e. signal nuclear recoils for this ternary classification model).
For example, if we wanted to calculate the nuclear recoil purity at a nuclear
recoil efficiency of $\varepsilon_{\mathcal{R}}=0.5$, we could determine, $\rm
p_{e}^{\prime}$, the median value of $\rm p_{e}$ for all He or F recoils. We
would then compute
($N_{\text{F}}^{\prime}+N_{\text{He}}^{\prime})/(N_{\text{e}}^{\prime}+N_{\text{F}}^{\prime}+N_{\text{He}}^{\prime})$
where $N_{\text{F}}^{\prime}$, $N_{\text{He}}^{\prime}$, and
$N_{\text{e}}^{\prime}$ represent the number of events of each of these three
classes after applying a selection of $\rm p_{e}<\rm p_{e}^{\prime}$.

Since signal purity depends on the relative composition of recoil
constituents, we opt to mostly use the notion of a rejection factor, $R$, to
quantify background rejection performance. We define the sample size-
independent electron rejection factor $R_{\text{e}}$ as the ratio of the total
number of electrons in a sample $N_{\text{e}}$ to the number of electrons
remaining after a selection, $N_{\text{e}}^{\prime}$

$\displaystyle R_{\text{e}}=\frac{N_{\text{e}}}{N^{\prime}_{\text{e}}}.$ (3.6)

For the $7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ sample shown in Figure
4, when using length as the event classification variable we find rejection
factors of $135\pm 2$ and $870\pm 30$ for $\varepsilon_{\text{He}}=0.5$ and
$\varepsilon_{\text{F}}=0.5$, respectively. Using 3DCNN output probability as
our event classification variable, we find that only five electrons remain
both for He recoils at $\varepsilon_{\text{He}}=0.5$ and F recoils at
$\varepsilon_{\text{F}}=0.5$, leading to an electron rejection factor of
$141,400\pm 63,000$ for both of these cases, a more than 3 order of magnitude
improvement over the traditional track length discriminant for He recoils.

##  4 Shallow learning classifiers

Here we introduce the nine event-shape discriminants used in Ref. [8] and
train both a BDT (Section 4.2) and an FCNN (Section 4.3) to combine these
observables into multivariate classification discriminants.

###  4.1 Defining electron rejection discriminants

Each of the following variables are computed using our simulated TPC’s native
$(100\times 100\times
100)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$ resolution, as
opposed to the reduced $(850\times 850\times
850)$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$$ segmentation used
for the 3DCNN classifier.

  1. 1.

Length along the principal axis (LAPA): We use a singular value decomposition
(SVD) to identify the principal axis of the 3D track, and then take the
difference between the maximum and minimum of the track’s ionization
distribution coordinates projected onto this principal axis.

  2. 2.

Standard deviation of charge distribution (SDCD): The standard deviation of
the 3D position vectors of all charges in the event.

  3. 3.

Maximum charge density of the event ($\rho_{\text{Max}}$): This discriminant
is optimized by varying the bin width of the charge in cubic voxels and
recording the maximum amount of charge in a voxel. $\rho_{\text{Max}}$ was
separately optimized to maximize electron rejection for weakly-directional
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ F recoils and for directional
$12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ He recoils as described in
Ref. [8].

  4. 4.

Charge uniformity (QUnif): The standard deviation of the distribution of mean
distances between each charge and all other charges in a recoil event.

  5. 5.

Cylindrical thickness (CylThick): Sum of the squared transverse distances from
each charge to the principal axis of the track.

  6. 6.

Number of clusters (NClust): Number of clusters determined by the DBSCAN [27]
clustering algorithm. Two of the input parameters to this algorithm were
adjusted to find the pair that optimizes electron rejection both for weakly-
directional $7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ F recoils and for
directional $12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ He recoils as
described in Ref. [8].

  7. 7.

Clustering threshold (CThres): The threshold that the fraction of total event
charge in the largest cluster of an event must be above. Clustering is
performed DBSCAN and the threshold value is optimized separately to maximize
electron rejection for weakly-directional
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ F recoils and for directional
$12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ He recoils.

Given the separate directional and weakly-directional optimizations for the
$\rho_{\text{Max}}$ and ClustThres observables, we have a total of nine
discriminants to work with.

###  4.2 Boosted decision tree

We use the XGBClassifier model in the XGBoost [28] software library as our BDT
classifier for recoil classification. For this BDT analysis, we first
randomize the order of our data and then partition the data into a $23.8\%$ /
$4.0\%$ / $72.2\%$ training/validation/testing sample split. The validation
data sample is not used for the BDT analysis; however we still opt to create a
validation set so that we have an identical testing set to what’s later used
in the FCNN analysis. For each event, we compute the nine discriminants and
store the results as a $(9\times 1)$ vector to feed into the BDT classifier.
We assign class labels of 0 for e recoils, 1 for F recoils, and 2 for He
recoils.

$\rho_{\text{Max,F}}$ | SDCD | CThres${}_{\text{F}}$ | $\rho_{\text{Max,He}}$ | CThres${}_{\text{He}}$ | CylThick | QUnif | LAPA | NClust
---|---|---|---|---|---|---|---|---
0.437 | 0.231 | 0.169 | 0.111 | 0.032 | 0.009 | 0.007 | 0.003 | 0.002
Table 3:  Normalized feature importance of the nine discriminants in the
trained BDT.

Given the relatively quick training and evaluation time of XGBoost compared to
our 3DCNN classifier, we perform a course, two-step partial grid search to
optimize some of the hyperparameters in the XGBClassifier model. We make the a
priori decision to choose the set of hyperparameters
$(N_{\text{trees}},D_{\text{tree}},\ell)$ that maximizes the electron
rejection factor $R_{\text{e}}$ for
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ F recoils at
$\varepsilon_{\text{F}}=0.5$, where $N_{\text{trees}}$, $D_{\text{tree}}$, and
$\ell$ represent the number of trees, the maximum depth of a tree, and the
learning rate of the model, respectively. The first step of our grid search is
to fix $\ell$ at 0.2, and then train and evaluate our BDT using all 16
combinations of $N_{\text{trees}}\in\\{50,100,200,500\\}$ and
$D_{\text{tree}}\in\\{2,3,5,7\\}$. We find that $N_{\text{trees}}=200$ with
$D_{\text{tree}}=3$ to be the best combination at $\ell=0.2$, leading to
$R_{\text{e}}=101,000\pm 38,000$. Since this combination of $N_{\text{trees}}$
and $D_{\text{tree}}$ does not sit at the boundary of our grid, we keep this
combination and move on to our next step of optimizing $\ell$. In our second
step of the grid search, we train our classifier using
$\ell\in\\{0.1,0.15,0.2,0.25,0.3\\}$ and find our initial choice of $\ell=0.2$
to be the best, so we freeze our optimal set of hyperparameters for our BDT
classifier at $(N_{\text{trees}}=200,D_{\text{tree}}=3,\ell=0.2)$.

We trained our BDT model using 53,959 F recoils, 22,265 He recoils and
2,992,534 e recoils in our training sample222Data for the shallow learning
classifiers uses the same simulation as the 3DCNN, but was processed
independently, leading to small differences in the training sample statistics
compared to Table 1. We ensured the number of each recoil species in the test
sample here matches exactly with Table 1. Furthermore, when binning in integer
steps of electron equivalent energy, the number of each recoil species in the
test sample here is identical to the number of each recoil species in the
3DCNN test sample within each energy bin.. Similar to the 3DCNN classifier,
the BDT outputs ternary class probabilities, $\rm p_{e,BDT}$, $\rm p_{F,BDT}$,
and $\rm p_{He,BDT}$, allowing for the same model output interpretation as the
3DCNN (Section 3.5).

A benefit of using a BDT classifier is the ability to track the decisions made
throughout all branches in the trained classifier model. Table 3 shows the
normalized feature importance of each of the nine discriminants. Here we
define the feature importance for a particular observable to be the sum of the
information gain [29] of all splits where the observable is used, normalized
so that the sum of the feature importance of all nine observables is 1. The
two observables with the highest feature importance, $\rho_{\text{Max,F}}$ and
SDCD, are the same two discriminants that perform best in terms of electron
rejection factor versus $E_{\text{ionization}}$ at $50\%$ F recoil efficiency
in Ref. [8].

###  4.3 Fully connected neural network

We next train and evaluate the performance of an FCNN using identical
training, validation, and test sets to those generated for the BDT classifier
to derive a fair comparison between the two classifiers. Our FCNN consists of
an input layer with 9 nodes representing each of the electron rejection
observables we’ve defined, 3 hidden layers, and an output layer with 3 nodes
that represent our three output classes (see Table 4). We use a ReLU
activation with batch normalization after the input and hidden layers, and
apply a randomized dropout of 0.05 in all layers except the output layer. We
train this FCNN classifier by forming minibatches of 512 events, where each
event consists of a ($9\times 1$) PyTorch tensor of normalized $z$-scores of
each of the nine observables and the associated class label. We then follow
the same training procedure (steps 3-9 shown in Section 3) with the same early
stopping trigger that we used for the 3DCNN classifier training, except we use
a learning rate of 0.001.

Layer | Input shape | Output shape | Dropout
---|---|---|---
Input | ($9\times 1$) | ($32\times 1$) | 0.05
Hidden 1 | ($32\times 1$) | ($64\times 1$) | 0.05
Hidden 2 | ($64\times 1$) | ($128\times 1$) | 0.05
Hidden 3 | ($128\times 1$) | ($32\times 1$) | 0.05
Output | ($32\times 1$) | ($3\times 1$) | —-
Table 4:  FCNN classifier architecture. As mentioned in the text, the input
vectors are formed with normalized $z$-scores of each of the nine electron
rejection observables. We use a ReLU activation function and a randomized
dropout of 0.05 at the input and each of the hidden layers. When test sample
events are passed through the trained network, a softmax function is applied
to the model output to give class probabilities.

We performed a moderate amount of hand-tuning of the hyperparameters for the
FCNN, rather than a full hyperparameter search, and find that in general the
FCNN gives better electron rejection performance than the BDT classifier,
especially for He recoils. Similar to the 3DCNN, the output of each event
passed through the FCNN classifier is a $(3\times 1)$ vector containing
outputs associated with each of the three recoil species, so we apply a
softmax to this vector to obtain class probabilities $\rm p_{e,FCNN}$, $\rm
p_{F,FCNN}$, and $\rm p_{He,FCNN}$ that may be interpreted analogously to the
output class probabilities given by the 3DCNN and BDT.

##  5 Event identification performance results

We now have a set of class probabilities from the 3DCNN, BDT, and FCNN
classifiers that can be used to quantify and compare event selection
performance. We also include the multivariate observable from Ref. [8] in our
comparisons where available.

For a future directional DM detector, it is useful to consider electron
background rejection both as a function of ionization energy at fixed nuclear
recoil selection efficiency and as a function of nuclear recoil efficiency at
fixed energies. Comparing electron rejection performance versus energy
provides insight toward the energy threshold required for near background-free
operation at a given exposure, while comparing versus nuclear recoil
efficiency allows us to compare detector exposure required for a fixed signal
yield. Figures 5 and 6 show the results for both of these approaches.

Figure 5:  Comparison of electron-rejection performance of all discriminants
versus electron-equivalent energy at $50\%$ fluorine recoil efficiency (left)
and $50\%$ helium recoil efficiency (right). The top and bottom row plots
share a common horizontal axis. The black dashed lines show the combined
observable introduced in Ref. [8]. Data within a given energy bin in the
$R_{\text{e}}$ plots are staggered slightly to separate overlapping points.
The bottom bar plots show the number of remaining electrons after a given
recoil selection in each integer energy bin. The transparent bars indicate the
total number of electron recoils in each energy bin.

Figure 5 shows $R_{\text{e}}$ performance versus energy at a fixed $50\%$
nuclear recoil efficiency for F recoils (left) and He recoils (right). In
cases where all electrons are rejected, we plot $R_{\text{e}}$ as $1\times
10^{6}$ with zero uncertainty. The black dashed line labeled “No ML” is the
combined observable from Ref. [8]. For the three machine learning
discriminants and the traditional track length discriminant, we apply the
following procedure to determine $R_{\text{e}}$ in each energy bin:

  1. 1.

Determine the values, $\rm p_{\mathcal{R}}^{\prime}$, $\rm
p_{\mathcal{R},BDT}^{\prime}$, $\rm p_{\mathcal{R},FCNN}^{\prime}$, and track
length, $L_{\mathcal{R}}^{\prime}$, that correspond to a 50$\%$ efficiency of
recoil $\mathcal{R}$. In the left panel of Figure 5, $\mathcal{R}=\text{F}$,
and in the right panel, $\mathcal{R}=\text{He}$.

  2. 2.

$N_{\text{e}}$ is the total number of electrons within the energy bin and
$N_{\text{e}}^{\prime}$ for the 3DCNN, BDT, FCNN, and track length
classifiers, is the number of electrons in the energy bin satisfying
selections of $\rm p_{\mathcal{R}}>p_{\mathcal{R}}^{\prime}$, $\rm
p_{\mathcal{R},BDT}>p_{\mathcal{R},BDT}^{\prime}$, $\rm
p_{\mathcal{R},FCNN}>p_{\mathcal{R},FCNN}^{\prime}$, and
$L_{\mathcal{R}}<L_{\mathcal{R}}^{\prime}$, respectively.

  3. 3.

$R_{\text{e}}$ is computed for the 3DCNN, BDT, FCNN, and track length
classifiers with Eq. (3.6) using the $N_{\text{e}}^{\prime}$ associated with
the chosen classifier.

In general, we find that the 3DCNN outperforms all other methods of electron
background rejection, especially for He recoils. The CYGNUS collaboration
mentions in their feasibility study that electron rejection will effectively
determine their energy threshold [5], and argue that with a 6 year exposure
and flat electron background energy spectrum of
$1\text{\times}{10}^{4}\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}^{-1}\mathrm{y}\mathrm{e}\mathrm{a}\mathrm{r}^{-1}$,
they will be essentially background-free at energies corresponding to an
electron rejection factor in excess of $6\times 10^{4}$. With
$\varepsilon_{\text{F}}=0.5$, we find $R_{\text{e}}>6\times 10^{4}$ starting
around $10\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ using length as a
discriminant, compared to somewhere between
$6\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ and
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ with the 3DCNN, meaning the
use of the 3DCNN for electron rejection lowers the “background-free” energy
threshold by more than $3\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ over
using length. At $\varepsilon_{\text{He}}=0.5$ Ref. [8] found that
$R_{\text{e}}$ doesn’t exceed $6\times 10^{4}$ until around
$14\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ using length as a
discriminant, while with the 3DCNN classifier it exceeds $6\times 10^{4}$ at
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$. We note here that given the
logarithmic scale shown in Figure 5, electron rejection performance drops off
exponentially with energy, so a factor of two reduction of the “background-
free” energy threshold for He recoils when using the 3DCNN classifier over
using length is a significant improvement. The otherwise state-of-the-art FCNN
and BDT combinations of observables have “background-free” thresholds of
around $9\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ for He recoils so the
3DCNN reduces this threshold by around
$2\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$. Even a seemingly small
reduction in energy threshold has a significant impact on DM reach, given the
that expected DM recoil energy spectrum falls steeply with energy.

Figure 6:  Comparison of electron-rejection performance of all discriminants
versus nuclear recoil efficiency for
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ fluorine recoils (left) and
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ helium recoils (right). The
top and bottom row plots share a common horizontal axis. The black dashed line
in the upper left hand plot shows the combined observable from Ref. [8]. Data
points at a given recoil efficiency in the $R_{\text{e}}$ plots are staggered
slightly to separate overlapping points. The bottom bar plots show the number
of remaining electrons after a given recoil selection in each efficiency bin.

Figure 6 shows $R_{\text{e}}$ and $N_{\text{e}}^{\prime}$ after selections at
various F recoil efficiencies (left) and He recoil efficiencies (right) for
these same discriminants. We choose to evaluate $R_{\text{e}}$ and
$N_{\text{e}}^{\prime}$ at a fixed electron-equivalent energy range of
$$6.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$<E<$7.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$$
because this is the lowest energy range exceeding the CYGNUS background-free
criteria of $R_{\text{e}}>6\times 10^{4}$ at both $\varepsilon_{\text{F}}=0.5$
and $\varepsilon_{\text{He}}=0.5$. We use the same procedure to compute
$R_{\text{e}}$ and $N_{\text{e}}^{\prime}$ as in Figure 5, except now the
values of $\rm p_{\mathcal{R}}^{\prime}$, $\rm p_{\mathcal{R},BDT}^{\prime}$,
$\rm p_{\mathcal{R},FCNN}^{\prime}$, and $L^{\prime}$ are separately computed
for each $\varepsilon_{\mathcal{R}}$. Since the $\rho_{\text{Max,F}}$,
$\text{CThres}_{\text{F}}$, and NClust variables were optimized specifically
for $7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ F recoils, we expect this
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ bin will produce the best
relative $R_{\text{e}}$ of the shallow learning classifiers compared to the
3DCNN. We note here that the combined observable from Ref. [8] was computed at
$12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ for He recoils versus
$\varepsilon_{\text{He}}$. Given that the 3DCNN rejects all but two electrons
down to $8\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ at
$\varepsilon_{\text{He}}=0.5$ (Figure 5), we expect the 3DCNN to reject
essentially all $12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ electron
recoils, so we do not investigate the
$12\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ scenarios discussed in [8].

Figure 7:  Comparison of He and F recoil rejection at fixed F and He signal
efficiencies. $R_{\text{He}}$ and $R_{\text{F}}$ are rejection factors for He
and F recoils, respectively and are both defined analogously to
$R_{\text{e}}$. The bins in the upper right plot with no data points plotted
have all F recoils rejected.

Comparing each multivariate classifier constructed from combining the nine
event-shape variables, we find that introducing machine learning for
multivariate analysis leads to an improvement in electron rejection above
$\varepsilon_{\text{F}}\approx 0.2$. We find further improvement with the use
of the 3DCNN, where it improves $R_{\text{e}}$ by more than 2 orders of
magnitude over selections on length at higher F recoil efficiencies. For He
recoils, we find that the 3DCNN classifier significantly outperforms all other
models, often leading to an order of magnitude increase in $R_{\text{e}}$
compared to the BDT and FCNN combinations of the nine electron rejection
observables. Furthermore, we find that selecting on 3DCNN output probability
at $85\%$ He recoil efficiency leads to a higher $R_{\text{e}}$ than selecting
on length at $10\%$ He recoil efficiency, meaning that usage of the 3DCNN
classifier could allow an experiment to run with a factor of 8.5 smaller
exposure and obtain the same He recoil signal, at an improved background level
over using track length to classify events.

An additional benefit of training the machine learning models as ternary
classifiers is the ability to use the three output class probabilities to not
only reject electron backgrounds, but also classify nuclear recoil species.
Figure 7 summarizes the nuclear recoil classification ability of the three
machine learning classifiers. Above
$6\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ all three classifiers maintain
50$\%$ F recoil efficiency while rejecting greater than $99\%$ of He recoils.
Above $5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$, each classifier starts
to do a better job selecting pure samples of He recoils than F recoils,
especially the 3DCNN classifier. Selecting pure samples of He recoils is
desirable, as low energy He recoils tend to have better angular resolution
than F recoils of the same energy, so the ability to select a background-free,
almost entirely pure sample of He recoils down to
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ at
$\varepsilon_{\text{He}}=0.5$ is significant. To quantify this, there are
3,721 He recoils, 1 F recoil, and 5 e recoils in our sample that satisfy $\rm
p_{He}>p_{He}^{\prime}$, where $\rm p_{He}^{\prime}$ is the 3DCNN He class
probability corresponding to $\varepsilon_{\text{He}}=0.5$ in the
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ energy bin. This means we have
a He recoil purity in excess of $99.8\%$ for
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ He recoils at
$\varepsilon_{\text{He}}=0.5$. Though the fraction of each event type present
in our simulation doesn’t model a physically expected event composition, these
numbers suggest that use of the 3DCNN will lead to a high He recoil purity for
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ He recoils at
$\varepsilon_{\text{He}}=0.5$ in a directional DM experiment.

##  6 Conclusion

We have introduced a deep learning approach to event selection in 3D recoil-
imaging TPCs and compared its performance to other state-of-the-art machine
learning-based multivariate classifiers. Training a 3D convolutional neural
network on the 3D charge distributions of recoil tracks, we found a
significant improvement in electron rejection performance between
$0.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ and
$10.5\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ compared to shallow
learning methods that form multivariate classification variables from
predefined event-shape observables. Notably, the improved performance of the
3DCNN holds when using input data with a reduced spatial segmentation of
$(850\times 850\times
850)$$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$ compared to
$(100\times 100\times
100)$$\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}^{3}$ segmentation used to
compute the observables fed into the shallow learning classifiers. We find
that the 3DCNN classifier outperforms other methods by a larger margin for He
recoils than for F recoils indicating that the nine event-shape observables
used in the BDT and FCNN classifiers are better suited for electron rejection
in F recoil samples than in He recoil samples. Since the 3DCNN classifier
decides which features for event classification are most useful on its own,
using a classifier like this is advantageous, as it seems to have found more
useful sets of features to classify nuclear recoils (especially He recoils)
than our predefined observables. Defining the background-free energy threshold
of our simulated detector to be the lowest energy corresponding to
$R_{\text{e}}\geq 6\times 10^{4}$, we find our threshold to be between
$6\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ and
$7\text{\,}\mathrm{k}\mathrm{e}\mathrm{V}_{ee}$ for both F and He recoils at
$50\%$ nuclear recoil efficiency when using the 3DCNN output probability to
select for events. Using this 3DCNN to classify events effectively reduces the
energy threshold of our simulated detector by over $30\%$ for F recoils and
around $50\%$ for He recoils compared to classification using track length; a
significant improvement over traditional keV-scale recoil identification
techniques considering the steeply falling energy spectrum expected for WIMP
recoils. The 3DCNN classifier is also able to assign classification
probabilities to multiple recoil species simultaneously, making it a flexible
approach for event classification that should be of general interest for
directional DM searches.

## Acknowledgments

This work was supported by the U.S. Department of Energy (DOE) via Award
Number DE-SC0010504.

## References

  * Spergel [1988] David N. Spergel.  Motion of the earth and the detection of weakly interacting massive particles.  _Phys. Rev. D_ , 37:1353–1355, Mar 1988.  doi: 10.1103/PhysRevD.37.1353.  URL https://link.aps.org/doi/10.1103/PhysRevD.37.1353.
  * Mayet et al. [2016] F. Mayet et al.  A review of the discovery reach of directional Dark Matter detection.  _Phys. Rept._ , 627:1–49, 2016.  doi: 10.1016/j.physrep.2016.02.007.
  * Vahsen et al. [2021] Sven E. Vahsen, Ciaran A. J. O’Hare, and Dinesh Loomba.  Directional Recoil Detection.  _Ann. Rev. Nucl. Part. Sci._ , 71:189–224, 2021.  doi: 10.1146/annurev-nucl-020821-035016.
  * O’Hare et al. [2022] C. A. J. O’Hare et al.  Recoil imaging for dark matter, neutrinos, and physics beyond the Standard Model.  In _2022 Snowmass Summer Study_ , 3 2022.
  * Vahsen et al. [2020] S. E. Vahsen et al.  CYGNUS: Feasibility of a nuclear recoil observatory with directional sensitivity to dark matter and neutrinos.  8 2020.
  * Breiman et al. [1984] Friedman Breiman, L., Olshen J.H., R.A., and C.J. Stone.  _Classification And Regression Trees_.  Routledge, 1st edition, 1984.  doi: https://doi.org/10.1201/9781315139470.
  * Billard et al. [2012] J Billard, F Mayet, and D Santos.  Low energy electron/recoil discrimination for directional dark matter detection.  _Journal of Cosmology and Astroparticle Physics_ , 2012(07):020–020, Jul 2012.  ISSN 1475-7516.  doi: 10.1088/1475-7516/2012/07/020.  URL http://dx.doi.org/10.1088/1475-7516/2012/07/020.
  * Ghrear et al. [2021] Majd Ghrear, Sven E. Vahsen, and Cosmin Deaconu.  Observables for recoil identification in high-definition Gas Time Projection Chambers.  _JCAP_ , 10:005, 2021.  doi: 10.1088/1475-7516/2021/10/005.
  * Baldi et al. [2014] Pierre Baldi, Peter Sadowski, and Daniel Whiteson.  Searching for exotic particles in high-energy physics with deep learning.  _Nature Communications_ , 5, 2014.
  * Baldi et al. [2016] Pierre Baldi, Kevin Bauer, Clara Eng, Peter Sadowski, and Daniel Whiteson.  Jet substructure classification in high-energy physics with deep neural networks.  _Phys. Rev. D_ , 93:094034, May 2016.  doi: 10.1103/PhysRevD.93.094034.  URL https://link.aps.org/doi/10.1103/PhysRevD.93.094034.
  * Sadowski et al. [2017] P Sadowski, B Radics, Ananya, Y Yamazaki, and P Baldi.  Efficient antihydrogen detection in antimatter physics by deep learning.  _Journal of Physics Communications_ , 1(2):025001, 2017.  URL http://stacks.iop.org/2399-6528/1/i=2/a=025001.
  * Sadowski and Baldi [2018] Peter Sadowski and Pierre Baldi.  Deep learning in the natural sciences: applications to physics.  In _Braverman Readings in Machine Learning. Key Ideas from Inception to Current State_ , pages 269–297. Springer, 2018.
  * Ziegler et al. [2010] James F. Ziegler, M. D. Ziegler, and J. P. Biersack.  SRIM - The stopping and range of ions in matter (2010).  _Nuclear Instruments and Methods in Physics Research B_ , 268(11-12):1818–1823, June 2010.  doi: 10.1016/j.nimb.2010.02.091.
  * Deaconu [2015] Cosmin Ştefan Deaconu.  _A model of the directional sensitivity of low-pressure CF 4 dark matter detectors_.  PhD thesis, Massachusetts Institute of Technology, United States, January 2015.
  * [15] Stephen Biagi.  Degrad - transport of electrons in gas mixtures.  URL https://degrad.web.cern.ch/degrad/.
  * Veenhof and Schindler [2018] R. Veenhof and H. Schindler.  Garfield++, 2018.  URL https://garfieldpp.web.cern.ch/garfieldpp/.
  * [17] Stephen Biagi.  Magboltz-transport of electrons in gas mixtures.  Available at https://magboltz.web.cern.ch/magboltz/.
  * Jaegle et al. [2019] I. Jaegle et al.  Compact, directional neutron detectors capable of high-resolution nuclear recoil imaging.  _Nucl. Instrum. Meth. A_ , 945:162296, 2019.  doi: 10.1016/j.nima.2019.06.037.
  * Paszke et al. [2019] Adam Paszke et al.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach et al., editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.  URL http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * LeNail [2019] Alexander LeNail.  Nn-svg: Publication-ready neural network architecture schematics.  _Journal of Open Source Software_ , 4(33):747, 2019.  doi: 10.21105/joss.00747.  URL https://doi.org/10.21105/joss.00747.
  * Ioffe and Szegedy [2015] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift, 2015.
  * Klambauer et al. [2017] Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter.  Self-Normalizing Neural Networks.  _arXiv e-prints_ , art. arXiv:1706.02515, June 2017.
  * Srivastava et al. [2014] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15(56):1929–1958, 2014.  URL http://jmlr.org/papers/v15/srivastava14a.html.
  * Kingma and Ba [2017] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Rumelhart et al. [1986] David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning Representations by Back-propagating Errors.  _Nature_ , 323(6088):533–536, 1986.  doi: 10.1038/323533a0.  URL http://www.nature.com/articles/323533a0.
  * Prechelt [2012] Lutz Prechelt.  _Early Stopping — But When?_ , pages 53–67.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2012.  ISBN 978-3-642-35289-8.  doi: 10.1007/978-3-642-35289-8_5.  URL https://doi.org/10.1007/978-3-642-35289-8_5.
  * Hahsler et al. [2019] Michael Hahsler, Matthew Piekenbrock, and Derek Doran.  dbscan: Fast density-based clustering with R.  _Journal of Statistical Software_ , 91(1):1–30, 2019.  doi: 10.18637/jss.v091.i01.
  * Chen and Guestrin [2016] Tianqi Chen and Carlos Guestrin.  XGBoost: A scalable tree boosting system.  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pages 785–794, New York, NY, USA, 2016. ACM.  ISBN 978-1-4503-4232-2.  doi: 10.1145/2939672.2939785.  URL http://doi.acm.org/10.1145/2939672.2939785.
  * Quinlan [1986] J. R. Quinlan.  Induction of decision trees.  1(1):81–106, mar 1986.  ISSN 0885-6125.  doi: 10.1023/A:1022643204877.  URL https://doi.org/10.1023/A:1022643204877.


# High Resolution non-Markovianity in NMR

Nadja K. Bernardes1,∗, John P. S. Peterson2, Roberto S. Sarthour2, Alexandre
M. Souza2, C. H. Monken1, Itzhak Roditi2, Ivan S. Oliveira2, Marcelo F.
Santos1 1Departamento de Física, Universidade Federal de Minas Gerais, Belo
Horizonte, Caixa Postal 702, 30161-970, Brazil  2Centro Brasileiro de
Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, Rio de Janeiro, 22290-180,
Brazil<EMAIL_ADDRESS>
###### Abstract

Memoryless time evolutions are ubiquitous in nature but often correspond to a
resolution-induced approximation, i.e. there are correlations in time whose
effects are undetectable. Recent advances in the dynamical control of small
quantum systems provide the ideal scenario to probe some of these effects.
Here we experimentally demonstrate the precise induction of memory effects on
the evolution of a quantum coin (qubit) by correlations engineered in its
environment. In particular, we design a collisional model in Nuclear Magnetic
Resonance (NMR) and precisely control the strength of the effects by changing
the degree of correlation in the environment and its time of interaction with
the qubit. We also show how these effects can be hidden by the limited
resolution of the measurements performed on the qubit. The experiment
reinforces NMR as a test bed for the study of open quantum systems and the
simulation of their classical counterparts.

## Introduction

Markov processes are defined as those in which the future of a system bears no
correlation to its past. Such behavior is found in the description of
phenomena as far apart as nuclear fission, the decay of an electronically
excited atom, the growth of bacteria colonies or rabbit populations and the
compound interests that dictate our debts, to name a few kampen . In many
situations, however, these time evolutions prove to be no more than a
resolution-induced approximation, i.e. the interaction with external degrees
of freedom (here called environment) produces correlations in the dynamics,
but in a way that is too weak to be observed. Most of the times, such small
deviations are either harmless and/or useless. However, there are meaningful
counterexamples of both cases. For instance, many modern applications, such as
the generation of large cryptographic keys or statistical sampling, rely on
random number generators chaitin . Testing the independence of events in these
devices is essential to guarantee the safety or fairness of these applications
where even the tinniest correlations may be exploited. At the same time,
controllable non-Markovian evolutions in small quantum systems have been
thought of as promising tools for quantum information processing vasile ;
matsuzaki ; chin ; elsi ; Bogna .

Markovian processes exist both in discrete time steps or as a continuous time
evolution. A broadly used tool to discretize continuous Markovian dynamics are
the so-called collisional models. In these models, the continuous stochastic
process is replaced by a stroboscopic sequence formed by a series of discrete
time steps defined, each one, by a collision between the system of interest
and one of the particles forming its environment. The most famous application
of collisional models in physics is the description of Brownian motion. In
open quantum systems, collisional models have been extensively used as a test
bed to investigate the details of system-environment interactions and the
conditions under which non-Markovian dynamics arises rau ; ziman1 ; ziman2 ;
giovannetti ; tomas ; ciccarello ; vacchini ; budini ; Paternostro ; nadja1 .
Provided system and environment are initially uncorrelated, the evolution of
the system $\rho$ after $j$ collisions with the environment $\rho_{env}$ is
given by

$\rho\rightarrow\rho^{\prime}=\text{Tr}_{env}\left[U_{j}U_{j-1}...U_{2}U_{1}(\rho\otimes\rho_{env})U_{1}^{\dagger}U_{2}^{\dagger}...U_{j-1}^{\dagger}U_{j}^{\dagger}\right],$
(1)

where $\text{Tr}_{env}$ denotes the partial trace over the environmental
degrees of freedom and $U_{i}$ is a global unitary operation due to the $i$-th
collision. These models are simple to describe, easy to compute and have been
recently simulated in optical setups where one degree of freedom of the photon
is used as the quantum system and extra ones play the role of the environment
Liu ; Liu2 ; tang ; Steve1 ; Fabio2 ; chiuri ; Xu ; Fabio1 ; nadja2 .

In this paper, we implement a collisional model in Nuclear Magnetic Resonance
(NMR) to investigate with high resolution the subtle non-Markovian evolution
of a qubit subjected to the action of a structured and fully quantized
external environment. This experimental technique has been very successful in
the study of quantum information processes and in proofs of quantum principles
as may be seen on t1 ; t2 ; t3 ; t4 . Here, the environment is composed of
independent qubits that are initially correlated among themselves and collide,
one at a time, with the target. We characterize non-Markovianity by
calculating the Trace Distance
$D(\rho_{1}(t),\rho_{2}(t))=\frac{1}{2}\text{Tr}|\rho_{1}(t)-\rho_{2}(t)|$
between two different states of the qubit as a function of time. It is well
established that in a quantum stochastic evolution,
$D(\rho_{1}(t+\tau),\rho_{2}(t+\tau))>D(\rho_{1}(t),\rho_{2}(t))$ implies non-
Markovian behavior BLP . The experiment features only two time steps because
these are enough to establish all the properties we want to discuss. We
measure $D(\rho_{1}(t),\rho_{2}(t))$ after each collision as a function of the
degree of correlation in the initial state of the environment and also as a
function of the duration of each collision. Both quantities influence how much
the dynamics deviates from a Markov process.

## Results

We encoded the system-environment state in a sample of trifluoroiodoethylene
(C2F3I). The NMR experiment was performed using a 500 MHz Varian spectrometer
with the used sample diluted ($\sim 1\%$) in deuterated acetone (containing
$97\%$ of deuterium). A single molecule of C2F3I contains three atoms of
fluorine with nuclear spin-$1/2$, each of them representing a qubit. One of
the atoms will be regarded as the system and the other two as the environment.
These spins interact with each other and with applied magnetic fields as well.
We designed the experiment to implement the operations represented in the
quantum circuit of Fig. 1: the first qubit, top line in the circuit,
represents the system and the other lines are the qubits of the environment.
The two qubits of the environment are initially prepared in a correlated state
$\rho_{env}=\frac{q}{2}(|00\rangle\langle 00|+|11\rangle\langle
11|)+\frac{1-q}{2}(|01\rangle\langle 01|+|10\rangle\langle 10|)$ where “0”
means spin up and “1” means spin down. In each collision, the state of the
system undergoes a random-walk type of evolution, changing conditionally to
the environment: environmental spin up induces a rotation in a certain
direction (say, $x$), whereas environmental spin down induces a rotation in
the orthogonal direction ($y$). This situation can be summed up by the unitary
transformation $U=e^{i\eta\sigma_{x}}\otimes|0\rangle\langle
0|+e^{i\eta\sigma_{y}}\otimes|1\rangle\langle 1|$ where $|j\rangle\langle j|$
represents the internal state of an environmental particle and the
exponentials dictate the system’s rotations.

Figure 1: a) The nuclear spins of the atoms of fluorine of the C2F3I molecule
are the qubits. Two of them represent the environment and the other one is the
system. Their interaction is showed in this pictorial representation of our
collisional model. The action of the environment on the system’s qubit is
represented by the two boxes and the demon. In the first collision, the demon
will toss a coin and will decide if he will operate either
$e^{i\eta\sigma_{x}}$ or $e^{i\eta\sigma_{y}}$ on the system $\rho(0)$.
Afterwards he can still decide if the operation done in the second collision
is correlated or not with the previous one, resulting in the state $\rho(2)$.
b) A quantum circuit that describes the experiment. The environment is
prepared in state $\rho_{env}$ and will serve as the control qubit for the
controlled operation that happen on the state of the system $\rho$. After that
a quantum state tomography is performed on the first particle in order to
determine the state of the system.

In a NMR experiment the signal of each fluorine spin can be distinguished and
singled manipulated. The natural Hamiltonian of the total system may be
described as

$\mathcal{H}=\sum_{n}\hslash\left(\omega_{0n}-\omega_{r}\right)I_{z}^{n}+\sum_{k\neq
m}\hslash 2\pi\mathcal{J}_{km}I_{z}^{k}\otimes I_{z}^{m},$ (2)

where $\omega_{0n}$ and $\omega_{r}$ are the natural resonance frequency and
the frequency of the rotating reference frame of the $n$-th spin,
respectively. The first term represents the interaction of the spins with the
static magnetic field applied along the $z$ direction and the second term
represents the interaction between the three spins, with $\mathcal{J}_{km}$
being the exchange integral. Here, $I_{\alpha}^{i}$ is the spin operator for
the $i$-th spin, which for spin-$1/2$ corresponds to the Pauli matrix
($\sigma_{\alpha}$) divided by 2. The physical parameters of our molecule are
given in Methods. Since the individual chemical shifts of each spin place them
far apart, i.e. their resonance frequencies are separated enough
($\left|\omega_{0k}-\omega_{0m}\right|\gg\mathcal{J}_{km}$), the Ising
coupling approximation was considered here key-2 . In the NMR experiment the
spins also interact with oscillating magnetic fields that can be turned on and
off at any frequency and different amplitudes, named radio frequency pulses
($rf$). These are responsible for the system control, and its interaction can
de described by the Hamiltonian

$\mathcal{H}_{rf}=\hslash\omega_{1}(t)\left[I_{x}^{n}\cos\phi+I_{y}^{n}\sin\phi\right],$
(3)

where $\omega_{1}(t)$ defines the pulse modulation and duration, and $\phi$
stands for its phase. The radio frequency pulses are used to perform rotations
on the individual spins and by properly choosing these rotations, the quantum
circuit from Fig. 1 can be implemented, as described in Methods.

The experiment was performed for two different initial states of the system,
$\rho_{1}(0)=|0\rangle\langle 0|$ and $\rho_{2}(0)=|1\rangle\langle 1|$
(eigenstates of $\sigma_{z}$) and for the environment prepared in different
states with different degrees of correlation $q$ ($q=0$, 0.15, and 0.25). The
theoretical analysis shows that for $q$ equals $0$ and $0.15$ these are the
states that maximize $\Delta D$ for the chosen dynamics. For $q=0.25$ there is
a better pair of initial states but the eventual gain in $\Delta D$ (a factor
of $\sim 3$) still lies within the current error bar. Using the $rf$ pulses,
the first and the second collisions were implemented, resulting in the system
state $\rho_{1}(1)$ ($\rho_{2}(1)$) and $\rho_{1}(2)$ ($\rho_{2}(2)$) ,
respectively. After each run, full state tomography key-1  was performed in
the system.

Figure 2 shows the change in distance between the two initial states of the
system after one and two collisions $\Delta
D=D(\rho_{1}(2),\rho_{2}(2))-D(\rho_{1}(1),\rho_{2}(1))$, as a function of the
strength of each collision $\eta$ and for different degrees of correlation $q$
of the environmental state. For large enough interactions (larger $\eta$) and
anti-correlation in the environment (smaller $q$), the collisions clearly
generate a non-Markovian dynamics in the system ($\Delta D>0$). The
phenomenon, witnessed by the increase in $\Delta D$ also reflects a backflow
of information from the environment to the system as time progresses.

It is also worth noticing that, as the non-Markovian effects induced by the
collisions get smaller, either due to weaker collisions (small values of
$\eta$) and/or environmental correlations ($q\rightarrow 1/2$), the
uncertainties in $\Delta D$ due to the lack of resolution of the measurements
become of the same order of $\Delta D$. This is true even for the strongest
anti-correlation possible, $q=0$. In this case, $\Delta D$, a sufficient but
not necessary condition for non-Markovianity, does not provide a conclusive
answer anymore, even though theory still predicts its vanishing only
asymptotically with both parameters $\eta\rightarrow 0$ or $q\rightarrow 1/2$.
Our resolution to conclusively guarantee non-Markovianity was limited to
$\Delta D\sim 3.5\times 10^{-4}$.

Figure 2: The trace distance,
$D(\rho_{1}(t),\rho_{2}(t))=\frac{1}{2}\text{Tr}|\rho_{1}(t)-\rho_{2}(t)|$,
was employed to witness non-Markovianity. The increasing of
$D(\rho_{1}(t),\rho_{2}(t))$ with time is a trace of non-Markovianity. a) The
plots show the dependence of $\Delta D$ with the strength of each collision
$\eta$ for different values of the correlation parameter $q$ (0 for the solid
blue line, 0.15 for the dashed orange line and 0.25 for the dash-dotted green
line). $q=0$ means that the two collisions were totally anti-correlated. b)
Detail of the $\Delta D$ curve for $\eta<0.01$ and $q=0$. In both plots the
error bars were estimated taking into account every possible error source, as
discussed in Ref. cesar . This led to the error determining the Bloch vector
of the quantum state of $\delta r=5\times 10^{-4}$, independent of the
direction. By propagating this error, the error in
$D(\rho_{1}(t),\rho_{2}(t))$ is given by $\delta D=\delta r/\sqrt{2}$.

## Discussions

Non-markovian dynamics of collisional models has been experimentally
investigated in linear optical setups Liu ; Liu2 ; tang ; Steve1 ; Fabio2 ;
chiuri ; Xu ; Fabio1 ; nadja2  where the environment is simulated by semi-
classical degrees of freedom. Usually, such experiments are performed in the
photon counting regime and strong system-environment coupling. Our results put
in evidence, in a controllable way, the role played by environmental
correlations in the non-Markovian dynamics of quantum systems and also show
how such effects may be hidden by the external noise brought in by unavoidable
measurement uncertainties. We have shown how NMR can be exploited to study
quantum dynamics of open systems using a bona fide quantum environment. NMR
also presents a great advantage over previously used platforms when it comes
to the versatility to design and control different environmental states and
the precision to probe very small memory effects in the dynamics of quantum
systems. For example, here we have observed variations of the Trace Distance
of the order of $0.05\%$, a gain of at least two orders of magnitude over the
typical optical experiments based on measurements of the polarization of
photons. In photon polarization measurements, the signal-to-noise ratio is
less than $\sqrt{N}$, where $N$ is the number of photons in each measurement.
Such gain in precision can be essential, for example, if one would like to use
a similar setup to test very small time correlations of black boxes such as
quantum or classical random number generations, cryptographic machines and
alike. Furthermore, one can also explore the capacity to dynamically change
the environment during the interaction with the system and to directly
interfere with the system during its time evolution in order to test and
investigate very recent theoretical results in the field, such as the ones
respectively published in ciccarello  and simone .

## Methods

### Model

We consider a qubit $\rho$ that interacts with an environment from which it is
initially decoupled. The interaction consists of consecutive collisions with
qubits of the environment. During each collision, the system undergoes with
equal probability an evolution $e^{i\eta\sigma_{y}}$ or $e^{i\eta\sigma_{x}}$
controlled by the internal state of the corresponding environmental particle.
The map that describes the state of the system after collision one is
$\Lambda_{10}(\cdot)=\frac{1}{2}[e^{i\eta_{1}\sigma_{x}}(\cdot)e^{-i\eta\sigma_{x}}+e^{i\eta_{1}\sigma_{y}}(\cdot)e^{-i\eta_{1}\sigma_{y}}]$.
All the effects we are interested in can be observed with only two collisions
described, in our case, by the general map
$\Lambda_{20}(\cdot)=\frac{q}{2}[e^{i\eta_{2}\sigma_{x}}e^{i\eta_{1}\sigma_{x}}(\cdot)e^{-i\eta_{1}\sigma_{x}}e^{-i\eta_{2}\sigma_{x}}+e^{i\eta_{2}\sigma_{y}}e^{i\eta_{1}\sigma_{y}}(\cdot)e^{-i\eta_{1}\sigma_{y}}e^{-i\eta_{2}\sigma_{y}}]+\frac{1-q}{2}[e^{i\eta_{2}\sigma_{y}}e^{i\eta_{1}\sigma_{x}}(\cdot)e^{-i\eta_{1}\sigma_{x}}e^{-i\eta_{2}\sigma_{y}}+e^{i\eta_{2}\sigma_{x}}e^{i\eta_{1}\sigma_{y}}(\cdot)e^{-i\eta_{1}\sigma_{y}}e^{-i\eta_{2}\sigma_{x}}]]$,
where $q$ defines their degree of correlation: for $q=1$ ($q=0$) the
collisions are totally correlated (anti-correlated). One possible way to
generate this evolution is to consider two qubits as the environment prepared
in the following state $\rho_{env}=\frac{q}{2}(|00\rangle\langle
00|+|11\rangle\langle 11|)+\frac{1-q}{2}(|01\rangle\langle
01|+|10\rangle\langle 10|)$. (Note that a classically correlated state or the
state
$\sqrt{\frac{q}{2}}(|00\rangle+|11\rangle)+\sqrt{\frac{1-q}{2}}(|01\rangle+|10\rangle)$
will render the same system evolution). The interaction in each collision is
described by a unitary transformation given by
$U_{\eta}=e^{i\eta\sigma_{x}}\otimes|0\rangle\langle
0|+e^{i\eta\sigma_{y}}\otimes|1\rangle\langle 1|$, which works as a
conditional operation. The resulting state after one collision will be
$\rho(1)=\Lambda_{10}(\rho(0))=\text{Tr}_{env}\left[(U_{\eta_{1}}\otimes\mathds{1})(\rho(0)\otimes\omega_{env})(U_{\eta_{1}}^{\dagger}\otimes\mathds{1})\right]=\frac{1}{2}(e^{i\eta_{1}\sigma_{x}}\rho(0)e^{-i\eta_{1}\sigma_{x}}+e^{i\eta_{1}\sigma_{y}}\rho(0)e^{-i\eta_{1}\sigma_{y}})$,
where $\text{Tr}_{env}$ denotes the partial trace over the environmental
degrees of freedom. And after the second collision the evolution of the system
reads
$\rho(2)=\Lambda_{20}(\rho(0))=\text{Tr}_{env}\left[U_{\eta_{2}}U_{\eta_{1}}(\rho(0)\otimes\omega_{env})U^{\dagger}_{\eta_{1}}U^{\dagger}_{\eta_{2}}\right]=\frac{q}{2}[e^{i\eta_{2}\sigma_{x}}e^{i\eta_{1}\sigma_{x}}(\rho(0))e^{-i\eta_{1}\sigma_{x}}e^{-i\eta_{2}\sigma_{x}}+e^{i\eta_{2}\sigma_{y}}e^{i\eta_{1}\sigma_{y}}(\rho(0))e^{-i\eta_{1}\sigma_{y}}e^{-i\eta_{2}\sigma_{y}}]+\frac{1-q}{2}[e^{i\eta_{2}\sigma_{y}}e^{i\eta_{1}\sigma_{x}}(\rho(0))e^{-i\eta_{1}\sigma_{x}}e^{-i\eta_{2}\sigma_{y}}+e^{i\eta_{2}\sigma_{x}}e^{i\eta_{1}\sigma_{y}}(\rho(0))e^{-i\eta_{1}\sigma_{y}}e^{-i\eta_{2}\sigma_{x}}]]$.
For simplicity, we chose these two specific operators, however, similar
results would be observed had we used two different non-commuting operators.
The parameters $\eta_{1}$ and $\eta_{2}$ are related to the strength of
interaction and the time for one collision. We consider
$\eta_{1}=\eta_{2}=\eta$, and considering a weak interaction and a short
collision time, $\eta$ are assumed to be sufficiently small ($\eta<0.1$) such
that $e^{i\eta\sigma_{x,y}}\sim(1+i\eta\sigma_{x,y})$.

In the first run of the experiment, the system interacts (collides) with qubit
one of the environment and, in the second run, the system repeats this
collision and subsequently collides with qubit two of the environment. After
each run, quantum state tomography key-1  is performed on the system in order
to fully determine its state. Two orthogonal initial states of the system,
$\rho_{1}(0)=\left|0\right\rangle\left\langle 0\right|$ and
$\rho_{2}(0)=\left|1\right\rangle\left\langle 1\right|$, are prepared so that
the distance $D(\rho_{1}(j),\rho_{2}(j))$ between them can be calculated after
each run $j$ (orthogonality sets $D(\rho_{1}(0),\rho_{2}(0))=1$).

Being a qubit, the state of the the system can also be represented by a vector
$\vec{r}=(\langle\sigma_{x}\rangle,\langle\sigma_{y}\rangle,\langle\sigma_{z}\rangle)$
in the Bloch Sphere
($\langle\sigma_{i}\rangle=\textrm{Tr}\\{\sigma_{i}\rho\\}$) in which case the
distance between $\rho_{1}$ and $\rho_{2}$ reduces to the geometric distance
between the respective vectors
$D(\rho_{1},\rho_{2})=\frac{1}{2}|\vec{r}_{1}-\vec{r}_{2}|$. This quantity can
be easily calculated for the initial states $\rho_{1}(0)=|0\rangle\langle 0|$
and $\rho_{2}(0)=|1\rangle\langle 1|$ and for the rotations used in the
experiment and gives
$D(\rho_{1}(1),\rho_{2}(1))=\frac{1}{2}[3+\cos(4\eta)]^{\frac{1}{2}}$ and
$D(\rho_{1}(2),\rho_{2}(2))=\frac{1}{2}\left\\{16[\cos\eta\sin\eta(\cos^{2}\eta-q\sin^{2}\eta)]^{2}+[(q+1)\cos(4\eta)-q+1]^{2}+4\sin^{2}\eta\cos^{2}\eta[(q+1)\cos(2\eta)-q+1]^{2}\right\\}^{\frac{1}{2}}$
after the first and second collision, respectively. Expanding these results
for $\eta\ll 1$ we get $\Delta D=(1-4q)\eta^{2}+O(\eta^{4})$, making it clear
why, for the range of $\eta$ used in the experiment, $\Delta D$ must be
positive, i.e. the dynamics must be non-Markovian, for any $q<1/4$ (highly
anti-correlated environmental qubits). This result corresponds to the solid
curves in Fig. 2.

### Experimental details

The NMR experiment was performed encoding the system of interest and
environment in a trifluoroiodoethylene (C2F3I) sample diluted ($1\%$)
deuterated acetone (containing $97\%$ of deuterium). The three nuclei
spin-$1/2$ of the fluorine atoms were associated to system and environment
qubits. The total system is described by the Hamiltonian

$\mathcal{H}=\sum_{n}\hslash\left(\omega_{0n}-\omega_{r}\right)I_{z}^{n}+\sum_{k\neq
m}\hslash 2\pi\mathcal{J}_{km}I_{z}^{k}\otimes I_{z}^{m},$ (4)

where $\omega_{0n}$ and $\omega_{r}$ are the natural resonance frequency and
the frequency of the rotating reference frame of the $n$-th spin, $I_{z}^{n}$
is the $z$-component of the spin angular moment of nucleus $n$. The second
term represents the interaction - being $\mathcal{J}_{km}$ the exchange
integral - between the three spins. Fig. 3 shows the values of
$(\omega{}_{0n}-\omega_{r})$ and $\mathcal{J}_{km}$.

Figure 3: Table of $(\omega_{0k}-\omega_{r})/(2\pi)$ in the main diagonal
elements and $\mathcal{J}_{km}$ as the other elements, all values are given in
Hz.

### Preparation of the initial state and unitary operations

Figure 4: Quantum circuit representing the NMR pulse sequence for implementing
the first collision. The boxes with the symbols $R_{\alpha}(\theta)$ indicate
that a rotation of the angle $\theta$ on that particular spin was performed
around the $\alpha$ direction and
$R_{\bar{\alpha}}(\theta)=R_{\alpha}(-\theta)$. The boxes with the symbol
$\tau_{1}$ represent a free evolution of the system and environment for a
period of time $\tau_{1}$. After the first collision, the resulting quantum
state of system and environment is represented by $\rho_{t}(1)$. $\tau_{1}$ is
the time needed such that a rotation of $\pi/2$ occurs in the system and here
it will be given by $\tau_{1}=1/(4\mathcal{J}_{s,env1})=0.00358$ s. Figure 5:
Quantum circuit representing the NMR pulse sequence for implementing the
second collision. The input state is the resulting state from the first
collision $\rho_{t}(1)$. The boxes with the symbols $R_{\alpha}(\theta)$
indicate that a rotation of the angle $\theta$ on that particular spin was
performed around the $\alpha$ direction and
$R_{\bar{\alpha}}(\theta)=R_{\alpha}(-\theta)$. The boxes with the symbol
$\tau_{1}$ represent a free evolution of the system and environment for a
period of time $\tau_{1}$. $\tau_{1}$ is the time needed such that a rotation
of $\pi/2$ occurs in the system and here it will be given by
$\tau_{1}=1/(4\mathcal{J}_{s,env2})=0.00525$ s.

The NMR implementation of the collisions between the system and particles of
the environment are shown in Figs. 4,5. The boxes with the symbols
$R_{\alpha}(\theta)$ indicate that a rotation of the angle $\theta$ on that
particular spin was performed around the $\alpha$ direction. The big boxes are
free evolution of the system for the time $\tau$. The angles $\theta$ are
related to the parameter $\eta$ of the collisional model, and they are of the
order of 1°. For the steps used in the experiment the differences between them
were approximately of 0.5°. The system we have used is homonuclear, which
means that all of the spins are of the same species (fluorine in this case),
and implementing single rotations were hard because they are close in
frequency. Then, the length of the pulses for exciting only one spin (qubit)
are long. This causes undesired evolutions of the quantum state due to the
interactions between the spins, leading the whole system to evolve while the
individual operations are applied. Furthermore, the $rf$ pulses are also
imperfect and they may affect others spins, besides the one which it is
intended for. This means that at the end of the quantum algorithm lots of
errors due to these undesirable evolutions and pulse imperfections will have
been accumulated. In order to correct the state of the system the method
described in key-2  was used.

Figure 6: Quantum circuit representing the NMR pulse sequence for the state
preparation of the initial state of system and environment given by
$\rho(0)\otimes\rho_{env}$. The boxes with the symbols $R_{\alpha}(\theta)$
indicate that a rotation of the angle $\theta$ on that particular spin was
performed around the $\alpha$ direction and
$R_{\bar{\alpha}}(\theta)=R_{\alpha}(-\theta)$. The boxes with the symbol
$\tau_{2}$ represent a free evolution of the system and environment for a
period of time $\tau_{2}$. The box with the symbol Grad stands for the
magnetic field gradient that is applied in the $z$-direction.

For the experiment, the values of $\eta$ were very small and this created a
huge source of errors in the operations, since they had to be performed in
tiny steps and, therefore, were not so different. In fact, the rotations
angles were so small that the controlled operations, which simulate the
collisions, were very hard to be performed properly. In order to achieve the
necessary precision in the experiment the rotations of $\theta$ had to be
precisely adjusted for each value of $\eta$ of the specific collision. For
adjusting the angles values, two parameters may be varied, the amplitude and
duration of the radio frequency pulses. However, varying these two parameters
the sequence of operations needed for the corrections changes as well.
Therefore, the correct pulse sequences were determined by combining the two
analysis and testing in the spectrometer. A good precision could be achieved
then and the small rotations, of the order of 0.30 degrees, could be
implemented. For the state preparation the pulse sequence shown in Fig. 6 was
performed, for more details see key-3 ; key-4 . The time $\tau_{2}$ of the
free evolution here has a strict relation to the correlation parameter $q$ and
is presented in Fig. 7. A the end of the circuit a magnetic field gradient in
the $z$-direction is applied, which works as a transversal relaxation time
killing the off-diagonal terms of the density matrix of system and
environment. Finally, after each run, full state tomography was performed in
the state of the system key-1 .

Note that this experiment was done much faster ($\sim 0.03$ s) than the
typical relaxation times of the elements of the C2F3I molecule, which is of
the order of 1 s.

Figure 7: Parameter $\tau_{2}$ has a strictly relation to the correlation
parameter $q$. The values used in the experiment are presented here.

### Acknowledgments

We would like to thank the support from the Brazilian agencies CNPq and CAPES.
J.P.S.P., R.S.S, A.M.S, I.R, and I.S.O would like to thank the support of
FAPERJ. M.F.S. would like to thank the support of FAPEMIG, project PPM IV.
This work is part of the INCT-IQ from CNPq.

## References

  * (1) van Kampen, N. G. Stochastic Processes in Physics and Chemistry (Elsevier, Amsterdam, 2007).
  * (2) Chaitin, G. J. Information randomness and incompleteness: papers on algorithmic information theory (World Scientific, Singapore, 1990).
  * (3) Vasile, R., Olivares, S., Paris, M. G. A. and Maniscalco, S. Continuous-variable quantum key distribution in non-Markovian channels. Phys. Rev. A 83, 042321 (2011).
  * (4) Matsuzaki, Y., Benjamin, S. C. and Fitzsimons, J. Magnetic field sensing beyond the standard quantum limit under the effect of decoherence. Phys. Rev. A 84, 012103 (2011).
  * (5) Chin, A. W., Huelga, S. F. and Plenio, M. B. Quantum Metrology in Non-Markovian Environments. Phys. Rev. Lett. 109, 233601 (2012).
  * (6) Laine, E.-M., Breuer, H.-P. and Piilo, J. Nonlocal memory effects allow perfect teleportation with mixed states. Sci. Rep. 4, 4620 (2014).
  * (7) Bylicka, B., Chruscinski, D. and Maniscalco, S. Non-Markovianity and reservoir memory of quantum channels: a quantum information theory perspective. Sci. Rep. 4, 5720 (2014).
  * (8) Rau, J. Relaxation Phenomena in Spin and Harmonic Oscillator Systems. Phys. Rev. 129, 1880 (1963).
  * (9) Ziman M. and Buzek V. All (qubit) decoherences: Complete characterization and physical implementation. Phys. Rev. A 72, 022110 (2005).
  * (10) Ziman, M., Stelmachovic, P. and Buzek, V. Description of Quantum Dynamics of Open Systems Based on Collision-Like Models. Open Syst. Inf. Dyn. 12, 81 (2005).
  * (11) Giovannetti, V. and Palma, G. M. Master Equations for Correlated Quantum Channels. Phys. Rev. Lett. 108, 040401 (2012).
  * (12) Rybár, T., Filipov, S. N. , Ziman, M. and Buzek, V. Simulation of indivisible qubit channels in collision models. J. Phys. B: At. Mol. Opt. Phys. 45, 154006 (2012).
  * (13) Cicarello, F., Palma, G. M. and Giovannetti, V. Collision-model-based approach to non-Markovian quantum dynamics. Phys. Rev. A 87, 040103(R) (2013).
  * (14) Vacchini, B. Non-Markovian master equations from piecewise dynamics. Phys. Rev. A 87, 030101(R) (2013).
  * (15) Budini, A. A. Embedding non-Markovian quantum collisional models into bipartite Markovian dynamics. Phys. Rev. A 88, 032115 (2013).
  * (16) McCloskey, R. and Paternostro, M. Non-Markovianity and system-environment correlations in a microscopic collision model. Phys. Rev. A 89, 052120 (2014).
  * (17) Bernardes, N. K., Carvalho, A. R. R., Monken, C. H., Santos, M. F. Environmental correlations and Markovian to non-Markovian transitions in collisional models. Phys. Rev. A 90, 032111 (2014).
  * (18) Liu, B.-H. et al. Experimental control of the transition from Markovian to non-Markovian dynamics of open quantum systems. Experimental control of the transition from Markovian to non-Markovian dynamics of open quantum systems. Nat. Phys. 7, 931 (2011).
  * (19) Liu, B.-H. et al. Photonic realization of nonlocal memory effects and non-Markovian quantum probes. Sci. Rep. 3, 1781 (2013).
  * (20) Tang, J.-S. T et al. Measuring non-Markovianity of processes with controllable system-environment interaction. Europ. Phys. Lett. 97, 10002 (2012).
  * (21) Fanchini, F., et al. Non-Markovianity through accessible information. Phys. Rev. Lett. 112, 210402 (2014).
  * (22) Jin, J., et al. All-optical non-Markovian stroboscopic quantum simulator. Phys. Rev. A 91, 012122 (2015).
  * (23) Chiuri, A., Greganti, C., Mazzola, L., Paternostro, M. and Mataloni, P. Linear Optics Simulation of Quantum Non-Markovian Dynamics. Sci. Rep. 2, 968 (2012).
  * (24) Xu, J.-S., et al. Experimental recovery of quantum correlations in absence of system-environment back-action. Nat. Comm. 4, 2851 (2013).
  * (25) Orieux, A., et al. Experimental on-demand recovery of entanglement by local operations within non-Markovian dynamics. Sci. Rep. 5, 8575 (2015).
  * (26) Bernardes, N. K., Cuevas, A. , Orieux, A., Monken, C. H., Mataloni, P., Sciarrino, F. and Santos, M. F. Experimental observation of weak non-Markovianity. Sci. Rep. 5, 17520 (2015).
  * (27) Batalhao, T. B., Souza, A. M., Sarthour, R. S., Oliveira, I. S., Paternostro, M., Lutz, E. and Serra, R. M. Irreversibility and the Arrow of Time in a Quenched Quantum System. Phys. Rev. Lett. 115, 190601 (2015).
  * (28) Auccaise, R., Araujo-Ferreira, A. G., Sarthour, R. S., Oliveira, I. S., Bonagamba, T. J. and Roditi, I. Spin Squeezing in a Quadrupolar Nuclei NMR System. Phys. Rev. Lett. 114, 043604 (2015).
  * (29) Batalhao, T. B., Souza, A. M., Mazzola, L., Auccaise, R., Sarthour, R. S., Oliveira, I. S., Goold, J., De Chiara, G., Paternostro, M. and Serra, R. M. Experimental Reconstruction of Work Distribution and Study of Fluctuation Relations in a Closed Quantum System. Phys. Rev. Lett. 113, 140601 ( 2014).
  * (30) Anvari Vind, F. , Foerster, A., Oliveira, I. S., Sarthour, R. S., Soares-Pinto, D. O., Souza A. M. and Roditi, I. Experimental realization of the Yang-Baxter Equation via NMR interferometry. Sci. Rep. 6, 20789 (2016).
  * (31) Breuer, H.-P., Laine, E.-M. and Piilo, J. Measure for the Degree of Non-Markovian Behavior of Quantum Processes in Open Systems. Phys. Rev. Lett. 103, 210401 (2009).
  * (32) Leskowitz G. M. and Mueller, L. J. State interrogation in nuclear magnetic resonance quantum-information processing. Phys Rev. A 69, 052302 (2004).
  * (33) Ryan, C. A., Negrevergne, C., Laforest, M., Knill, E. and Laflamme, R. Liquid-state nuclear magnetic resonance as a testbed for developing quantum control methods. Phys. Rev. A 78, 012328 (2008).
  * (34) Gershenfeld N. A. and Chuang I. L. Bulk Spin-Resonance Quantum Computation. Science 275, 350 (1997).
  * (35) Cory, D. G., Fahmy, A. F. and Havel, T. F. Ensemble quantum computing by NMR?spectroscopy. Proc. Natl. Acad. Sci. USA 94, 1634 (1997).
  * (36) Mukherjee, V., Giovannetti, V., Fazio, R., Huelga, S. F., Calarco, T. and Montangero, S. Efficiency of quantum controlled non-Markovian thermalization. New J. Phys. 17 063031, (2015).
  * (37) Raitz, C., Souza, A. M., Auccaise, R., Sarthour, R. S. and Oliveira, I. S. Experimental implementation of a nonthermalizing quantum thermometer. Quantum Inf. Process 14, 37 (2015).