# Interpretable, calibrated neural networks for analysis and understanding of
inelastic neutron scattering data

Keith T. Butler1,2111email<EMAIL_ADDRESS>Manh Duc Le3, Jeyan
Thiyagalingam1,4, Toby G. Perring3 1SciML, Scientific Computing Department,
STFC Rutherford Appleton Laboratory, Harwell Campus, Didcot, OX11 0QX, UK
2Department of Materials Science and Engineering, University of Oxford, 21
Banbury Rd, Oxford OX2 6HT  3ISIS Neutron and Muon Source, STFC Rutherford
Appleton Laboratory, Harwell Campus, Didcot, OX11 0QX, UK  4Department of
Engineering Science, University of Oxford, Parks Road, Oxford, OX1 3PJ

###### Abstract

Deep neural networks provide flexible frameworks for learning data
representations and functions relating data to other properties and are often
claimed to achieve ‘super-human’ performance in inferring relationships
between input data and desired property. In the context of inelastic neutron
scattering experiments, however, as in many other scientific scenarios, a
number of issues arise: (i) scarcity of labelled experimental data, (ii) lack
of uncertainty quantification on results, and (iii) lack of interpretability
of the deep neural networks. In this work we examine approaches to all three
issues. We use simulated data to train a deep neural network to distinguish
between two possible magnetic exchange models of a half-doped manganite. We
apply the recently developed deterministic uncertainty quantification method
to provide error estimates for the classification, demonstrating in the
process how important realistic representations of instrument resolution in
the training data are for reliable estimates on experimental data. Finally we
use class activation maps to determine which regions of the spectra are most
important for the final classification result reached by the network.

##  1  Introduction

Neural networks have witnessed a renaissance in the last decade, dramatically
improving on existing state of the art performance in fields from image
processing to automatic language translation. In the sciences, it has been
proposed that we are on the cusp of a “fourth paradigm” of data-driven
discovery [1, 2], an idea that is supported by the recent explosion in
research using machine learning (ML), not least in materials science and
condensed matter physics [3, 4, 5, 6].

Experiments at large-scale facilities such as neutron sources would seem to be
a natural fit for the application of data-driven methods for analysis [7].
Indeed there have been a number of recent publications applying machine
learning to many aspects of neutron science, from pulse discrimination in
scintillator detectors [8], to enhancing resolution [9] and exploration of
collected data [10]. ML techniques have been used for the analysis of diffuse
neutron scattering in spin-ice systems [11], small angle neutron scattering
[12, 13] and for constructing inter-atomic potentials to compare to inelastic
neutron scattering data [14]. However, questions still remain to be answered
before machine learning becomes a fully integrated and trusted part of the
neutron data analysis workflow; for example can we place uncertainty estimates
on our predictions and can we explain why a neural network predicted what it
did? A particularly challenging case is the interpretation of inelastic
neutron scattering (INS) data from single crystal samples.

INS is a powerful technique to investigate the wave vector and frequency
(equivalently excitation energy) dependency of the spectrum of excitations in
condensed matter. The neutron scattering cross-section is particularly simple:
it directly yields two-particle correlation functions, for example the
particle-particle correlation function, or in the case of magnetic systems the
spin-spin correlation function, and by virtue of the fluctuation-dissipation
theorem, the generalised response functions of the system [15, 16].
Consequently INS has played a central role determining, for example, phonon
and spin wave dispersion relations and the spectra of magnetic fluctuations in
materials, which in turn can play central roles in the thermal and charge
transport properties, and the mechanisms of exotic phenomena such as high
temperature superconductivity and heavy fermion behaviour [17, 18, 19, 20].

The richest information comes from experiments with single crystal samples,
where the spectra can be explored as a function of all four components of wave
vector and frequency. The latest INS spectrometers for such experiments at
pulsed neutron sources (e.g. [21, 22, 23]) and equivalent instruments at
reactors (e.g. [24]) allow the incident neutron energy to be fixed, and
typically have $\approx$3 steradians of position sensitive area detectors
divided into $\sim 10^{5}$ detector elements, in each of which the energy of
every scattered neutron (and hence the energy transferred to excitations) is
resolved into $\sim 200-500$ channels. With just one crystal orientation,
therefore, data are collected in $>10^{7}$ voxels in a three dimensional (3D)
space (corresponding to the the two coordinates on the area detector, and
energy transfer). In the case of quasi-2D magnetic materials, where the
magnetic coupling is significant only within planes of atoms but negligible
between planes, these three coordinates can be transformed into those of wave
vector within the plane and energy transfer, thereby in principle enabling the
entire excitation spectrum to be measured in parallel. An extension of this
approach is to combine $\sim 100-500$ measurements between which the crystal
is rotated by a few tenths of a degree. The crystal rotation provides another
independent coordinate so that the full 4D wave vector - excitation energy
space of excitations can be measured, in $\sim 10^{9}-\sim 10^{10}$ voxels
([25] and references therein).

The sizes of the datasets collected on these instruments are challenging to
analyse and interpret with conventional but well-established and well-
understood methods and software. The motivation for this work has been to
investigate to what extent ML methods can reduce the labour of inelastic
neutron data analysis while still returning reliable and interpretable
information. In this contribution we re-examine the data from a well-
understood INS measurement of spin waves in a single crystal of
Pr(Ca0.9Sr0.1)2Mn2O7 [26], a moderately complex magnetic material, using
recent advances in deep learning. Different models of the exchange
interactions in Pr(Ca,Sr)2Mn2O7 were proposed based on the observed magnetic
and crystallographic structure, but INS measurements of the spin wave spectrum
and theoretical modelling were key to distinguishing between them. We show
that ML techniques can be used to distinguish between the models and can
highlight the diagnostic regions of the spectrum responsible for this
categorisation.

In Section 2 we briefly describe the interest in and theoretical background to
the physics of the manganites, and Pr(Ca0.9Sr0.1)2Mn2O7 in particular,
together with the different exchange models. In Section 3 we discuss the
experimental data and how it was modelled using linear spin wave theory to
create training datasets for our neural networks, taking account of the
instrument resolution. In Section 4, the main part of the paper, we introduce
the neural network classifier, incorporate uncertainty quantification into the
network and show how class activation maps may be used to understand its
inferences. Finally, we discuss future challenges and summarise our
conclusions.

##  2 Exchange interaction models in half-doped manganites

Perovskite manganese oxides RE(1-x)AxMnO3 (RE = rare earth, A = Ca,Sr,Ba,Pb)
and their layered analogues have been widely studied not just because of the
colossal magnetoresistance (CMR) they can show - up to many orders of
magnitude change in resistivity in applied magnetic fields of a few Tesla -
but also because of the rich physics that arises from the coupling of the
charge, spin, lattice and orbital degrees of freedom [27]. Typically the most
pronounced CMR is with light hole doping $x\approx 0.2-0.4$ where the material
is balanced between being in a ferromagnetic (FM) metallic phase and an
antiferromagnetic (AFM) charge and orbital ordered (CO) insulating phase. One
of the commonest AFM CO phases is found near half doping, $x=0.5$, and the
nature of this phase was the topic of the INS study in ref. [26].

The structure of the perovskite manganese oxides consists of manganese ions,
each surrounded by six oxygen ions forming a distorted octahedron, and these
octahedra are in turn connected at their vertices to create a (nearly) cubic
network, with a rare earth or alkaline earth cation filling the gap between
the octahedra. Pr(Ca,Sr)2Mn2O7 is a half-doped bilayer manganite consisting of
pairs of MnO6 octahedral layers, which are separated by a layer of the
interstitial atoms that reduce the magnetic exchange coupling between bilayers
by two orders of magnitude compared to the intra-bilayer couplings.

The famous and long-standing _Goodenough_ model for the AFM CO phase is shown
in Figure 1(a). It is based on the Goodenough-Kanamori-Anderson rules [28, 29]
for magnetic superexchange together with a checkerboard ordering of equal
numbers of formally $3+$ and $4+$ valence Mn ions. The octahedra around the
Mn3+ are elongated due to the singly occupied higher energy $e_{g}$ orbitals.
Accommodating this elongation within the structure results in a herringbone
pattern of the $e_{g}$ orbitals and the experimentally observed CE-type
magnetic structure, consisting of ferromagnetic ordering along zig-zag chains
with antiferromagnetic ordering between neighbouring zig-zags.

A very different and somewhat controversial model came to prominence in 2002,
following a seminal neutron diffraction experiment [30]. In this model, known
as the Zener polaron model, adjacent Mn3+ and Mn4+ spins form dimers with
strong ferromagnetic intra-dimer coupling, but much weaker inter-dimer
coupling (Figure 1(b)). The two models have the same underlying periodicities
and it is difficult to distinguish between them using diffraction techniques
as this relies on the presence or absence of (and the differences between the
intensities of) weak superlattice Bragg peaks. However, measuring the spin
wave spectrum offers a way to distinguish between the models. The limiting
case of strong intra-dimer coupling was straightforwardly eliminated on the
grounds of the periodicities of the spin wave dispersion relations [26]
compared to the data. A more realistic scenario is that of weaker
dimerisation. This _dimer_ model was eventually eliminated in ref. [26], but
it was far from straightforward, because with an appropriate set of magnetic
exchange parameters the dimer model can closely reproduce the spin wave
dispersion relations in the Goodenough model, with only subtle differences in
the intensities throughout most of the Brillouin zone. Determining which of
the Goodenough and dimer models correctly described the magnetic interactions
was the primary result of ref. [26] which concluded that the Goodenough model
best fits the INS data.

Figure 1:  Two magnetic exchange models in a single sheet of Mn ions in a
half-doped manganite showing the CE-type magnetic ordering. (a) Goodenough
model (b) dimer model, similar to Goodenough but neighbouring Mn3+/Mn4+ sites
along the ferromagnetic zig-zag chains are loosely bound in a dimer.

Both the Goodenough and dimer models can be described by a spin Hamiltonian
taking the standard Heisenberg form

$\mathcal{H}=\sum_{i,j}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+D\sum_{i}(S_{i}^{z})^{2}$
(1)

where the atom pairs $ij$ are not the same in the two models, so the topology
of the exchange interactions $J_{ij}$ will differ. In this formulation, a
positive $J_{ij}$ indicates antiferromagnetic interactions (this is opposite
to the convention used in ref [26]). The final term is a single-ion anisotropy
term which with $D>0$ tends to keep the spins in the $a-b$ plane consistent
with observations. Figure 1 shows the exchange topology of each model, and
Table 1 summarises the published exchange interactions [26] deduced for the
two models.

In the Goodenough model the nearest-neighbour FM interaction $J_{F1}$ along
the zig-zags is supplemented by two next-nearest-neighbour interactions
$J_{F2}$ and $J_{F3}$, which need not be FM but must be permitted to reproduce
experimentally observed periodicities in the experimental data. In addition
there is an AFM nearest-neighbour exchange $J_{A}$ between the zig-zags. In
the dimer model, the nearest-neighbour FM interactions within a zig-zag chain
are replaced by a stronger intra-dimer FM interaction $J_{FS}$ and a weaker
inter-dimer FM interaction $J_{FW}$. Between the zig-zag chains are additional
inter-dimer interactions denoted $J_{A1}$ and $J_{A2}$ as they link
antiparallel spins. The symmetry of the model does not require them to be
identical, nor is it necessary for both of them to be AFM. Finally, in both
models there is an inter-layer interaction coupling the two sheets of the
bilayer along the crystallographic $c$ direction, denoted $J_{\perp}$. The
magnetic coupling between bilayers is negligible; consequently the spin wave
spectrum is dispersive only in the $a-b$ plane and so can be treated as
quasi-2D.

Goodenough model | dimer model
---|---
| Ref (meV) | Range (meV) |  | Ref (meV) | Range (meV)
$J_{F1}$ | -11.39(5) | [-20, 0] | $J_{FS}$ | -14.20(8) | [-20, 0]
$J_{A}$ | 1.50(2) | [0, 3] | $J_{FW}$ | -8.43(6) | [-20, 0]
$J_{F2}$ | -1.35(7) | [-3, 3] | $J_{A1}$ | 1.52(1) | [-3, 3]
$J_{F3}$ | 1.50(5) | [-3, 3] | $J_{A2}$ | 1.52(1) | [-3, 3]
$J_{\perp}$ | 0.88(3) | [0, 3] | $J_{\perp}$ | 0.92(3) | [-20, 0]
$D$ | 0.074(1) | [0, 0.2] | $D$ | 0.073(1) | [0, 0.2]
Table 1:  Spin wave exchange parameters (meV) for the Goodenough and dimer
models for Pr(Ca0.9Sr0.1)2Mn2O7 with standard errors in parentheses as
determined by [26], and the range of values used to generate random datasets
for training the neural network classifier in this work. Positive values
indicate antiferromagnetic exchange. $D$ denotes the single-ion anisotropy
term, with positive values here indicating an easy-plane anisotropy.

##  3 Data

###  3.1 Experimental data

The experimental INS data was acquired using the MAPS spectrometer at the ISIS
Neutron and Muon Source, and the conventional analysis of that data was
previously published [26]. Datasets were collected with several different
incident neutron energies and monochromating (Fermi) chopper speeds
($E_{i}=25$, 35, 50, 70, 100, 140 meV, with corresponding speeds $f=300$, 200,
200, 250, 300, 400 Hz) to span the full bandwidth of the spin wave spectrum in
a series of energy transfer windows with appropriate wave vector and energy
resolution. The co-aligned array of single crystals of Pr(Ca0.9Sr0.1)2Mn2O7
was mounted with the incident beam parallel to the $c$ axis so that the $a-b$
plane is imaged in the detectors. All measurements were carried out in a
closed cycle refrigerator at 4 K. For each dataset, an estimate of non-
magnetic scattering and background was made from several cuts as a function of
energy transfer ($\hbar\omega$) at different wave vectors, from which sections
were stitched together to construct a single cut that included no spin wave
scattering.

Figure 2:  2D representation of the experimental data (left panel). The data
are arranged column-wise in terms of incident neutron energy ($E_{i}$), with
$E_{i}$ in meV given at the top of the column. The data are arranged row-wise
into bins of energy transfer $\hbar\omega$=0.10-0.16$E_{i}$, 0.16-0.22$E_{i}$
etc. The middle and right panels show the calculated spectra using the
Goodenough (middle) and dimer (right) exchange models at the published
parameters noted in Table 1.

To use with the convolutional neural networks described below, the background-
subtracted data were arranged into a gray scale image, obtained by first
vertically stacking ten constant-energy transfer slices per incident energy
spectrum, each stack uniformly dividing the energy transfer range between
0.1$E_{i}$ and 0.7$E_{i}$ into slices with thickness 0.06$E_{i}$. We chose to
exclude data close to the elastic line (i.e. zero energy transfer) as the
signal there is dominated by non-magnetic scattering, and also to exclude the
spectrum above 0.7$E_{i}$ as the signal there has a significant non-magnetic
background. The constant-energy transfer slices are gray scale 40$\times$40
pixel images over the wave vector range $-1<h,k<1$, which covers the full
periodicity of the spin waves in the reciprocal lattice. The ten slices for
each of the six $E_{i}$ form 40$\times$400 pixel strips, which were then
stacked horizontally producing a 240$\times$400 pixel image, shown in Figure
2. This image was used as input to the trained neural network classifier in
order to determine which exchange model (Goodenough or dimer) better describes
the measured data.

###  3.2 Training Data

The convolutional neural networks described in the next section were trained
on simulated data, in the same 240$\times$400 gray scale image format as the
experimental data. For the training datasets we used linear spin wave theory,
as implemented in the _SpinW_ code [31], to simulate the spin wave spectrum
corresponding to sets of distinct spin Hamiltonian parameters
$\\{J_{ij},D\\}$. Each constant energy transfer slice (40$\times$40 pixel sub-
image) is computed independently by convolving the spin wave theory spectrum
with the spectrometer resolution function for each detector element - energy
transfer bin in the data, using a Monte Carlo integration method [32] (named
“Tobyfit” after the first program to implement it), as implemented in the
neutron scattering data analysis code _Horace_ [25]. The resulting simulation
is sliced to produce a stack of constant energy transfer slices in precisely
the same way as the experimental data are sliced. This resolution function
calculation properly accounts for the functional form of the contributions
from individual instrument components. However, the resolution convolution is
quite computationally intensive, requiring approximately 12 CPU-hr per
parameter set (240$\times$400 image). The largest amount of time is spent by
the spin wave calculation itself, that is, within _SpinW_ , because the
resolution convolution requires the calculated structure factor at a large
number of wave vector and energy transfer $(\mathbf{Q},\hbar\omega)$
positions, each of which requires the computation and diagonalisation of a
Hamiltonian matrix and additional matrix multiplications to obtain the spin-
spin correlation function and hence neutron scattering intensity. In order to
speed this up we used another application, _Brille_ [33], to pre-compute these
structure factors in a dense grid within the first Brillouin zone, and then to
linearly interpolate within this grid for the required
$(\mathbf{Q},\hbar\omega)$ points needed by the resolution convolution. This
reduced the calculation time to $\approx$1 CPU-hr per parameter set.

We also explored an approximate and faster resolution convolution method,
motivated by the general desire to reduce the number of CPU hours needed when
using ML techniques. Instead of using expensive Monte Carlo integration, the
resolution function covariance matrix is pre-computed for each incident energy
spectrum in a grid in $(\mathbf{Q},\hbar\omega)$ space and then used to define
a 4D Gaussian profile with which the spin wave theory dispersion is convolved.
With an 80$\times$80 grid per 40$\times$40 pixel slice the calculation times
were reduced to $\approx$10 CPU-min per parameter set. However, the resulting
calculated spectra yielded larger classification errors as described below.
The coarseness of the grid resulted in an underestimate of the resolution
broadening, and consequently the calculation is “sharper” than the data and
the energy transfer of the top (bottom) of the spin wave dispersion bands were
calculated to be too low (high) in energy transfer. The problem can be
rectified by using finer grids but this then results in similar calculation
times to the “Tobyfit” method with _Brille_ interpolation. Furthermore, the
use of a pre-computed grid also results in aliasing effects for certain grid
sizes which the Monte Carlo method is not susceptible to. The important lesson
from this experience is that accurate account of the resolution function needs
to be taken for robust interpretation of INS, despite the concomitant cost in
computing resources.

We generated 3322 images for each of the Goodenough and dimer models, split
6000:644 between training and validation datasets, which required a total of
$\sim$7000 CPU-hr. The parameters for these datasets are randomly generated by
independently selecting from a uniform distribution within the limits noted in
Table 1 for each $J_{ij}$ and model. These limits bracket likely maximum
ranges of estimates of the parameters from the bandwidths of the spin wave
branches and values of the exchange constants determined in various other
manganites. This mimics the procedure that would be followed if the training
datasets were being generated for an unknown system.

The data can be downloaded from an online repository [34]. Figure 2 shows the
calculated spectra using the Monte Carlo resolution convolution method for
both the Goodenough and dimer models alongside the measured data.

##  4 A network for phase discrimination

###  4.1 Can a neural network learn to distinguish phases?

We first want to establish that we can train a neural network (NN) capable of
learning to distinguish between the Goodenough and dimer models based on the
INS spectrum. This is a classic inverse problem; while the forward mapping
from the spin Hamiltonian for each model to the INS spectrum is well known,
the inverse is not.

Analysing the spin wave spectrum pixel-by-pixel using a multi-layer perceptron
(MLP) type NN architecture would lead to a rapidly exploding number of
parameters. In addition, as an MLP treats the input image as a vector of data,
it can be very sensitive to small shifts in spectral feature positions which
may not be characteristic of the different exchange models. While these shifts
do not pose a problem for the training of the network which is done on
generated data, it may become an issue when confronted by experimental data,
because slight miscalibrations of the spectrometer may produce data which is
offset compared to the ideal simulations. Furthermore, due instrumental
resolution broadening, neighbouring pixels in a spectrum (image) are
correlated, and a MLP would have to learn these correlations explicitly.
Instead, we employ the popular convolutional neural network (CNN) architecture
[35, 36], whereby the image (spectrum) is first passed through a series of
filters, before feeding the values of these filters into a MLP architecture.
During training the weights of the filters and the MLP are updated. The
filters can be seen as learning a reduced representation of the spectrum
(including implicit correlations between neighbouring pixels), while the MLP
learns the function to map these features to the corresponding exchange
interaction model.

Figure 3: Network architectures. (a) The convolutional neural network (CNN)
architecture; at each (green) convolution layer the number of filters is
indicated below and the size of the convolution matrix above the blocks, and
at each fully connected (orange) layer the number of input dimensions is
indicated. After each block of convolution layers we perform a MaxPool (brown
block with pool size indicated) and batch normalisation (see text). The large
dense layer (344448 element orange block) is a vector obtained by flattening
the outputs of the final layer of convolution filters, this layer is then
connected to a 64 node layer and finally to the classification head. (b) The
same CNN as (a), but rather than flatten the final layer of filters a global
pooling of each of the 32 filters is used to generate 32 values which are
connected directly to the classification head; this architecture is used to
construct class activation maps. (c) The variational autoencoder architecture;
at each convolution (green) layer the number of filters is indicated, we also
perform $(2\times 2)$ striding (the convolution filter skips across the input
in steps of two, thereby halving the size of the output) at each convolutional
layer, the central latent (orange) layer consists of $20$ Gaussian
distributions as characterised by a mean and standard deviation.

There are a wide array of CNN architectures already trained and available for
image recognition problems, for example the alexnet [37] and resnext [38]
architectures which are both previous winners of the imagenet [39] data
classification challenge. In this work we build a relatively simple CNN
architecture, as the data that we are classifying and the number of potential
classes are quite small. As we will demonstrate later, choosing simpler
architectures can help with interpretation of the final CNN inference.

The architecture we built is represented schematically in Figure 3 (a). The
convolution layers are coupled with a pooling layer so that the image is
subsequently down-sized as it passes through the convolutions and a batch
normalisation is applied after pooling. In batch normalisation the outputs
from a layer are normalised to fall between 0 and 1 across a pre-defined batch
size of training samples. In MaxPool the size of a filter is reduced between
layers by pooling sections of the filter and taking the maximum value of each
pool. At the final convolutional layer the values of each filter
$\mathcal{F}^{j}_{xy}$ (where $x$ and $y$ are spatial coordinates and $j$ is
the filter index) are flattened to produce a $117\times 92\times
32=344448$-element vector 222A convolutional filter with size $n$, padding $p$
and stride $s$ changes an image with input dimension $d$ to $(d-n+2p)/s+1$; in
our classification CNNs $p=0$ and $s=1$ for all cases. That is, after the
first $5\times 5$ convolutional layer, the image (initially $240\times 400$)
size is $236\times 396$ and after the first $1\times 2$ MaxPool its size is
$234\times 197$. The number of filters we used at each layer is specified in
Figure 3. The final convolution layer has 32 filters and the image has been
downsized to $117\times 92$ pixels. which is fed into a fully-connected layer
(MLP) with $d=64$ output nodes. Classification is then performed by a weighted
summation of these 64 values $v_{i}$ based on trainable weights $w^{c}_{i}$
and a trainable bias $b_{c}$, which serve as input to a sigmoid function
$\phi$. That is, the score for class $c$ is

$\displaystyle Y^{c}=\phi(\sum_{i=1}^{i=d}w^{c}_{i}v_{i}+b_{c})$ (2)

with the class with the highest score selected by the network.

The simulated training data are split into training and validation sets as
described in the data generation section. To save memory the training dataset
is divided into batches of 32 images. In a single training iteration, the
network is fed these 32 images to calculate the loss function (we use the mean
binary cross-entropy loss of the classification of all images in the batch as
the loss function). The gradient of each adjustable weight in the network with
respect to the loss is then calculated using backpropagation and the weight is
adjusted in the direction of steepest descent. This iteration is repeated for
each batch of 32 images, until all batches have been passed through the
network once, which is termed an epoch. At the end of each epoch the loss from
the separate validation dataset is calculated. If this validation loss is much
larger than the loss calculated from the training set after adjusting the
weights, the network is overfitting the training data such that it cannot
generalise its inference to the similar data in the validation set. The
training is set to terminate if the validation accuracy has not improved for
20 epochs, up to a maximum of 500 epochs. The full details of the network
training and the code for the architecture are available in an online
repository [34]. The final trained CNN achieves $>96\%$ accuracy in predicting
the correct exchange model for data in the validation set for both Goodenough
and dimer datasets. The training curves are included as supporting
information.

###  4.2 How much can we trust the predictions?

When a network has been trained to separate data into two classes, it can be
difficult to know how reliable a prediction on a new piece of data is. Binary
classification (Equation 2) applies a sigmoid function, which intentionally
favours output values close to either 0 or 1, forcing separation between
classes; however, this means that an individual sample may be classified
strongly into a class despite a degree of ambiguity in the result. For
example, when we trained the previously described network on data simulated
with the Monte Carlo resolution convolution, and then feed it the experimental
data spectrum, it gives class scores corresponding to
$[Y^{\mathrm{dimer}},Y^{\mathrm{Goodenough}}]$ of $[0.012,0.988]$. However,
when the same network was trained using data simulated using the faster pre-
computed grid of covariance matrix resolution convolution method, it gives
$[0.999,0.001]$. That is, it not only wrongly classifies the experimental data
as being in the dimer class, it does so strongly. It is thus incorrect and
misleading to interpret this output vector of the sigmoid classification as a
measure of confidence for that classification. To address this concern we have
adapted the recently published deterministic uncertainty quantification (DUQ)
scheme to work with our classification network [40].

There are a range of methods available for uncertainty quantification in NN
predictions, for example Bayesian neural nets (BNNs) [41]. Pure BNNs are
intractable for exact inference, although a number of approximations have been
proposed, including the simple-to-implement MC Dropout approach [42]. In
practice most of these Bayesian approaches are outperformed by the Deep
Ensembles (DE) approach [43], in which multiple networks with the same
architecture are trained from different initial values of weights and biases
and dataset orderings, resulting in an distribution of answers. The drawback
of the DE approach is the requirement to train and run multiple networks and
the associated linear growth of computational cost. We instead choose the DUQ
approach, which can provide uncertainty estimates in a single forward pass of
the network, but shows performance in out-of-distribution data detection
(identifying cases where predictions would be an extrapolation from training
data, rather than an interpolation) on a par with DE approaches [40]. We note
that several new approaches show good potential for uncertainty quantification
in NNs, such as the stochastic differential equation approach [44] and radial
BNNs [45], but a full comparison of all approaches is beyond the scope of the
current study.

Figure 4: Schematic representation of the deterministic uncertainty
quantification (DUQ) method. The input initially passes through a series of
convolutions (orange block) to extract features. In standard logistic
regression the outputs from the convolutions are classified by summing the
weights connecting each filter $f_{i}$ to the class $\mathcal{C}$ of interest,
as in equation 2. DUQ clusters the examples based on distances $K_{c}$ in a
high dimensional space of the outputs from the convolutions from the centre
$\mathbf{e}_{c}$ of clusters of training examples, according to equation 3.

The principle of DUQ as applied here is as follows. Let
$M_{\Theta}:\mathbf{x}\in\mathbb{R}^{n}\rightarrow\mathbf{v}\in\mathbb{R}^{d}$
be the transformation learnt by the network which maps an input vector
$\mathbf{x}$ (the input image, size $240\times 400$ pixels) with dimension
$n=96000$ to an output features vector $\mathbf{v}=\\{v_{1},\ldots,v_{d}\\}$
with dimension $d=64$ (the number of output nodes connected to the
classification head), and $\Theta$ are the parameters (weights and filters) of
the network. In DUQ [40] (Figure 4), instead of the final dense layer (the
classification head) which maps $\mathbf{v}$ to the class scores
$\mathbf{Y}=\\{Y^{1},\ldots,Y^{c}\\}$, a different mapping is learnt which
transforms $\mathbf{v}$ into a vector space with dimension $m$ where the
samples of the same class $c$ are clustered together. This is embodied by a
_per-class_ learnable weight matrix $\mathbf{W}_{c}\in\mathbb{R}^{m\times d}$
with $m=64$ in our case. The correlation $K_{c}$ between this
$\mathbf{W}_{c}\mathbf{v}$ vector and the class centroid $\mathbf{e}_{c}$,
given by

$K_{c}(\mathbf{v},\mathbf{e}_{c})=\exp\left[-\frac{||\mathbf{W}_{c}\mathbf{v}-\mathbf{e}_{c}||^{2}_{2}}{2n\sigma^{2}}\right]\
,$ (3)

is then a measure of the uncertainty in the classification, with the hyper-
parameter $\sigma$ being a characteristic length scale, also learned during
the training. The network then assigns a given input sample $\mathbf{x}$ to
the class with the largest correlation $K_{c}$. Thus instead of classifying
the sample based on the highest score $Y^{c}$ from equation 2, it is
classified as belonging to the centroid closest to it, with a confidence
$K_{c}$ that corresponds to the Euclidean distance between the new point and
the centroid. In this way, data that are out of the training distribution are
far from any of the trained centres of gravity and the confidence for
classification is low.

The network architecture we now use to classify the experimental neutron
scattering data is the same one demonstrated in the previous section to
reliably predict the correct exchange model with the validation datasets, and
is depicted in Fig 3 (a). The only difference is that the vector at the 64
node layer is then used as the input to the DUQ analysis. This CNN was trained
for 100 epochs on data prepared as described in the Training Data section.

Networks were trained on the simulated datasets for both the Goodenough and
dimer models, generated by the resolution convolution methods described
earlier: one using the fast but approximate method using a pre-computed grid
of covariance matrices ($M^{\mathrm{pre}}_{\Theta}$), and the other using the
more expensive but accurate Monte Carlo integration
($M^{\mathrm{MC}}_{\Theta}$) method. We then use these networks to classify
the background-subtracted experimental spectrum. The DUQ networks yields a
correlation value between 0 and 1 to indicate the distance in hyperspace
between the weight vector associated with the input and the output classes.

The network trained on pre-computed resolution data
$M^{\mathrm{pre}}_{\Theta}$ gives an output corresponding to
$[K_{\mathrm{dimer}},K_{\mathrm{Goodenough}}]$ of $[0.73,0.49]$ , showing that
the network correlates this input more with the dimer model than with the
Goodenough model. This classification is incorrect; moreover, the correlation
values suggest the the network does not distinguish between the two models
with great certainty and the experimental data lie out of distribution of the
training data.

To demonstrate this point further we have calculated the correlation values
$K_{\mathrm{Goodenough}}$ that the network assigns to the Goodenough category
for every spectrum in the training set for $M^{\mathrm{pre}}_{\Theta}$. The
correlation values are presented in Figure 5, demonstrating clearly that on
data corresponding to the training distribution the network will
overwhelmingly predict values very close to 1.0 for the (ground truth)
Goodenough class and close to 0.0 for the dimer class. This further emphasises
that the result from $M^{\mathrm{pre}}_{\Theta}$ on the background subtracted
experimental data indicates that these data are not within the training
distribution.

Figure 5: The correlation values $K_{c}$ for the $c=\mathrm{Goodenough}$
category assigned by the $M^{\mathrm{pre}}_{\Theta}$ DUQ classifier. The plot
is a histogram of values of the Goodenough class output given by the neural
network on the full training datasets of spectra generated using the
Goodenough (orange) and dimer (green) models (3322 spectra each).

In contrast, when the network trained on Monte Carlo resolution-convolved data
$M^{\mathrm{MC}}_{\Theta}$ was tested with the background subtracted
experimental spectrum, the output corresponding to the classes
$[K_{\mathrm{dimer}},K_{\mathrm{Goodenough}}]$ was $[0.05,0.99]$. The
$M^{\mathrm{MC}}_{\Theta}$ network not only predicts the correct result, but
does so with a high degree of confidence, indicating that the training data
with the Monte Carlo resolution convolution method covers a distribution
encapsulating the experimental data with the background subtracted.

This result highlights two important points. First the inclusion of accurate
representations of instrument and experiment resolutions is critical if one
wishes to train a neural network on simulated data and later to apply it to
analyse experimental results. Second having a method to quantify the
confidence of the machine learning is extremely important. In a standard CNN
with logistic regression rather than DUQ the $M^{\mathrm{pre}}_{\Theta}$ would
have classified the spectrum as belonging to the dimer class, with no
indication that there was uncertainty in this classification or that the data
were out of distribution.

To further demonstrate the utility of the DUQ method for identifying out of
distribution data we have tested the $M^{\mathrm{MC}}_{\Theta}$ network on
some “Franken-data”. The Franken-data are generated by stitching together
spectra from several different materials measured with the same instrument
settings (and hence instrument resolution) as the PCSMO data. The resulting
image, constructed to have the same dimensions as the PCSMO data, looks like a
reasonable INS spectrum at first glance. The $M^{\mathrm{MC}}_{\Theta}$
network gives correlation values of $[0.35,0.85]$, where the closeness of the
two values is indicative of a sample outside of the training distribution,
indicating the utility of DUQ for identifying cases where the classification
should not be trusted.

###  4.3 Why does the network predict what it does?

We have established that a CNN is capable of performing the inverse problem of
inferring the magnetic exchange model from the simulated spectrum. However, a
highly pertinent question is: can the network tell us what influenced this
decision? If we can address this question we can gain confidence in the
network, and can also use it in a predictive fashion, so that in advance of
analysing data the CNN could be used to work out what are the most important
regions for attention in the data. Moreover, knowing what in the data
influenced the decision is what helps yields scientific insight.

To ask the network which are the important regions for making a distinction we
employ class activation maps (CAMs) [46, 47]. CAMs relate regions of an input
signal to the output classification. The mechanism for producing a CAM is
outlined in Figure 6. In this example, the final convolutional layer is
directly connected to the classification head (a dense fully-connected layer),
so the extent to which a filter $\mathcal{F}^{j}_{xy}$ is responsible for a
classification $Y^{c}$ is determined by the weight $w_{j}^{c}$ connecting the
$j^{\mathrm{th}}$ filter to the $c^{\mathrm{th}}$ class. The activation map is
then obtained by:

$A^{c}_{xy}=\frac{1}{n_{f}}\sum_{j=1}^{j=n_{f}}w_{j}^{c}\mathcal{F}^{j}_{xy}\
.$ (4)

That is, the image $A^{c}_{xy}$ is the average (denoted $\langle\ldots\rangle$
in Figure 6) of the $n_{f}$ filters weighted by their contributions to the
classification $c$.

Figure 6:  Class activation map: The input initially passes through a series
of convolutions (orange block) to extract features. The final layer of filters
($\mathcal{F}_{xy}^{j}$ summed to give $f^{j}$) are fed into the
classification head as the global average of each filter to do classification
(red block). To build an activation map, we take the cross product of the
filter weights with the activation linking the average of that filter to the
class of interest. For a single image we take the average of all of these
cross products to build a map of why that particular class was activated,
equation 5.

This approach however, cannot work with the classification architecture we
introduced in section 4.1 where there is an intervening fully connected layer
between the final convolutional layer and the classification head, because it
requires the weights $w_{j}^{c}$ implying that the $j^{\mathrm{th}}$ filter
directly impact the class $c$. Thus, we use the Grad-CAM [47] method where
instead of the weights $w_{j}^{c}$, the gradients $\frac{dY^{c}}{df^{j}}$ are
used where $f^{j}=\sum_{xy}\mathcal{F}^{j}_{xy}$ is the sum of all pixels
(global pool) of a filter. These gradients can be calculated through any
number of intermediate layers between the last convolutional layer and the
classification head using backpropagation. Thus, the activation map we used is
given by:

$A^{c}_{x,y}=\frac{1}{n_{f}}\sum_{j=1}^{j=n_{f}}\frac{dY^{c}}{df^{j}}\mathcal{F}^{j}_{xy}\
.$ (5)

Applying this Grad-CAM to the original architecture, depicted in Figure 3 (a),
yielded a class activation map which does not reveal any activation patterns
which correspond to physical intuition. Thus, whilst the extra complexity from
the intermediate densely connected layer between the final convolution layer
and the classification head allows the classification network to be flexible
and to produce an accurate determination between the Goodenough and dimer
models, it also obscures the relationships between the spatial inputs of the
model and the outputs.

We next build a simpler network architecture (Figure 3 (b)), in this case
rather than flattening the final convolutional layer we take the global
average of each of the filters and make an input vector of size 32 (one for
each filter) and connect this directly to the classification head of the
network. This means that the filters now directly contribute to the
classification and the more activated a filter is, the more it contributes.
However, when presented with the experimental spectrum this network classifies
with very little certainty, as the new network is less flexible and less
robust to experimental noise.

To overcome the sensitivity to experimental noise we build a variational
autoencoder (VAE) to clean the signal [48]. Previous work has shown the
application of autoencoders for cleaning diffuse (i.e. energy integrated)
neutron scattering signals [11]. We build a VAE with two convolutional layers
in the encoder, with strided convolutions in order to compress the image size
- the architecture is shown in Figure 3 (c). At the “bottle-neck” or “latent”
layer we have 20 units. The network is then trained on 3500 examples of
simulated dimer and Goodenough datasets (1750 of each) with Monte Carlo
resolution convolution, the training is run for 100 epochs to maximise the
evidence lower bound (ELBO) loss[48, lin2019balancing]. For a VAE the network
trains to reproduce the input, so no labels are required. ELBO balances a loss
on the reconstruction (how similar output is to input) with the divergence
between the inferred distribution in the latent space and the prior
distribution (normal distribution with mean 0 and variance 1) of the latent
space. Because the data pass through the reduced dimensionality of the latent
layer the information content is necessarily reduced and only the stronger
features of the initial signal are retained. In this way noise can be removed
from the original, the result of cleaning the signal with the autoencoder is
shown in Figure 7. The full neural network architecture and training routine
are available in the repository associated with this paper [34].

Figure 7:  Upper: the effect of cleaning the experimental spectrum with a
variational autoencoder, left is the experimental spectrum, right is the same
spectrum passed through the trained autoencoder. Lower: Class activation maps.
Left shows the class activation map for the autoencoder cleaned signal, right
we have taken the mean of the activations per slice of the spectrum and
plotted these to show which slices contributed most strongly to the
classification of Goodenough.

Figure 7 shows the experimental spectrum, with the background subtracted (top
left) and passed through a cleaning autoencoder (top right), above a class
activation map (bottom left), which highlights the regions of this spectrum
that were important for the classification as a Goodenough structure. The
neural network classifier now correctly infers that the input signal
corresponds to the Goodenough model, with the class scores $[0,1]$
corresponding to $[Y^{\mathrm{dimer}},Y^{\mathrm{Goodenough}}]$. Looking at
the activation maps, in particular the averages per slice of the input (bottom
right), the classification is most heavily influenced by the signal in the
slice corresponding to 35.6 to 40.4 meV in the spectrum collected at
$E_{i}=70$ meV. This is the region at the top of the lower spin wave band and
is precisely the same region as was used in the original paper to make the
distinction. The activation map also has higher input from this same region at
the top of the lower band in the spectra collected at $E_{i}=100$ and 140 meV,
however this is not as strong as the contribution from the $E_{i}=70$ meV
spectrum.

These results demonstrate the importance of balancing complexity versus
interpretability, when one is interested in understanding why a network
provides the results that it does. In the case of the highly accurate network,
which could take into account experimental noise we are able to obtain the
correct result, however the layers of abstraction within the network make it
difficult to explain these results in an interpretable way. When we construct
a simpler network we have to remove some of the experimental noise in order to
make it work properly, however the result can be interpreted in a way which
corresponds to human understanding. This is an example of the well-known
trade-off between interpretability and completeness of explanations [49].

We have implemented the CAM networks in an interactive jupyter notebook [34].
In the notebook both the networks described above can be loaded and tested on
the data before and after cleaning with the autoencoder. The notebook also
features a tool for zooming in on slices of energy transfer space to look at
the difference between dimer and Goodenough model spectra.

##  5 Future challenges

While we have addressed the classification problem in this work, another
particular area of interest which we have not covered is the use of NNs for
regression to obtain a set of exchange parameters for a given model and given
measured data. With conventional fitting it might take a researcher some
months (work in practice sometimes spread over years) to obtain a set of
exchange parameters for a given model and measured dataset depending on the
complexity of the data. This is due to the need first for a careful
examination of the data for diagnostic regions and then a exhaustive search of
the parameter space. Despite the required training time (perhaps a few weeks
on a cluster), NNs offer a way to obtain these exchange parameter values on
much shorter timescales.

The challenge in applying these functionalities (classification and
regression) more generally are twofold. First, for users to take full
advantage of the new methods, the workflows must be relatively easy to use and
this will involve chaining together several disparate software packages in
several languages. At the ISIS Neutron and Muon Source this will be addressed
by the PACE project which aims to bring together data analysis tools for INS
[50].

Second, as we noted in the introduction, one key feature of single-crystal INS
data from modern spectrometers is their large size. The Pr(Ca,Sr)2Mn2O7
dataset explored in this work is two orders of magnitude smaller than those
typical of full 4D wave vector - energy transfer space maps of excitations,
and yet we had to reduce it further to a 400$\times$240 pixel image for input
to the neural network (the full data consist of 32 million
$(\mathbf{Q},\hbar\omega)$ bins although many of these will have zero counts).
Thus INS data pose almost the opposite problem to that often tackled by
machine learning, in that other applications can rely on a large number of
datasets with each dataset being relatively small in size; for INS there are
typically few _experimental_ datasets, necessitating the use of synthetic data
for training, but each dataset is almost unmanageable in size.

New techniques will have to be developed to address this challenge, and it is
likely that machine learning methods will play a role. For example, although
the total number of data bins is large, in many cases a significant fraction
of these contain scattering that is not of relevance, for example, a slowly
varying function of wave vector and energy transfer that forms a “background”
underlying what is known _a priori_ to be sharply defined set of excitations
of interest, such as the spin waves in the present work. If a reliable way can
be found to identify, that is, segment this ”background”, then we can treat
only the signal from the non-”background” data. This signal might be amenable
to Fourier analysis or to treatment using sum rules, or may be transformed
into a compressed representation which might be more appropriate for input to
a neural network.

##  6 Conclusion

We have applied deep neural networks to analyse INS data. This is the first
time (to our knowledge) that NNs have been applied directly for INS data
analysis. We demonstrated that a CNN architecture could be trained to classify
spin Hamiltonians based on simulated data. We then demonstrated how
uncertainty quantification can be included in the classification using the DUQ
method. DUQ highlights that an accurate calculation of the experimental
resolution is required to train a network using synthetic data that then works
well on real data. DUQ was also able to detect when data from a different
material to the training data but measured with the same spectrometer settings
was passed to the network.

We built class activation maps to determine which parts of the input spectrum
are important for the network classification and found that while more complex
networks can classify correctly with noisy data, the results are not
physically interpretable. Simpler (fewer layer) networks, however, yield more
interpretable class activation maps at the cost of requiring cleaner data
which had been passed through an autoencoder. In this case, the maps showed
that the most diagnostic feature of the measured spectrum is the region at the
top of the lower excitation band, which accords with previous work [26].

We believe that this work demonstrates how trustable, interpretable ML can be
a powerful tool in the analysis of INS data. We have highlighted future
directions that can be explored to further harness the potential of ML in this
field. We hope that this work is an early step in using ML methods not just
for processing, but for better understanding INS experiments.

## Acknowledgements

Experiments at ISIS were supported by a beamtime allocation RB1310483 from the
Science and Technology Facilities Council. TGP thanks co-authors A.T.Boothroyd
and D.Prabhakaran of ref. [26] for permission to use the datasets from those
experiments. This work was partially supported by Wave 1 of The UKRI Strategic
Priorities Fund under the EPSRC Grant EP/T001569/1, particularly the “AI for
Science” theme within that grant and The Alan Turing Institute. The simulated
datasets were generated using computing resources provided by STFC Scientific
Computing Department’s SCARF cluster.

## Data Access Statement

All of the training data, trained neural networks and code for generating the
training data for this study are openly available at
https://zenodo.org/record/4270057#.X65xr1BpHIU.

A git repository containing the code used to build and train the neural
networks, as well as notebooks to recreate the DUQ and CAM experiments is
available at https://github.com/keeeto/interpretable-ml-neutron-spectroscopy

The software environment required to run the codes in the git repository is
available in the form of docker images from https://hub.docker.com/u/mducle
with instructions in the github repository. Additionally there are Conda
environment files provided in the repository.

## Author Contributions

KTB, MDL and TGP jointly conceived, planned and steered the project. KTB
built, trained and applied the neural networks; MDL produced the simulated
training datasets; TGP provided the experimental data. KTB, MDL and TGP wrote
the manuscript together in an iterative fashion to ensure the contents present
neural networks and related methods in a manner that is accessible and useful
to the condensed matter physics community. JT was involved in conception and
establishing of the project and facilitated the work in the paper.

## References

  * [1] Hey T, Tansley S, Tolle K et al. 2009 The fourth paradigm: data-intensive scientific discovery vol 1 (Microsoft Research Redmond, WA) ISBN 978-0-9825442-0-4 URL https://www.microsoft.com/en-us/research/publication/fourth-paradigm-data-intensive-scientific-discovery/
  * [2] Agrawal A and Choudhary A 2016 Apl Materials 4 053208
  * [3] Carrasquilla J and Melko R G 2017 Nature Physics 13 431–434
  * [4] Butler K T, Davies D W, Cartwright H, Isayev O and Walsh A 2018 Nature 559 547–555
  * [5] Bednik G 2019 Physical Review B 100 184414
  * [6] Morita K, Davies D W, Butler K T and Walsh A 2020 arXiv preprint arXiv:2005.05831
  * [7] Hey T, Butler K, Jackson S and Thiyagalingam J 2020 Philosophical Transactions of the Royal Society A 378 20190054
  * [8] Doucet E, Brown T, Chowdhury P, Lister C, Morse C, Bender P and Rogers A 2020 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 954 161201
  * [9] Islam F, Lin J Y Y, Archibald R, Abernathy D L, Al-Qasir I, Campbell A A, Stone M B and Granroth G E 2019 Rev. Sci. Instrum. 90 105109 (Preprint https://arxiv.org/abs/1906.09482) URL https://doi.org/10.1063/1.5116147
  * [10] Hui Y and Liu Y 2019 Science and Information Conference (Springer) pp 257–271 URL https://arxiv.org/abs/1710.05994
  * [11] Samarakoon A M, Barros K, Li Y W, Eisenbach M, Zhang Q, Ye F, Sharma V, Dun Z L, Zhou H, Grigera S A, Batista C D and Tennant D A 2020 Nature Communications 11(1) 892 URL https://doi.org/10.1038/s41467-020-14660-y
  * [12] Archibald R K, Doucet M, Johnston T, Young S R, Yang E and Heller W T 2020 Journal of Applied Crystallography 53
  * [13] Demerdash O N, Shrestha U, Petridis L, Mitchell J C, Smith J and Ramanathan A 2019 Frontiers in molecular biosciences 6 64
  * [14] Qian X and Yang R 2018 Physical Review B 98 224108
  * [15] Lovesey S 1984 Theory ofneutron scattering from condensed matter (Clarendon Press, Oxford)
  * [16] Squires G 1978 Introduction to the theory of thermal neutron scattering (Cambridge University Press, Cambridge)
  * [17] Chen T, Chen Y, Kreisel A, Lu X, Schneidewind A, Qiu Y, Park J, Perring T G, Stewart J R, Cao H et al. 2019 Nature materials 18 709
  * [18] Wang M, Zhang C, Lu X, Tan G, Luo H, Song Y, Wang M, Zhang X, Goremychkin E, Perring T et al. 2013 Nature communications 4 1–10
  * [19] Kieslich G, Skelton J M, Armstrong J, Wu Y, Wei F, Svane K L, Walsh A and Butler K T 2018 Chemistry of Materials 30 8782–8788
  * [20] Li X, Liu P F, Zhao E, Zhang Z, Guidi T, Le M D, Avdeev M, Ikeda K, Otomo T, Kofu M et al. 2020 Nature Communications 11 1–9
  * [21] Bewley R, Taylor J and Bennington S 2011 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 637 128 – 134 ISSN 0168-9002
  * [22] Abernathy D L, Stone M B, Loguillo M J, Lucas M S, Delaire O, Tang X, Lin J Y Y and Fultz B 2012 Review of Scientific Instruments 83 015114 URL https://doi.org/10.1063/1.3680104
  * [23] Kajimoto R, Nakamura M, Inamura Y, Mizuno F, Nakajima K, Ohira-Kawamura S, Yokoo T, Nakatani T, Maruyama R, Soyama K, Shibata K, Suzuya K, Sato S, Aizawa K, Arai M, Wakimoto S, Ishikado M, Shamoto S i, Fujita M, Hiraka H, Ohoyama K, Yamada K and Lee C H 2011 Journal of the Physical Society of Japan 80 SB025 URL https://doi.org/10.1143/JPSJS.80SB.SB025
  * [24] Ollivier J and Mutka H 2011 Journal of the Physical Society of Japan 80 SB025 ISSN 0031-9015
  * [25] Ewings R, Buts A, Le M, van Duijn J, Bustinduy I and Perring T 2016 Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 834 132 – 142 ISSN 0168-9002
  * [26] Johnstone G E, Perring T G, Sikora O, Prabhakaran D and Boothroyd A T 2012 Phys. Rev. Lett. 109(23) 237202 URL https://link.aps.org/doi/10.1103/PhysRevLett.109.237202
  * [27] Tokura Y 2006 Reports on Progress in Physics 69 797–851 URL https://doi.org/10.1088%2F0034-4885%2F69%2F3%2Fr06
  * [28] Goodenough J B 1955 Phys. Rev. 100(2) 564–573 URL https://link.aps.org/doi/10.1103/PhysRev.100.564
  * [29] Kanamori J 1959 Journal of Physics and Chemistry of Solids 10 87 – 98 ISSN 0022-3697 URL http://www.sciencedirect.com/science/article/pii/0022369759900617
  * [30] Daoud-Aladine A, Rodríguez-Carvajal J, Pinsard-Gaudart L, Fernández-Díaz M T and Revcolevschi A 2002 Phys. Rev. Lett. 89(9) 097205 URL https://link.aps.org/doi/10.1103/PhysRevLett.89.097205
  * [31] Toth S and Lake B 2015 Journal of Physics: Condensed Matter 27 166002
  * [32] Perring T G 1991 High energy magnetic excitations in hexagonal cobalt Ph.D. thesis University of Cambridge
  * [33] Tucker G S 2020 https://brille.github.io/stable/index.html
  * [34] Butler K T and Le M D 2020 https://github.com/keeeto/interpretable-ml-neutron-spectroscopy
  * [35] LeCun Y, Boser B, Denker J S, Henderson D, Howard R E, Hubbard W and Jackel L D 1989 Neural computation 1 541–551
  * [36] Fukushima K and Miyake S 1982 Competition and cooperation in neural nets (Springer) pp 267–285
  * [37] Krizhevsky A, Sutskever I and Hinton G E 2012 Advances in neural information processing systems pp 1097–1105
  * [38] Xie S, Girshick R, Dollár P, Tu Z and He K 2017 Proceedings of the IEEE conference on computer vision and pattern recognition pp 1492–1500
  * [39] Deng J, Dong W, Socher R, Li L J, Li K and Fei-Fei L 2009 2009 IEEE conference on computer vision and pattern recognition (Ieee) pp 248–255
  * [40] van Amersfoort J, Smith L, Teh Y W and Gal Y 2020 International Conference on Machine Learning URL https://arxiv.org/abs/2003.02037
  * [41] Neal R M 2012 Bayesian learning for neural networks vol 118 (Springer Science & Business Media)
  * [42] Gal Y and Ghahramani Z 2016 international conference on machine learning pp 1050–1059
  * [43] Lakshminarayanan B, Pritzel A and Blundell C 2017 Advances in neural information processing systems pp 6402–6413
  * [44] Kong L, Sun J and Zhang C 2020 arXiv preprint arXiv:2008.10546
  * [45] Farquhar S, Osborne M A and Gal Y 2020 stat 1050 7
  * [46] Zhou B, Khosla A, Lapedriza A, Oliva A and Torralba A 2016 Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) URL https://www.cv-foundation.org/openaccess/content_cvpr_2016/html/Zhou_Learning_Deep_Features_CVPR_2016_paper.html
  * [47] Selvaraju R R, Cogswell M, Das A, Vedantam R, Parikh D and Batra D 2017 Proceedings of the IEEE international conference on computer vision pp 618–626
  * [48] Kingma D P and Welling M 2013 arXiv preprint arXiv:1312.6114 URL https://arxiv.org/abs/1312.6114
  * [49] Gilpin L H, Bau D, Yuan B Z, Bajwa A, Specter M and Kagal L 2018 2018 IEEE 5th International Conference on data science and advanced analytics (DSAA) (IEEE) pp 80–89
  * [50] Buts A, Battam N, Ewings R A, Fair R L, Jackson A, Le M D, Marooney C, Perring T G, Saunders H and Tucker G S 2020 PACE: Proper Analysis of Coherent Excitations https://github.com/pace-neutrons


# Long-range adiabatic quantum state transfer through tight-binding chain as
quantum data bus

Bing Chen1, Wei Fan1, Yan Xu1,3, Zhao-yang Chen2, Xun-li Feng3 and C. H. Oh3
1School of Science, Shandong University of Science and Technology, Qingdao
266510, China  2 Department of Mechanical Engineering, University of
California, Berkeley, CA 94720, USA  3 Center for Quantum Technologies and
Physics Department, Faculty of Science, National University of Singapore, 2
Science Drive 3, Singapore 117542

###### Abstract

We introduce a scheme based on adiabatic passage which allows for long-range
quantum communication through tight-binding chain with _always-on_
interaction. By adiabatically varying the external gate voltage applied on the
system, the electron can be transported from the sender’s dot to the aim one.
We numerically solve the schrödinger equation for a system with given number
of quantum dots. It is shown that this scheme is a simple and efficient
protocol to coherently manipulate the population transfer under suitable gate
pulses. The dependence of the energy gap and the transfer time on system
parameters is analyzed and shown numerically. Our method provides a guidance
for future realization of adiabatic quantum state transfer in experiments.

###### pacs:

03.67.Hk, 03.65.-w, 73.23.Hk

##  I Introduction

Quantum state transfer (QST), as the name suggests, refers to the transfer of
an arbitrary quantum state from one qubit to another, which is a central task
in quantum information science. For the solid-state based quantum computing at
the large scale, it is very crucial to have a solid system serving as such
quantum data bus, which can provide us with a quantum channel for quantum
communication. During the last years many efforts have been made in different
fields to design a feasible proposal for perfect QST. One kind of proper QST
proposals is based on solid-state system with _always-on_ interaction Bose1 ;
Song ; Christandle1  . The communication is achieved by simply placing a
quantum state at one end of the chain and waiting for an optimized time to let
this state propagate to the other end with a high fidelity. The other kinds of
proposals have paid much attention to adiabatic passage for coherent QST in
time-evolving quantum systems, which is a powerful tool for manipulating a
quantum system from an initial state to a target state. This way of population
transfer has the important property of being robust against small variations
of the Hamiltonian and the transport time, which is crucial in experiment
since the system parameters are often hard to control. The typical scheme for
coherently spatial population transfer has been independently proposed for
neutral atoms in optical traps Eckert  and for electrons in quantum dot (QD)
systems CTAP  via a dark state of the system, which is termed coherent
tunneling via adiabatic passage (CTAP) following Ref. CTAP . In such a scheme,
the tunneling interaction between adjacent quantum units is dynamically tuned
by changing either the distance or the height of the neighboring potential
wells following a counterintuitive scheme which is a solid-state analog of the
well-known stimulated Raman adiabatic passage (STIRAP) protocol STIRAP  of
quantum optics. Since then, the CTAP technique has been proposed in a variety
of physical systems for transporting single atoms atom1 ; atom2 , spin states
spin , electrons electron1 ; electron2  and Bose-Einstein condensates BEC1 ;
BEC2 ; BEC3 . It has also been considered as a crucial element in the scale up
to large quantum processors LR1 ; LR2 .

Recently, Ref. chen1  presented a scheme to adiabatically transfer an electron
from the left end to the right end of a three dot chain using the ground state
of the system. This technique is a copy of the frequency chirping method CF1 ;
CF2 , which is used in quantum optics to transfer the population of a three-
level atom of the Lambda configuration. The scheme chen1  is presented as an
alternative to a well known transfer scheme (CTAP) CTAP . However, different
from CTAP process, the protocol in Ref. chen1  considers a three QD array with
_always-on_ interaction which can be manipulated by the external gate voltage
applied on the two external dots (sender and receiver). Through maintaining
the system in the ground state, it shows that it is a high-fidelity process
for a proper choice of system parameters and also robust against experimental
parameter variations. The obvious extension of this work is to consider the
passage through more than one intervening dot chen2 . In this paper we will
consider a quasi-one-dimensional chain of QDs, which is schematically
illustrated in Fig. 1. The central tight-binding chain serves as the quantum
channel and two external QDs are attached to the media chain. The sender
(Alice) and the receiver (Bob) control one external QD each and Alice
transmits information via a qubit using the chain to Bob by adiabatically
changing the gate voltages. Different from previously discussed schemes, we
will consider a fixed $N$-site coupled QDs media chain and QST can be realized
in required transfer distance by modulating the positions where QD A and B are
connected to the chain. In particular, the nearest-neighbor hopping amplitudes
are set to be uniform. We first theoretically elaborate the adiabatic QST in
this scheme. Taking a 50-dot structure as an example, we show that the
electron can be robustly transported from Alice to Bob through the media
chain, by slowly varying the gate voltages.

The paper is organized as follows. In Sec. II the model is setup and we
describe the adiabatic transfer of an electron between QDs. In Sec. III we
show numerical results that substantiate the analytical results. The last
section is the summary and discussion of the paper.

##  II Model Setup

Figure 1: (a) Schematic illustrations of adiabatic QST in multi-dot array. The
system is controlled by external gates voltage $\mu_{i}(t)$ $(i=A,$ $B)$. By
adiabatically varying the gates voltage, one can achieve long-range QST from
the QD $A$ to QD $B$ of the chain. (b) Gate voltages as a function of time (in
units of $\tau$). $\mu_{A}(t)$ is the solid line and $\mu_{B}(t)$ is the dash
line.

Consider a quasi-one-dimensional chain of QDs, realized by the empty or singly
occupied states of a positional eigenstate, see Fig. 1. The whole quantum
system consists of two sites ($A$ and $B$) and a simple tight-binding $N$-site
chain. The sender (Alice) and the receiver (Bob) can only control the external
gates voltage $\mu_{\alpha}(t)$ $(\alpha=A,$ $B)$. The total Hamiltonian
$\mathcal{H}(t)=\mathcal{H}_{M}+\mathcal{H}_{I}+\mathcal{H}_{C},$contains
three parts, the medium Hamiltonian

$\mathcal{H}_{M}=-J\sum_{j=1}^{N-1}\left|j\right\rangle\left\langle
j+1\right|+\text{h.c.}$ (1)

describing the tight-binding chain with uniform nearest neighbor hopping
integral $-J$ $(J>0),$ the coupling Hamiltonian

$\mathcal{H}_{I}=-J_{0}\left(\left|A\right\rangle\left\langle
l\right|+\left|B\right\rangle\left\langle
l^{\prime}\right|\right)+\text{h.c.},$ (2)

describing the connections between QDs $A$, $B$, and the chain with hopping
integral $-J_{0}$ $(J_{0}>0)$, and the operating Hamiltonian

$\mathcal{H}_{C}=\mu_{A}(t)\left|A\right\rangle\left\langle
A\right|+\mu_{B}(t)\left|B\right\rangle\left\langle B\right|$ (3)

describing the adiabatic manipulation of the Hamiltonian parameters. In
$\mathcal{H}(t)$, $\left|j\right\rangle$ represents the Wannier state
localized in the $j$-th quantum site for $j=A,$ $1,$ $2,...,$ $N,$ $B$. In
$\mathcal{H}_{I}$, $l$ and $l^{\prime}$ denote the sites of medium connecting
to the QDs $A$ and $B$. The distance between $A$ and $B$ is
$D=l^{\prime}+1-l$. In this proposal, we just consider one connection way:
$l^{\prime}=N+1-l$, that is the quantum state is transferred between site $A$
and its _mirror-conjugate_ site $B$. In term $\mathcal{H}_{C}$, $\mu_{A}(t)$
and $\mu_{B}(t)$ are site energies (externally controlled), which are
modulated by a Gaussian pulses (shown in Fig. 1(b))

$\displaystyle\mu_{A}(t)$ $\displaystyle=$
$\displaystyle-\mu_{0}\exp\left[-\frac{1}{2}\alpha^{2}t^{2}\right],$
$\displaystyle\mu_{B}(t)$ $\displaystyle=$
$\displaystyle-\mu_{0}\exp\left[-\frac{1}{2}\alpha^{2}\left(t-\tau\right)^{2}\right],$
(4)

where $\mu_{0}$ is the peak voltage of the pulse; $\tau$ and $\alpha$ the
total adiabatic evolution time and standard deviation of the control pulse. To
realize high fidelity transfer in this scheme, the peak voltage $\mu_{0}$ must
be much larger than hopping integral, i.e. $\mu_{0}\gg J,J_{0}$. The reason is
that small peak values improve adiabaticity, but lead to a low fidelity
because the final instantaneous eigenstate is not the desired one. According
to Ref. chen2 , the transfer is optimized when we choose $\alpha=8/\tau$.
Throughout this paper, all energies ($J_{0}$ and $\mu_{0}$) are scaled in
units of $J$, and evolution time $\tau$ is in units of $1/J$.

In this proposal, we focus our study on the ground state
$\left|\psi_{g}(t)\right\rangle$ of Hamiltonian $\mathcal{H}(t)$ to induce
population transfer from state $\left|A\right\rangle$ to
$\left|B\right\rangle$. For single electron transfer, a state in the single
particle Hilbert space is assumed as
$\left|\psi_{k}\right\rangle=f_{A}^{k}\left|A\right\rangle+\sum_{j=1}^{N}f_{j}^{k}\left|j\right\rangle+f_{B}^{k}\left|B\right\rangle$,
where $k$ denotes the momentum. Duo to the translational symmetry of the
present system, the instantaneous Hamiltonian’s eigen equation for
$f_{j}^{k},$ $j\in\left[1,N\right]$ is easily shown to be

$\displaystyle-J\left[f_{j-1}^{k}+f_{j+1}^{k}\right]$ $\displaystyle=$
$\displaystyle\left[\varepsilon_{k}-V_{A}\left(\varepsilon_{k}\right)\delta_{j,l}-V_{B}\left(\varepsilon_{k}\right)\delta_{j,N+1-l}\right]f_{j}^{k},$
$\displaystyle J_{0}f_{A}^{k}$ $\displaystyle=$ $\displaystyle
V_{A}\left(\varepsilon_{k}\right)f_{l}^{k},$ (5) $\displaystyle
J_{0}f_{B}^{k}$ $\displaystyle=$ $\displaystyle
V_{B}\left(\varepsilon_{k}\right)f_{N+1-l}^{k},$

here $\varepsilon_{k}$ is the eigenenergy and the term

$V_{i}\left(\varepsilon_{k}\right)=\frac{J_{0}^{2}}{\varepsilon_{k}+\mu_{i}},i\in\left[A,B\right]$
(6)

on the right hand side is contributed by the interactions between the sites
$A$, $B$, and medium chain which is dependent on the eigenenergy
$\varepsilon_{k}$. The $\delta$-type potential forms a confining barrier to
the transportation of single electron in the chain and forms a bounded state
of single electron, similar to those proposed in Ref. Sun . In this work, we
focus our attention on the bound state, which is the ground state of the total
system to realize the long-range QST.

Starting from $t=0$, we have $\mu_{A}(0)=-\mu_{0}$ and $\mu_{B}(0)\approx 0$.
The solution to Eq. (5) is

$\left|\psi_{g}(0)\right\rangle=\mathcal{N}^{-1/2}\left[\sum_{j=1}^{N}e^{-\kappa\left|j-l\right|}\left|j\right\rangle+\frac{\lambda_{+}}{J_{0}}\left|A\right\rangle+\frac{\lambda_{-}}{J_{0}}e^{-\kappa(N+1-2l)}\left|B\right\rangle\right],$
(7)

where $\lambda_{\pm}=\left(\sqrt{\mu_{0}^{2}+4J_{0}^{2}}\pm\mu_{0}\right)/2$
and $\kappa=\ln\left(\lambda_{+}/J\right)$;
$\mathcal{N}=\sum_{j=1}^{N}e^{-2\kappa\left|j-l\right|}+\left(\mu_{0}^{2}+2J_{0}^{2}\right)/J_{0}^{2}$
is the normalization factor. By choosing a sufficiently large value of
$\mu_{0}$, the ground state $\left|\psi_{g}(0)\right\rangle$ can be reduced to
$\left|\psi_{g}(0)\right\rangle\approx\left|A\right\rangle$.

With the same results, in the time limit $t=\tau$, the parameter $\mu_{A}(t)$
goes to zero and $\mu_{B}(t)$ goes to $-\mu_{0}$. Duo to the reflection
symmetry (relabeling sites from right to left) of the system, for
$j=A,1,2,...,N,B,$

$\langle
j\left|\psi_{g}(\tau)\right\rangle=\langle\bar{j}\left|\psi_{g}(0)\right\rangle.$
(8)

We have used $\bar{j}=N+1-j$ to indicate the _mirror-conjugate_ site of $j$.
This leads to the final ground state
$\left|\psi_{g}(\tau)\right\rangle\approx\left|B\right\rangle$. To illustrate
with an example, the probability of $\left|B\right\rangle$ in
$\left|\psi_{g}(\tau)\right\rangle$ can achieve 99.7% when the parameters are
set to be $J_{0}/J=1$ and $\mu_{0}/J=20$. Preparing the system in state
$\left|\Psi\left(t=0\right)\right\rangle=\left|A\right\rangle\,$and
adiabatially changing $\mu_{A}(t)$ and $\mu_{B}(t)$, one can see that the
system will end up in $\left|B\right\rangle$,

$\left|\Psi\left(t=0\right)\right\rangle=\left|A\right\rangle\rightarrow\left|\Psi\left(t=\tau\right)\right\rangle=\left|B\right\rangle,$
(9)

hence we can see that a high fidelity transfer may still be possible, even
with imperfect controls.

In the absence of hopping term between the two external QDs $A$, $B$ and the
medium chain ($J_{0}=0$), the Hamiltonian $\mathcal{H}(t)$ can be diagonalized
as

$\mathcal{H}(t)=\sum_{k}-2J\cos k\left|k\right\rangle\left\langle
k\right|+\mu_{A}(t)\left|A\right\rangle\left\langle
A\right|+\mu_{B}(t)\left|B\right\rangle\left\langle B\right|$

with
$\left|k\right\rangle=\sqrt{2/\left(N+1\right)}\sin(kj)\left|j\right\rangle$,
where $k=n\pi/(N+1),n=1,2,...,N$. One can see that at $t=\tau/2$, the
eigenstates $\left|A\right\rangle$ and $\left|B\right\rangle$ are degenerate
which leads to a breakdown of adiabaticity. The presence of the hopping term
$J_{0}$ will open up energy gap at the crossing. To evaluate instantaneous
eigenvalues of the Hamiltonian is generally only possible numerically. In Fig.
2(a) we present the results showing the eigenenergy gap between the
instantaneous first-excited state and ground state undergoing evolution due to
modulation of the gate voltages according to pulse Eq. (4). The eigenvalues
shown in this figure exhibit pronounced avoided crossing and approach nonzero
minimum
$\Delta=\varepsilon_{1}\left(\tau/2\right)-\varepsilon_{g}\left(\tau/2\right)$.
This minimum energy gap plays a significant role in the transfer, because the
total evolution time should be larger enough compared to $\Delta$. In this
scheme, the energy gap $\Delta$ depends both on the transfer distance $D$ and
the coupling constance $J_{0}$. To study the relationship between the total
evolution time and system parameters is one of the important contributions of
this paper.

As an example, we show in Fig. 2 the effect two factors has on the energy gap
$\Delta$ for a system with $N=48$ QDs and coupling strength $J=1.0$. The
pulse’s parameters we choose are $\alpha=8/\tau$, and $\mu_{0}=20$. In Fig.
2(b), we plot the energy gap $\Delta$ as a function of transfer distance $D$
for $J_{0}=0.5J$, and $J_{0}=0.7J$. The horizontal line $\delta\approx
3J\pi^{2}/N^{2}$ indicates the minimum gap of medium chain $\mathcal{H}_{M}$.
One can see that the logarithmic scales chosen suggest that as transfer
distance $D$ increasing there is an exponential disappearance of the gap
$\Delta$. The smaller hopping constant $J_{0}$, the slower the decay of
$\Delta$. The other thing is that $\Delta$ is also determined by the coupling
strength $J_{0}$. Fig. 2(c) shows the numerically computed behavior of
$\Delta$ as a function of $J_{0}$ for $D=8$, $16$, and $24$. As $J_{0}$
increases, the gap $\Delta$ increase for short-distance transfer ($D<N/3$) and
decreases for long-distance transfer ($D>N/3$). The results shown in Fig. 2(c)
also suggest that decreasing coupling strength $J_{0}$ can obtain relatively
large $\Delta$ for long-range transfer. But it does not mean the weaker the
coupling $J_{0}$, the better the result of QST will be. The reason is that the
energy gap is not the sufficient and necessary condition for adiabatic
process. In this proposal, the negative effects of distance on the gap can be
partially compensated by $J_{0}$.

Figure 2: (a) The instantaneous eigenenergy (in units of $J$) of the lowest
two states $\psi_{g}$ and $\psi_{1}$ through the gate pulse shown in Fig.
1(b). The gap is minimum at $t=\tau/2$,
$\Delta=\varepsilon_{1}\left(\tau/2\right)-\varepsilon_{g}\left(\tau/2\right)$.
(b) Plot of gap $\Delta$ (in units of $J$) as a function of transfer distance
$D$ for $J_{0}=$0.5$J$, and 0.7$J$. The energy gap $\Delta$ is exponentially
decreasing as transfer distance $D$ increase for a given $J_{0}$. (c) The gap
$\Delta$ (in units of $J$) as a function of coupling constant $J_{0}$ for the
transfer distance $D=$8, 16, and 24. The horizontal line $\delta$ is the
minimum gap of the medium chain $\mathcal{H}_{M}$.

##  III Numerical Simulations

In this section let us firstly review the transfer process of this protocol.
At $t=0$ we initialize the device so that the electron occupies site-$1$,
i.e., the total initial state is $\left|A\right\rangle$. Provided we transform
the gate pulses adiabatically, then the adiabatic theorem states that the
system will stay in the same eigenstate. Therefore, the quantum state starting
in $\left|A\right\rangle$ will end up in $\left|B\right\rangle$.

The analysis above is based on the assumption that the adiabaticity is
satisfied. The adiabaticity parameter defined for this scheme is

$\mathcal{A}\left(t\right)=\frac{\left|\left\langle\psi_{g}(t)\right|\partial\mathcal{H}/\partial
t\left|\psi_{1}(t\right\rangle\right|}{\left|\varepsilon_{g}(t)-\varepsilon_{1}(t)\right|^{2}},$
(10)

where $\left|\psi_{g}(t\right\rangle$ ($\left|\psi_{1}(t\right\rangle$) is the
instantaneous ground state (first-excited state) of the Hamiltionian
$\mathcal{H}(t)$ and $\varepsilon_{g}(t)$ ($\varepsilon_{1}(t)$) is the
corresponding instantaneous eigenvalue of the state
$\left|\psi_{g}(t\right\rangle$ ($\left|\psi_{1}(t\right\rangle$). For
adiabatic evolution of the system we require $\mathcal{A}\left(t\right)\ll 1$
for all time, which greatly suppresses the quantum transition from the ground
state $\left|\psi_{g}(t)\right\rangle$ to the first-excited state
$\left|\psi_{1}(t)\right\rangle$. Fig. 3(a) and (b) show the numerical result
of energy difference $\varepsilon_{1}(t)-\varepsilon_{g}(t)$ and
$\mathcal{A}\left(t\right)\tau$ as a function of pulse time. Note that the
appearance time of maxima of $\mathcal{A}\left(t\right)$ is not at the middle
of the pulse sequence, but the energy difference
$\varepsilon_{1}(t)-\varepsilon_{g}(t)$ at this time is slightly larger than
the minimum gap $\Delta$. So we can use minimum gap $\Delta$ to estimate the
minimum pulse time required for high-fidelity transfer. Fig. 3(c) shows the
maximum adiabaticity $\max\\{\mathcal{A}\left(t\right)\tau\\}$ through the
protocol as a function of $J_{0}$ for $D=10$, $16$, and $20$. One sees that
there is an optimal value of $J_{0}$ which ensures the shortest time for
realizing perfect QST and the optimal value decreases as transfer distance $D$
increases. To sum up, the adiabatic regime necessary to obtain transport with
high fidelity can be concluded to the condition $\tau\gg
1/\Delta^{\ast}\left(D\right)$, where $\Delta^{\ast}\left(D\right)$ is the
minimum gap of the system when the coupling strength $J_{0}$ takes the optimal
value corresponding to transfer distance $D$.

Figure 3: (a) The energy difference $E_{1}(t)-E_{g}(t)$ (in units of $J$) and
(b) adiabaticity $\mathcal{A}\left(t\right)\tau$ (in units of $10^{3}/J$) as a
function of time in the time interval $t\in[0.3\tau,0.7\tau]$. Note that the
energy gap approximately equals $\Delta$ in a wide time range and the maxima
of adiabaticity is not at the middle of the pulse sequence. (c) Maximum
adiabaticity $\max\\{\mathcal{A}\left(t\right)\tau\\}$ through the protocol as
a function of $J_{0}$ for three different transfer distances.
$\max\\{\mathcal{A}\left(t\right)\tau\\}$ varies with $J_{0}$, and reaches
minimum value when the $J_{0}$ takes specific value.

The consequent time evolution of the state is given by the Schrödinger
equation (assuming $\hbar=1$)

$i\frac{d}{dt}\left|\Psi\left(t\right)\right\rangle=\mathcal{H}(t)\left|\Psi\left(t\right)\right\rangle.$
(11)

The time evolution creates a coherent superposition:

$\left|\Psi\left(t\right)\right\rangle=c_{A}(t)\left|A\right\rangle+\sum_{j=1}^{N}c_{j}(t)\left|j\right\rangle+c_{B}(t)\left|B\right\rangle,$
(12)

where $c_{j}(t)$ denotes the time-dependent probability amplitude for the
electron to be in $j$-th QD. We define the probability of finding the electron
on the medium chain as
$\left|c_{M}(t)\right|^{2}=\sum_{j=1}^{N}\left|c_{j}(t)\right|^{2}$ that obeys
the normalization
condition$\,\left|c_{A}(t)\right|^{2}+\left|c_{M}(t)\right|^{2}+\left|c_{B}(t)\right|^{2}=1$.
At time $\tau$ the fidelity of initial state transferring to the dot-$B$ is
defined as $F=\left|\langle
B\left|\Psi(\tau)\right\rangle\right|^{2}=\left|c_{B}(\tau)\right|^{2}$.

Figure 4: (Color online) QST from QD $A$ to $B$ attached to a tight-binding
chain of length with $N=48$ for two different transfer distance: $D=11$
(a)-(d) and $D=21$ (e)-(h). (a) The transfer fidelity as a function of
coupling strength $J_{0}$ for transfer distance $D=11$. The calculation was
done for various values of $J_{0}$ at fixed total evolution time $\tau=480/J$.
Note that there is an optimal value of $J_{0}$ to achieve high fidelity QST.
Population transfer as a function of time obtained for three coupling strength
$J_{0}$:(b) $J_{0}=0.8J$; (c) $J_{0}=0.91J$; (d) $J_{0}=1.0J$. Initially the
population is on QD $A$ (dashed red line) and finally mainly on QD $B$ (solid
blue line). The population on the media chain is shown as a dotted black line.
(e) The same as in (a), but for transfer distance $D=21$ and $\tau=1200/J$.
(f)-(h) The population behavior under the influence of $J_{0}$. The parameters
are like that in (b)-(d).

In order to proceed, we used standard numerical methods to integrate the
Schrödinger equation for probability amplitudes. Because the scheme relies on
maintaining adiabatic conditions, we examine the effect of system parameters
on the target state population. In Fig. 4, we consider the system with $N=48$
and show QST from QD $A$ to QD $B$ for two different transfer distance: $D=11$
and $D=21$. Firstly, we examine the effect of coupling strength $J_{0}$ on
transfer fidelity. It is seen that there is an optimal value of $J_{0}$ to
achieve high-fidelity transfer.

To illustrate the process of QST for $D=11$, we exhibit in Fig. 4(b)-(d) the
exact evolution of the probabilities of finding electron in QD $A$ (red dashed
line), $B$ (blue solid line), and media chain (black dotted line) as a
function of time for three different values of coupling strength $J_{0}$ but
the same remaining parameters ($\mu_{0}=20J$, $\tau=480/J$ and
$\alpha=8/\tau$). We get good results for the transfer if we choose
$J_{0}=0.89J$ as shown in Fig. 4(b). The populations on the QD $A$ and QD $B$
are exchanged in the expected adiabatic manner. If we choose parameters
deviated from the optimum value, $J_{0}=0.8J$ and $1.0J$, we find the results
in Fig. 4(c) and Fig. 4(d). We can see that a slight deviations from the
values will breaks adiabaticity and lead to major deteriorations of the
quality of transfer. When the transfer distance becomes large, we should
enlarge evolution time $\tau$ to enhance the adiabaticity. In Fig. 4(e),
transfer fidelity as a function of $J_{0}$ for transfer distance $D=21$ and
evolution time $\tau=1200/J$. The time evolution of the probabilities the same
as Fig. 4(b)-(d) for $D=21$ and three different $J_{0}$ are illustrated in
Fig. 4(f)-(g). We can still see from Fig. 2 that the optimum value of $J_{0}$
to achieve high-fidelity transfer decreases as the transfer distance increase
which is consistent with the results shown before.

Figure 5: (Color online) For a given chain length $N=48$ and a chosen fidelity
$F=0.995$, (a) the optimum value of coupling strength $J_{0}$ (in units of
$J$) and (b) the minimum time $\tau$ (in units of $1/J$) required for state
transfer varies as a function of transfer distance $D$. The other system
parameters we choose is $J=1$, $\mu_{0}=20$ and $\alpha=8/\tau$. It shows that
optimum value of $J_{0}$ decreases as $D$ increasing and time $\tau$ scales
linearly with transfer distance $D$.

In order to provide the most economical choice of parameters for reaching high
transfer efficiency, we perform numerical analysis, as shown in Fig. 5.
Specifically, we depict the optimum coupling $J_{0}$ (shown in Fig. 5(a)) and
the corresponding minimum state transfer time $\tau$ (shown in Fig. 5(b)) for
a given chain length $N=48$ and a given tolerable transfer fidelity
$F=99.5\%$, the minimum time varies as a function of transfer distance $D$.
One can see that the time $\tau$ required for high-fidelity transfer scales
linearly with transfer distance $D$.

##  IV Summary and Discussion

An efficient QST scheme should not only admit a state transfer of any quantum
state in a fixed period of time of the state evolution with high fidelity, but
also the transfer time can not grow fast as communication distance increases.
In this paper, we have introduced a long-range transport mechanism for quantum
information around a quasi-one-dimensional QDs network, based on adiabatic
passage. This scheme is realized by modulation of gate voltages applied on the
two external QDs which is connected to the tight-binding chain. Under suitable
system parameters, the electron can be transported from the sender QD to the
receiver one with high efficiency, carrying along with it the quantum
information encoded in its spin. Different from the CTAPn Scheme CTAPn , our
method is to induce population transfer through tight-binding chain by
maintaining the system in its ground state and this is more operable in
experiments. We have studied the adiabatic QST through the system by
theoretical analysis and numerical simulations of the ground state evolution
of tight-binding model. The result demonstrates that it is an efficient high-
fidelity process (99.5%) for an economical choice of system parameters.
Increasing the transfer distance, we found that the efficiency of QST is
inversely proportional to the distance of the two QDs.

In a real system, decoherence is the main obstacle to the experimental
implementation of quantum information Ivanov ; Kam . There are two sources of
quantum decoherence in QDs, one is due to charge dephasing brought by lead-QD
coupling and the other is due to the hyperfine interaction. For the former,
the coherence time of quantum dot is $\sim$1 ns, which plays a role in this
case. On the other hand, the maximum time in $N=48$ QD system needed for the
appearance of the better fidelity is roughly proportional to $10^{4}/J$. As a
simple estimate of the effects of decoherence, we compare this time with the
dephasing time, which leads to a limit of coupling strength of $J$ of $\sim$10
THz. The probability of realization of this idea in experiment can be
maximized by more precise manipulation technology and by cooling the system.
Furthermore, the development of cold atom physics provides us with an
alternative realization of our systems in experiment, because decoherence in
cold-atom system is much less destructive.

## Acknowledgments

We acknowledge the support of the NSF of China (Grant No.10847150 and
No.11105086), the National Research Foundation and Ministry of Education,
Singapore (Grant No. WBS: R-710-000-008-271), the Shandong Provincial Natural
Science Foundation (Grant No. ZR2009AM026 and BS2011DX029), and the basic
scientific research project of Qingdao (Grant No.11-2-4-4-(6)-jch). Y. X. also
thanks the Basic Scientific Research Business Expenses of the Central
University and Open Project of Key Laboratory for Magnetism and Magnetic
Materials of the Ministry of Education, Lanzhou University (Grant No.
LZUMMM2011001) for financial support.

## References

  * (1) S. Bose, Phys. Rev. lett. 91, 207901 (2003).
  * (2) Z. Song and C. P. Sun, Low Temperature Physics 31, 686 (2005).
  * (3) M. Christandl, N. Datta, A. Ekert and A.J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
  * (4) K. Eckert, M. Lewenstein, R. Corbalán, G. Birkl, W. Ertmer, and J. Mompart, Phys. Rev. A 70, 023606 (2004).
  * (5) A. D. Greentree, J. H. Cole, A. R. Hamilton, and L. C. L. Hollenberg, Phys. Rev. B 70, 235317 (2004).
  * (6) N. V. Vitanov, T. Halfmann, B. W. Shore, and K. Bergmann, Annu. Rev. Phys. Chem. 52, 763 (2001).
  * (7) K. Eckert, J. Mompart, R. Corbalan, M. Lewenstein, and G. Birkl, Opt. Commun. 264, 264 (2006).
  * (8) T. Opatrný, K. K. Das, Phys. Rev. A 79, 012113 (2009).
  * (9) T. Ohshima, A. Ekert, D. K. L. Oi, D. Kaslizowski, L. C. Kwek, e-print arXiv:quant-ph/0702019.
  * (10) P. Zhang, Q. K. Xue, X. G. Zhao, and X. C. Xie, Phys. Rev. A 69, 042307 (2004).
  * (11) J. Fabian and U. Hohenester, Phys. Rev. B 72, 201304(R) (2005).
  * (12) E. M. Graefe, H. J. Korsch, and D. Witthaut, Phys. Rev. A 73, 013617 (2006).
  * (13) M. Rab, J. H. Cole, N. G. Parker, A. D. Greentree, L. C. L. Hollenberg, and A. M. Martin, Phys. Rev. A 77, 061602(R) (2008).
  * (14) V. O. Nesterenko, A. N. Nikonov, F. F. de Souza Cruz, and E. L. Lapolli, Laser Phys. 19, 616 (2009).
  * (15) L. C. L. Hollenberg, A. D. Greentree, A. G. Fowler, and C. J. Wellard, Phys. Rev. B 74, 045311 (2006).
  * (16) A. D. Greentree, S. J. Devitt, and L. C. L. Hollenberg, Phys. Rev. A 73, 032319 (2006).
  * (17) B. Chen, W. Fan, and Y. Xu, Phys. Rev. A 83, 014301 (2011).
  * (18) B. Chen, W. Fan, and Y. Xu, Sci China Ser G-Phys Mech Astron, in press.
  * (19) J. Cheng and J.-Y. Zhou, Phys. Rev. A 64, 065402 (2001).
  * (20) D. Goswami, Phys. Rep. 374, 385 (2003).
  * (21) D. Z. Xu, H. Lan, T. Shi, H. Dong, and C. P. Sun, Sci China Ser G-Phys Mech Astron, 53(7): 1234-1238 (2010).
  * (22) L. C. L. Hollenberg, A. D. Greentree, A. G. Fowler, and C. J. Wellard, Phys. Rev. B 74, 045311 (2006).
  * (23) P. A. Ivanov, N. V. Vitanov, and K. Bergmann, Phys. Rev. A 70, 063409 (2004).
  * (24) I. Kamleitner, J. Cresser, and J. Twamley, Phys. Rev. A 77, 032331 (2008).
  * (25) Z. Song, P. Zhang, T. Shi and C.-P. Sun, Phys. Rev. B 71, 205314 (2005).