# Application of artificial neural networks for rigid lattice kinetic Monte
Carlo studies of Cu surface diffusion

Jyri Kimari<EMAIL_ADDRESS>Ville Jansson  Simon Vigonski  Ekaterina
Baibuz  Roberto Domingos  Vahur Zadin  Flyura Djurabekova  Helsinki Institute
of Physics and Department of Physics, P.O. Box 43 (Pietari Kalmin katu 2),
FI-00014 University of Helsinki, Finland  Institute of Technology, University
of Tartu, Nooruse 1, 50411 Tartu, Estonia  Instituto Politécnico de Nova
Friburgo – Universidade do Estado do Rio de Janeiro, Rua Sao Francisco Xavier,
524, 20550–900, Rio de Janeiro, RJ, Brazil

###### Abstract

Kinetic Monte Carlo (KMC) is a powerful method for simulation of diffusion
processes in various systems. The accuracy of the method, however, relies on
the extent of details used for the parameterization of the model. Migration
barriers are often used to describe diffusion on atomic scale, but the full
set of these barriers may become easily unmanageable in materials with
increased chemical complexity or a large number of defects. This work is a
feasibility study for applying a machine learning approach for Cu surface
diffusion. We train an artificial neural network on a subset of the large set
of $2^{26}$ barriers needed to correctly describe the surface diffusion in Cu.
Our KMC simulations using the obtained barrier predictor show sufficient
accuracy in modelling processes on the low-index surfaces and display the
correct thermodynamical stability of these surfaces.

###### keywords:

Copper , Kinetic Monte Carlo , Artificial neural networks , Machine learning ,
Surface diffusion , Migration barriers

††journal: Computational Materials Science

##  1 Introduction

Diffusion in crystalline material is an important phenomenon in many
situations. For example, computational studies of irradiation damage may have
a hard time finding agreement with experiments without accounting for bulk
defect diffusion and annihilation after the initial collision cascade [1].
Surface diffusion, on the other hand, is hypothesised to play a role in e.g.
the events preceding vacuum arc breakdowns in devices with high electric field
gradients [2]. Vacuum arcs hinder the operation of many such devices,
including particle accelerators, free-electron lasers and fusion reactors. The
immediate motivation for studying the Cu surface specifically is the Compact
Linear Collider (CLIC) [3], which has been proposed to be built in CERN to
succeed the Large Hadron Collider.

Diffusion is difficult to study with molecular dynamics (MD) because it is a
much slower process compared to the MD timestep, which has to be small enough
to capture the atomic vibrations. Kinetic Monte Carlo (KMC) is a method
suitable for longer time scale studies, such as simulations of the diffusion
in crystalline solids. In this method, diffusion is approximated as a series
of migration events (jumps) between potential energy minima. The jumps, while
actually determined by Newton’s equations of motion with atoms vibrating most
of the time near potential energy minima, only occasionally crossing the
energy barriers between them, can in aggregate be regarded as stochastic
events that occur at rates determined by the height of those barriers.

Such migration energy barriers must be known for each event that is to be
considered in the KMC simulation. Depending on the desired specificity to
distinguish between different events, the number of barriers that must be
known may be too high to be calculated in a feasible time. Machine learning
was proposed earlier as an alternative approach—the problem of too many
barriers may be solved by only calculating accurately a subset of the barriers
and obtaining the rest from a computationally inexpensive regression model.
Djurabekova et al. used artificial neural networks (ANN) to predict migration
barriers in Fe-Cu bulk [4, 5]. Castin et al. [6, 7, 8, 9, 10, 11, 12, 13],
Pascuet et al. [14] and Messina et al. [15] have since applied them to study
bulk diffusion in various Fe-based alloys. Machine learning for surface
diffusion barriers have been used at least by Sastry et al. [16], who applied
genetic programming for vacancy-assisted 100 surface migration barriers in Cu-
Co alloy, and Verma et al. [17], who built a cluster expansion model for Ag,
Al, Cu, Ni, Pd and Pt (100) surface migration.

In this paper, we will study the capabilities of ANNs to predict migration
barriers on Cu surfaces. As an outcome of this study, we present a
configuration of trained neural networks for parameterizing a rigid lattice
KMC model of the (arbitrarily oriented) Cu surface. The networks and
instructions to utilise them are published in the parallel Data in Brief
publication [18]. On the one hand, this work is an extension of the rigid
lattice surface diffusion model developed earlier in our group [19], towards a
more detailed description of the atomic environments and thus, hopefully, more
accurate dynamics in the simulations. For the numerous challenges in the
migration barrier calculations on the surface, we deploy solutions developed
by Baibuz et al. [20]. On the other hand, this work applies the methodology of
refs. [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15] to a new system: the face-
centered cubic (fcc) crystal surface. We will discuss various obstacles of the
methods and how to overcome them.

The KMC method, barrier calculations, and ANN models that we have used are
described in section 2. The results are presented in section 3 and discussed
in section 4. Conclusions about the applicability of ANNs for the considered
problems are drawn in section 5.

##  2 Methods

###  2.1 Kinetic Monte Carlo

We used the residence-time KMC algorithm with the rigid lattice approximation,
as implemented in the Kimocs code [19].

KMC models the time evolution of a system by choosing events to be carried out
one after another. The probability of an event to be chosen is proportional to
its rate $\Gamma$, which in our model is given by the Arrhenius equation:

$\Gamma=\nu\exp\left(\frac{-E_{\mathrm{m}}}{k_{\mathrm{B}}T}\right),$ (1)

where $\nu$ is the attempt frequency, $E_{\mathrm{m}}$ is the migration energy
barrier, $k_{\mathrm{B}}$ is the Boltzmann constant and $T$ is temperature. In
our current model, $E_{\mathrm{m}}$ is strongly dependent on the local atomic
environment of each jump, while $\nu$ is taken to be the same for all events.
The $\nu$ parameter is connected to the vibration frequency of atoms, but it
can also been seen as a scaling factor for the simulation time.

###  2.2 Barrier calculations

The barrier data set was calculated using the nudged elastic band (NEB) [21,
22] method implemented in the LAMMPS molecular dynamics program [23]. The NEB
method finds the minimum energy path (MEP) between an initial and a final
configuration. The barrier is the difference between the highest energy point
along this path $E_{\mathrm{max}}$ (saddle point, since it is the maximum
along the MEP and the minimum along other directions) and the initial state
energy $E_{\mathrm{i}}$:

$E_{\mathrm{m}}=E_{\mathrm{max}}-E_{\mathrm{i}}$ (2)

This is also illustrated in fig. 1.

Figure 1: A schematic illustration of a 1-dimensional potential energy surface
from a NEB calculation. The barrier $E_{\mathrm{m}}$ and the reverse barrier
$E_{\mathrm{m,reverse}}$ can be obtained from the same calculation.

The potential energy function we chose was based on the molecular
dynamics/Monte Carlo corrected effective medium (MD/MC-CEM) theory, developed
by Stave et al. [24]. The potential is reported to describe the properties of
the Cu surfaces well [25].

Figure 3: The local atomic environment (LAE) used in the process descriptor.
_(Top left panel)_ The octahedral 28-atom cluster containing the migrating
atom and its 1nn and 2nn sites. 1nn sites are coloured light yellow and the
2nn sites are light grey. _(Right panel)_ The indexing of the neighbour sites
from 0 to 25 within the cluster. The initial position of the migrating atom is
marked blue and the final position red. _(Bottom left panel)_ Process “1 0 0 1
1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 0 1 1 1 0 0”: the blue atom jumps to the
position circled in red. The LAE has been embedded in the 100 surface. The
figures were generated with Ovito [26].

To connect the barriers to the corresponding processes during KMC simulations,
a process descriptor is required. This same descriptor can also be used when
training the ANN regressor. We describe the migration processes by their local
atomic environments (LAE) before the jumps take place. We include in the LAE
the first and the second nearest neighbours (1nn and 2nn) of the initial and
the final position, as was done in the earlier parameterizations of Kimocs
[19]. In systems with the fcc structure, this definition of LAE covers 26
lattice sites in total (see fig. 3). We did not expand the LAE beyond the 2nn
sites to keep the input space at a manageable size. The impact of the size of
the local atomic environment on the accuracy of the surface migration barriers
is subject to future study.

In ref. [19], the LAE was described using a 4-dimensional vector
$(a,\,b,\,c,\,d)$, where $a$ and $c$ are the number of 1nn atoms of the
initial and the final position of the jump and $b$ and $d$ are the number of
2nn atoms of the initial and the final positions of the jump. We will
hereafter refer to this approach as the “4D description”, as it was named in
ref. [20]. KMC simulations that employ the 4D description have been found to
produce good agreements with both MD and experimental results [19, 27, 20, 28,
29, 30]. The 4D description does not include the information on the exact
locations of the 26 LAE atoms, but only reflects the stability of the initial
and the final positions by counting the number of neighbours in each neighbour
shell. This descriptor is also not necessarily extensible to systems with
multiple atomic species present. To save the location information, in our new
descriptor the occupation state of each site is encoded with either 0 (vacant)
or 1 (occupied). The descriptor is a 26-dimensional binary vector, referred to
as the 26D description. The total number of different LAEs that can be
distinguished with this descriptor is $2^{26}\approx 67$ million. This choice
of LAE encoding does not restrict the fcc surface orientation that can be
described, and it is equally well applicable for bulk diffusion processes or
systems with more than one atomic species, by denoting different elements with
1, 2, etc. This type of encoding is similar to what was used in refs. [6, 7,
8, 9, 10, 11, 14, 13, 15] to describe processes in bulk systems with the body-
centered cubic (bcc) structure. A similar descriptor was also used by Trushin
et al. [31] for a _self-learning_ 2D KMC model fcc surfaces, and by Latz et
al. [32] for the 3D expansion of that model. A self-learning KMC model, unlike
a machine learning parameterized one, does not build a regression function for
migration barriers, but simply calculates them on-the-fly and saves each
result in a searchable catalogue for later use. This kind of model needs a
descriptor for labelling the catalogue, whereas a machine learning model uses
the descriptor as a way to encode input for the regression function.

The 26D descriptor is not compact in the sense that some LAEs are mirror
images of each other, and thus have the same migration energy barrier:
multiple descriptors have the same expected output. Up to four different LAEs
may belong into these “families” of equivalent processes. This is problematic
from the machine learning point of view: even though the redundant symmetric
cases may be included in the training data for the ANN, training will never be
perfect and the network output will be different for each case. If the ANN
predicts different barriers for mirrored, but otherwise physically equivalent,
LAEs, diffusion may work differently e.g. on different 100 facets like (100)
and (010), or be biased towards some directions. This problem was solved by
removing the redundant processes from the training set systematically: only
one process represents each family in the training set. When calling the
network to produce the barrier of a given jump, the LAE is first transformed
to correspond to the “representative” that is shown (or would be shown, in the
case of extrapolation) to the network in the training.

In the remainder of this paper we will refer to the LAE configuration around
the initial and final jump positions up to the second nearest neighbour shells
as the LAE cluster. Even though, as it is said, the 26D LAE cluster contains
the information about the neighbours up to the second nearest neighbour shell,
for barrier calculations by NEB, we embed such a cluster in a larger lattice,
where all the sites are occupied by the same atoms and affect the calculations
in the systematic manner—calculating the migration barriers within an isolated
cluster of a few dozen atoms would not be realistic. As a first approach we
embed the LAE cluster in Cu bulk for the NEB calculations, similarly as we did
in refs. [19, 20]. This is valid especially in studies where only a limited
number of vacancies is present in the LAE. However, when calculating surface
migration barriers, the LAE is on average half-empty. Previously we discovered
[20] that in some cases embedding such an LAE, essentially a large vacancy
cluster (void), behaves differently inside bulk compared to the surface,
causing very strong forces and unphysically high barriers in the range of tens
of eV. We thus returned to the original recipe, calculating all the barriers
with the LAE clusters embedded in a surface (see fig. 3).

Embedding the cluster in a surface is more complicated than in bulk. It is not
immediately obvious which surface should be used, and how each LAE should be
oriented to be best fit into the surface. To answer these questions, we
developed an automated procedure, where

  1. 1.

Three different surfaces were considered: 100, 110 and 111. Within the local
environment limited to 2nn sites, also the higher index surfaces will resemble
one of these lowest-index cases, and thus the resulting model can be generally
applied on any surface orientation.

  2. 2.

For each of the three surface orientations, trial configurations were
generated by embedding the LAE cluster in every possible orientation (with the
constraint that the lattice points of the LAE have to match the lattice points
of the surface). The LAE was embedded deep enough that at most one layer of
atoms was above the surface of the surrounding lattice.

  3. 3.

The centre-of-mass of the LAE cluster was calculated in each of the trial
configurations. The configuration with the lowest centre-of-mass with respect
to the surface was considered to be the most stable one, and was selected to
be used in the NEB calculation to find the barrier in this LAE.

We used a simulation cell of approximately $45\times 45\times 25\,\AA^{3}$ in
size. Depending on the specific crystal orientation, the exact dimensions of
the cell varied. Expanding the cell beyond this size did not affect the
barrier value significantly. The boundary conditions were periodic in the
horizontal dimensions and two fixed atom layers at the bottom.

Eleven replicas were used in the NEB calculations. In addition to the force
given by the MD/MC-CEM potential, a tethering spring force was applied on each
atom in the same way as in [20]. The tethering force keeps the atoms close to
their initial lattice sites, but does not directly contribute to the potential
energy of the system. Previous KMC simulations with barrier sets for Cu and
Au, calculated using the tethering method, have also shown good agreements
with both MD and experimental results [20, 28, 29, 30].

In our case, if any atom slips to unintended lattice sites during the NEB
relaxation, the process will no longer correspond to its 26D descriptor and
thus the obtained descriptor-barrier pair cannot be used in the KMC
simulations. This is a challenge specific to the rigid lattice KMC. The
problem is more severe when performing NEB calculations on the surface, where
there are other adatoms around the jumping atom, which may also have very few
neighbours and thus are not strongly bound to the their lattice site. These
loosely bound adatoms tend to easily “follow” the jumping atom, or otherwise
leave their original lattice sites. The tethering force constant, as
implemented in LAMMPS, was set to 2.0 eV/Å2. For details about the tethering
force approach, see [20].

The first set of barriers was calculated for adatom migration processes on
flat 100, 110 and 111 surfaces, with the migrating atom surrounded by
different configurations of atoms in the same atomic layer. The set formed by
these LAEs will be referred to as _flat surface processes_ for the remainder
of this paper. The ANNs were first tested by training and predicting barriers
within this flat surface set, comprising 3168 processes.

The set was then expanded by calculating the barriers for all processes where
every atom within the LAE cluster and the migrating atom had at least three
1nn atoms around it when embedded on a surface. In total, 11 652 085 processes
were found in this category—approximately 17 % of the entire LAE space.
Attempting to calculate the barriers for the rest of the configurations, where
some atoms have less than three neighbours, is likely to produce artificial
results—even with the tethering force applied. In KMC simulations with only
1nn jumps permitted, it is still beneficial to include these unstable events
in the simulation as they can serve as intermediate steps in e.g. 2nn or 3nn
jumps that would be possible in a real system. An ANN regression model, that
has learnt the general tendencies of the LAE dependent migration barriers, is
likely to give more realistic values for the barriers of these events than NEB
calculations that are less than reliable far from potential energy minima.
Thus, ANNs are a good way to expand the NEB-calculated set of barriers to
allow for more realistic kinetics.

###  2.3 Artificial neural networks

ANNs are a class of machine learning methods that can be used for
classification and function regression. A regressor is trained with known
input-output pairs to learn the underlying function and also to predict the
output for previously unseen input. In this work, ANNs were fitted to predict
values of the migration barrier function
$E_{\mathrm{m}}=E_{\mathrm{m}}(\mathrm{LAE})$.

Two different ANN models were studied: multilayer perceptrons (MLP) and radial
basis function (RBF) networks. The MLP implementation was taken from the Fast
Artificial Neural Network (FANN) library [33], and the RBF networks are from
Python’s SciPy package [34].

####  2.3.1 Multilayer perceptrons

MLPs consist of one input layer, one output layer, and one or more hidden
layers between them. The layers hold nodes that are connected to each other.
In a fully connected feed-forward network without shortcut connections, all
nodes of each layer are connected to all nodes of the next layer (see fig. 4
for a schematic illustration). This was one of the network structures used in
this work. The other type of MLP structure used here was a cascade network
where hidden nodes are added one by one during training so that each new node
is connected to each of the old nodes. This way the network will have $N$
hidden nodes in $N$ hidden layers with shortcut connections between each
layer.

Figure 4: A schematic illustration of a fully connected feed-forward MLP with
7 input nodes, 10 hidden nodes and 3 output nodes.

The input is mapped to the output by passing it through the hidden nodes: each
node calculates the weighted sum of its inputs and passes this value to the
next nodes through an activation function. The weights are chosen in an
iterative training process to produce minimal error in the training set. The
training algorithm used in this work was the FANN library implementation of
the improved resilient back-propagation (iRPROP) algorithm proposed by Igel
and Hüsken [35].

The purpose of the activation function is to allow for non-linear regression.
For the MLPs with static layout, we used sigmoid (logistic) activation
functions in the hidden nodes. In the cascade networks the Cascade2 [36]
algorithm sets the hidden node activation functions automatically. In the
output node we found linear activations function to perform best in all MLP
networks. With the default sigmoid activation function (scaled to accommodate
for the full range of energy barrier values) in the output node, the training
would often stagnate, with the network weights drifting towards extreme values
before sufficient accuracy was achieved. Cascade networks would completely
diverge for this problem with a sigmoid output activation function.

To avoid using negative barrier values given by the linear activation, we take
the barrier to be $E_{\mathrm{m}}=\max(0,\,E_{\mathrm{ANN}})$ during the
simulations.

Another post-processing step was to forbid processes that would lead to
complete detachment of atoms from the substrate. Examples of such processes
are not included in the training set, so the ANN-predicted barriers for them
are greatly underestimated.

The number of input and output nodes is determined by dimensionality of the
problem at hand. In the case of learning the 1D barrier values that correspond
to the 26D input, the input layer will have 26 nodes, and the output layer
will have one node. The number of hidden layers and nodes is another matter of
optimisation. In this study, a single hidden layer with 35 nodes was found to
be optimal for the static (non-cascade) network. For the cascade networks,
30–70 nodes were needed for reaching a sufficient accuracy. Overfitting was
not observed: the root mean square error in the validation set was equal to
the error in the training set.

####  2.3.2 Radial basis function networks

RBF networks are somewhat similar to MLPs. The difference is that the
activation functions are radially symmetric functions that depend on the
distance $r$ of the input vector of the _prototype vector_ of each node. The
Gaussian function is an example of a radial basis function. In this work, we
used multiquadric functions

$f(r)=\sqrt{\left(\frac{r}{\varepsilon}\right)^{2}+1}$ (3)

as basis functions; $\varepsilon$ is a width parameter that was automatically
adjusted by the library function.

RBF networks usually have only one hidden layer, with one basis function
corresponding to each hidden node. The nodes have their own prototype vectors,
and there is one weight vector that has a dimensionality equal to the number
of hidden nodes. The weights are set by a matrix inversion that is faster than
the iterative training procedure used with MLPs, but very memory-intensive for
large training sets. For training sets beyond the ~3000 barrier flat surface
set, this could not be done with the available resources—only MLP networks
were trained with the full 11.7 million barrier data set.

The prototype vectors for the RBF networks can be chosen from among the input
vectors in the training set and the number of these vectors could in principle
be smaller than the total number of data points. In the SciPy library this
functionality is not implemented—instead, there is always one prototype vector
for each data point. As a result, the accuracy of the obtained RBF network
will always be perfect for the training set. For this reason, in this work,
only 50 % of the 3000 barrier data set was used to train the RBF networks in
order to get an estimate on how well they can actually learn the migration
energy function. To make a fair comparison to MLPs, these were also trained
using only 50 % of the flat surface data.

####  2.3.3 Error reduction techniques

When the barrier set was expanded from the initial ~3000 barriers to the 11.7
million barriers, the MLP accuracy decreased significantly, with the cascade
networks performing the best. Two additional techniques were used to reduce
the error in barrier prediction. Firstly, the barriers calculated on each
surface—100, 110, or 111—were treated as separate sets and different ANN
predictors were trained to predict barriers on each of these sets. The issue
of knowing which predictor to call to get the barrier for an arbitrary LAE
encountered during KMC was solved by introducing an ANN classifier. This
network uses the same LAE input encoding, but instead of energy output, it
outputs the _class_ , or the surface that the input LAE corresponds to. We
used 1-of-C encoding for the surfaces, meaning that the classifier has three
output nodes with each node corresponding to one of the surface classes. A
very good success rate was achieved for the classifier (see table 1). During
the KMC simulation, every encountered LAE is first passed to the classifier
and an appropriate regressor based on the classifier output is used to obtain
the migration barrier. In the case that more than one classifier output
signals non-zero value (the classifier is uncertain of the correct surface
orientation), the barrier will be calculated as a weighted average of multiple
regressors.

Table 1: Confusion matrix for the surface classifier ANN. The numbers on the diagonal represent correctly classified LAEs. The total success rate was 99.39 %. | Classified as:
---|---
Actual: |  100 |  110 |  111
100 | 1 115 537 | 6640 | 17 104
110 | 3310 | 1 470 607 | 21 242
111 | 4580 | 18 582 | 8 994 483

Using this technique, the RMS prediction error decreased by approx. 10 %
compared to using only a single predictor for the entire data set.

Secondly, groups of ANN regressors that were trained on different, random
subsets of the training data were combined into regressor ensembles. The
predicted barrier for each LAE was taken to be the average of the predictions
given by the individual networks. The individual regressors were originally
the products of 5-fold cross-validation, but combining them into ensembles
turned out to make better predictions than any of the components. The RMS
prediction error decreased by 13–19 % compared to using only a single
predictor for each surface. Combining the ensembles and the classifier reduced
the RMS error by a total of 21.8 %: from 0.110 eV to 0.086 eV.

####  2.3.4 Setup of the kinetic Monte Carlo simulations

In addition to calculating the RMS error estimate, the accuracy of the barrier
predictions was assessed by implementing the ANN barrier function into the
Kimocs code and running KMC simulations with the predicted barriers. Three
different simulation scenarios were considered:

  1. 1.

Flattening of a 12 monolayer (576 atoms) cuboid nanotip on the three lowest-
index fcc surfaces, 100, 110, and 111. The system dimensions are the same as
the test cases in ref. [19].

  2. 2.

Relaxation of nanoparticles with different initial shapes. This will show how
well the model captures the correct thermodynamics—every shape should relax
approximately to the shape of the Wulff construction. Each shape was relaxed
in three different sizes: ~900 atoms, ~1400 atoms, and ~2100 atoms. The Wulff
constructions were created with the Atomic Simulation Environment [37], using
the surface energies reported by the authors of the Cu interatomic potential
[25].

  3. 3.

Stability of ¡110¿ nanowires. Vigonski et al. found crossing Au nanowires to
break through surface diffusion in KMC simulations near the junction points of
the wires [28]. Cu has the same lattice structure as Au, and thus it is
interesting whether crossing Cu nanowires are found to behave similarly. Cu
nanowire networks have applications in the large-scale fabrication of
transparent and flexible electronics [38]. The ¡110¿ orientation was chosen
since wires oriented this way can be constructed to have the thermodynamically
favourable 111 surfaces on all sides. Thick (1.1 nm in radius) ¡110¿ Cu wires
were found to be indefinitely stable in the original 4D parameterization of
Kimocs [19], so the radius was reduced down to ~0.5 nm to observe
fragmentation. Both single and crossing wires were simulated.

##  3 Results

###  3.1 Regression accuracy

Fig. 5 shows the distribution of barriers in the initial set of ~3000 flat
surface processes.

Figure 5: Distribution of migration barriers of the 3168 flat surface
processes.

Fig. 6 shows the accuracy of the static MLP, cascade MLP and RBF networks in
this flat surface set. 50 % of the set was used in training, while the
correlation is plotted in the full flat surface set.

Figure 6: The static MLP _(left)_ , the cascade MLP _(centre)_ , and the RBF
networks _(right)_ accuracy for flat surface migration barrier prediction. All
of the networks were trained using 50 % of the data set. Colouring is
according to the point density on a logarithmic scale. The RMS errors were
0.036 eV (static MLP), 0.033 eV (cascade MLP) and 0.024 eV (RBF).

The migration energy distribution in the full ~11.7 million barrier set is
shown in fig. 7. The large quantity of 0 eV barriers are for processes that
are spontaneous either in one direction (e.g. jumping towards a much higher
number of 1nn atoms), or, less frequently, both directions. The latter
processes include some events on the 111 surface that happen via hexagonal
close-packed (hcp) sites. The full barrier set can be found as a Data in Brief
entry [18].

Figure 7: Distribution of migration barriers in the full ~11.7 million barrier
data set.

Comparisons between the 11.7 million barrier full set and the 4D-parameterized
barrier sets are shown in figs. 8 and 9. In these figures the error bars
indicate the highest and the lowest values of the barriers, which were
calculated for different permutations within the same 4D description. The dark
dot in the middle shows the value of the mean barrier calculated for all
permutations. The mean value is calculated as a simple arithmetic average.
Following the notation used in refs. [20, 39], we will refer to two of such
sets as _Cu set 1_ and _Cu set 2_. Correspondence between physical local
atomic environments and the 4D descriptors is not one-to-one, and thus the 4D
barrier sets have been constructed by selecting one representative case from
the family of environments that correspond to each 4D descriptor. _Cu set 1_
was constructed by randomly selecting one of the permutations of nearest
neighbors, while _Cu set 2_ was selected by determining the lowest energy
configurations of the initial and final sites of a jump. The tethering force
approach was applied to the barrier calculations in _Cu set 2_ , with the same
tethering force constant 2.0 eV/Å2 as in this work. The general trend of
barriers in our 26D set is similar to the 4D sets, especially to _Cu set 2_.
It can be seen that the groups of configurations that correspond to the same
4D descriptor can have deviations of a few eV in energy barrier between
themselves.

Figure 8: Comparison of barriers of the full barrier set to the
4D-parameterized _Cu set 1_ of ref. [20, 39] that was calculated without the
tethering force. Correspondence between the 26D and the 4D descriptions is
many-to-one; red bars show the minimum and the maximum values, which are
calculated for all possible permutations of the same 4D description. The blue
dots are the mean values of the barriers in that range. Barriers over 1 eV are
generally lower in the 26D set than in the 4D set. Figure 9: Comparison of
barriers of the full barrier set to the 4D-parameterized _Cu set 2_ of ref.
[20, 39] that was calculated with the additional tethering force. Red bars and
the blue dots are the same as in fig. 8. Most of the mean values correlate
well with the values of _Cu set 2_ , but there is a region where the 26D set
estimates barriers to be lower than the corresponding 4D barriers.

The combination MLP regressor accuracy in the full ~11.7 million barrier set
is shown in fig. 10. As explained in section 2.3.3, the final barrier
predictor consists of three ensembles (one for each low-index surface) of five
regressors each and a classifier to combine them. Training specialist networks
to predict migration barriers on different physical surfaces improved total
regression accuracy by 10 %, compared to training a single regressor to the
entire data set. We propose that this improvement can be explained by the
networks implicitly learning the effect of different surface relaxation and
surface stress present on the different surfaces. These surface effects affect
the energy of all surface atoms in the system, and thus modify the surface
migration energy barriers.

Figure 10: Accuracy of the combination MLP in the full ~11.7 million barrier
data set. The plotted result is from 16 individual networks: an ensemble of
five for each surface type (100, 110 and 111) and a classifier to combine the
ensembles. The RMS error is 0.086 eV.

Even though KMC does not have an explicit total potential energy parameter,
the information about relative potential energy change in a jump can be
inferred from the forward and the reverse barrier for that jump: $\Delta
E=E_{\mathrm{fin}}-E_{\mathrm{ini}}=E_{\mathrm{m}}-E_{\mathrm{m,reverse}}$
(see figure 1). For a KMC model to produce thermodynamically correct behaviour
tending towards potential energy minima, these $\Delta E$ values must be
sufficiently accurate. We examined how accurately our machine learning model
reproduces this thermodynamical information; fig. 11 shows the comparisons
between the predicted $\Delta E$ and the values given by NEB. The overall
correlation of the $\Delta E$ values is good, even though the model was fitted
to the migration energy barriers only, and not the energy differences. We
note, however, that in the region of barriers with near-zero values, the
prediction of the $\Delta E$ values is less reliable. However, these barriers
describe mainly unstable configurations that can be encountered due to
limitations of the rigid lattice approximation. Since the barriers are low,
they are not expected to affect the overall thermodynamic behaviour of the
system.

Figure 11: The correlation between true and predicted differences between the
final and the initial energy values ($E_{\mathrm{fin}}-E_{\mathrm{ini}}=\Delta
E$) of every jump. Roughly 6.7 % of the correlation points fall in the top
left and the bottom right quadrants, representing cases where the
energetically favourable jumps are predicted to be unfavourable and vice
versa.

###  3.2 KMC simulations

We performed a number of KMC simulations to verify the applicability of the
developed ANN for simulations of surface diffusion processes. We chose three
typical scenarios of surface evolution driven by the surface energy
minimisation principle. These are _(i)_ the flattening of a Cu nanotip on
surfaces with different crystallographic orientation; _(ii)_ equilibrium shape
of a Cu nanoparticle reached after relaxation of a particles with the
different initial shape; _(iii)_ stability of Cu nanowires. All these
processes can contribute to evolution of a rough surface with various surface
features, which can be observed on surfaces subject to high electric fields.
Since in the current study we aim to validate the proposed model, we focus on
the processes under equilibrium condition (no electric field effects are yet
taken into account). It is important to show that the model is stable and
predicts physically reasonable behaviour of surfaces.

####  3.2.1 Nanotip flattening

A 12 monolayer nanotip flattening on 100, 110, and 111 surfaces was simulated
with KMC using the trained combination MLP as a barrier predictor. The attempt
frequency was fitted to the 110 case at different temperatures ranging from
850 K to 1200 K. For each temperature we performed 20 statistically different
simulations. In some simulations at high temperatures, the surface started to
self-roughen. We excluded these cases from the statistical averaging, since it
was not clear how to account for the evolution of the height of the tip. See
section 3.2.2 for more detail. Hence the final number of cases used for
statistical averaging were as follows: 18 (1000 K), 20 (1050 K), 13 (1100 K),
16 (1150 K), 6 (1200 K), and full 20 at all lower temperatures. The attempt
frequency value $\nu=2.81\times 10^{14}$ s-1 produced the best fit for the
flattening time $t_{\mathrm{f}}$, comparing to the MD results ranging from 850
K to 1200 K, reported in ref. [19] (see fig. 12). As the flattening times span
over many orders of magnitude, the $\nu$ parameter was fitted to minimise

$\mathrm{loss}=\left(\log t_{\mathrm{f,ANN}}-\log
t_{\mathrm{f,MD}}\right)^{2}$

The comparison shows the full parameterization to follow the trend observed in
the MD simulations somewhat closer than the 4D parameterization, even though
the difference is not drastic.

Figure 12: Flattening times of a 12 monolayer nanotip on the 110 surface at
different temperatures. The KMC and MD results are by Jansson et al. [19], and
KMC with machine learning (KMC+ML) refers to this work. The solid lines are
fits of eq. (4) to the data points.

Arising from the event rate formula (1), the mean residence time $t_{0}$ and
the mean migration barrier $E_{\mathrm{a}}$ can be obtained by fitting

$t_{\mathrm{f}}=t_{0}\exp\left(\frac{E_{\mathrm{a}}}{k_{\mathrm{B}}T}\right)$
(4)

to the flattening time data over a range of temperatures. The fitting
parameters for the ANN parameterized KMC model, as well as for the KMC and the
MD results from ref. [19] are tabulated in table 2.

Table 2: Parameters of eq. (4) fitted to the flattening times of a nanotip on the 110 surface with three different models. Model |  $t_{0}$ (s) |  $E_{\mathrm{a}}$ (eV)
---|---|---
KMC+ML | $7.01\cdot 10^{-13}$ | 0.80
KMC [19] | $2.34\cdot 10^{-12}$ | 0.72
MD [19] | $7.33\cdot 10^{-14}$ | 1.00

With the fitted attempt frequency, the flattening of the tips on all three
surfaces was simulated at 1000 K. The mean results from 20 simulations are
tabulated in table 3. In both tables we see much better agreement between the
results obtained in the current work and those obtained with MD simulations.
Fig. 13 shows the initial and the typical final configurations of these
simulations.

Table 3: The mean flattening times of 12 monolayer, 576 atom cuboid nanotips at 1000 K on different surfaces. KMC+ML refers to this work. The error estimates for the KMC+ML data are the standard deviations of 20 simulations with different seeds in the case of the 100 and the 111 surfaces, and standard deviations of 18 simulations in the case of 110 surface—two simulations out of 20 resulted in the destabilisation of the surface. See section 3.2.2 for details. Note that as the $\nu$ parameter is fitted on simulations on the 110 surface across a temperature range from 850 K to 1200 K, the agreement to the MD result is not perfect at the single data point of 1000 K. Surface | KMC+ML (ns) | KMC [19] (ns) | MD [19] (ns)
---|---|---|---
100 | $5.2\pm 0.6$ | $31.0\pm 6.61$ | $1.62\pm 0.60$
110 | $7.5\pm 0.7$ | $9.25\pm 1.10$ | $9.29\pm 1.44$
111 | $4.8\pm 0.9$ | $18.8\pm 0.96$ | $6.01\pm 1.48$

(a) 100 surface

(b) 110 surface

(c) 111 surface

Figure 13: Examples of the initial and final frames of the cuboid nanotip
flattening simulations on the 100, 110 and 111 surfaces at 1000 K. The
simulation was stopped when the tip had reached half of its original height,
replicating the experiment in ref. [19]. The boundary conditions in both
horizontal directions were periodic. Colour coding (available online) is
according to the height coordinate of the atoms.

####  3.2.2 Instability of the {1 1 0} surface

At temperatures $T\geq 950$ K, the 110 sometimes becomes unstable before the
nanotip has time to flatten. In those cases, the surface near the nanotip
starts to deplete, growing larger and larger 100 and 111 facets (see fig. 14).
The atoms accumulate on top of the nanotip, which thus cannot be expected to
flatten. The probability of destabilisation increases with temperature: see
fig. 15 for the probability of roughening as a function of temperature.
Additional 100 simulations on the 110 surface were run for the statistics in
that figure, with a slightly different settings for the simulation box: the
horizontal dimensions were $118\times 130$ Å2, and surface was oriented to
align atomic ridges of the 110 surface with the one of the sides of the
simulation box.

Figure 14: A valley (connected through periodic boundaries) develops in some
simulations at 1000 K or higher temperatures. The system is essentially stuck
in this configuration, with the valley growing deeper with atoms accumulating
from the bottom to the top of the ridge as the simulation progresses. Colour
coding is according to the height coordinate of atoms. The simulation box
borders are shown in black; the box is replicated to illustrate the shape of
the valley. Figure 15: Probability for the 110 surface roughening before the
cuboid nanotip flattens, at different temperatures. 100 simulations were run
at each temperature.

####  3.2.3 Nanoparticle relaxation

Fig. 16 shows examples the shapes of the nanoparticles before and after KMC
relaxation. Each simulation was performed at 1000 K, each nanoparticle was
evolving for about 1–4 µs, depending on the size. The figure is organized in
panels separated by the thin vertical lines. Each panel is named according to
the initial shape of the nanoparticles. Each panel shows three different size
nanoparticles of the same initial shape, which are shown to the left of the
panel and complemented by the final shape of the same nanoparticle after
annealing on the right. The last panel shows the evolution starting from the
minimal surface energy case, which is given by the Wulff construction, as
generated by the Atomic Simulation Environment (ASE) [37]. At finite
temperature, the shape should be expected to fluctuate near the minimum energy
configuration, if the model produces correct surface energies. Since all
particles eventually relaxed into shapes close to the one given by the Wulff
construction, we conclude that our ANN model captures the surface energies
correctly through the migration barrier information, and hence the
thermodynamical evolution of the system given by the model is reliable.

(a) Sphere

$\rightarrow$

(b) 887 atoms

$\rightarrow$

(c) 1409 atoms

$\rightarrow$

(d) 2093 atoms

(e) Octahedron

$\rightarrow$

(f) 891 atoms

$\rightarrow$

(g) 1469 atoms

$\rightarrow$

(h) 2255 atoms

(i) Cube

$\rightarrow$

(j) 864 atoms

$\rightarrow$

(k) 1372 atoms

$\rightarrow$

(l) 2048 atoms

(m) Cuboctahedron

$\rightarrow$

(n) 923 atoms

$\rightarrow$

(o) 1415 atoms

$\rightarrow$

(p) 2057 atoms

(q) Wulff construction

$\rightarrow$

(r) 807 atoms

$\rightarrow$

(s) 1289 atoms

$\rightarrow$

(t) 1925 atoms

Figure 16: Examples of nanoparticle shapes before and after 1–4 microseconds
of KMC simulation at 1000 K. Regardless of the initial shape, the particles
relax to the same final shape that is very close to the minimum energy Wulff
construction.

####  3.2.4 Nanowire stability

We simulated Cu nanowires of 0.5 nm radius at 1000 K. Thin nanowires melt at
temperatures considerably lower than bulk Cu; Granberg et al. found Cu
nanowires of diameter 1.5 nm to melt at 1000 K, using a similar MC/MD-CEM
potential for which the bulk melting points was 1656.72 K [40]. Furthermore,
Cu nanowires are known to fragment by Rayleigh instability at much lower
temperatures than the melting point [41, 42]. Nevertheless, we chose a small
wire radius and a temperature of 1000 K to accelerate the simulation, with the
primary goal of analyzing the performance of the current parameterization of
KMC simulations of Rayleigh instability in thin wires as well as the
mechanisms of such instability.

In fig. 17 we show examples of the ¡110¿ nanowire junction simulations. In
these simulations, atoms diffuse along the wires towards the central junction,
tending to minimise the surface area of the structure. The sharp dips at the
junction are filled up by the atoms from nearby regions of the intersecting
nanowires. When a stable channel to deliver the atoms from the wire to the
knob at the junction is established, the atoms move from large surface-to-
ratio areas at a nanowire to the smaller surface-to-volume regions at the
knob. This leads to thinning of the wires and eventually breaking near the
center. The first breaking occurred right next to the junction in all of the
20 simulations. The mean time for the first breaking to occur in a system of
crossing nanowires was $8\pm 4$ ns. Individual nanowires were much more
stable, taking $250\pm 130$ ns to break by diffusion. Errors are standard
deviations observed in 20 simulations with different random seeds.

Figure 17: Evolution of a system of two crossing ¡110¿ nanowires. Left side of
the figure shows the initial configuration, and the right side shows the
configuration after both wires have split around the central junction. Top
insets show the cross-section of a single wire, and bottom insets show the
central junction.

Fragmentation of Cu nanowire networks at the intersection points was observed
experimentally by Mallikarjuna et al. [38] and Oh et al. [43]. Similar
behaviour was found in Au nanowires by Vigonski et al. [28] both in
experiments and in KMC simulations using the Kimocs code.

##  4 Discussions

The prediction performance of the ANNs in the set of the 3168 flat surface
process barriers (fig. 6) is extremely good. This evidence supports the
applicability of ANNs for surface migration barrier prediction and the
validity of the 26D parameterization for describing the LAEs. The accuracy
observed for the full ~11.7 million barrier set was considerably lower (fig.
10). Some improvement was achieved by dividing the barrier data into three
subsets by surface orientation—100, 110, and 111—and training specialist
networks to predict barriers on each surface. We propose that the improvement
is due to the networks implicitly learning the effect of the different surface
stresses on the migration barriers. Further accuracy was gained by training
ANN ensembles to each barrier subset, and averaging the predicted barrier over
them. Even with these additional techniques, an accuracy comparable to the
flat surface case was not reached.

The reason for the loss of accuracy in the full set is not entirely clear. The
full set is three orders of magnitude larger than the set of flat surface
barriers, with a more heterogeneous composition of the LAEs and a larger range
of barrier values. The data itself seems also to be of reasonable quality, as
the overall agreements with previous barrier sets (figs. 8 and 9) are good.

One possible explanation is, in addition to the different ranges, the
different distributions of barrier values in the sets. The full set
distribution (fig. 7) resembles half of a Gaussian distribution except for the
large amount of zero barriers, while the 3000-barrier flat set distribution
(fig. 5) has no such “discontinuities”. This kind of discontinuity in the
output value distribution may be one of the reasons for the difficulties in
the training. To verify this, in future works it could be considered to assign
some negative pseudo-barriers to the 0 eV cases for the duration of the
training to make the output distribution smoother. During KMC, negative
barriers given by the ANN would be set back to zero. Another method to
smoothen the distribution could be by introducing a finer classification
scheme: in addition to splitting the data set into subsets of 100, 110, and
111 processes, further subdivisions could be made based on the jumping
direction on the surface or on the stability of the LAE (e.g. how many
neighbours the least bound atom in the LAE has).

The complexity of the network architecture could be questioned as a potential
culprit for the lack of accuracy—are single layer ANNs and cascade networks
sophisticated enough to produce sufficient fitting? However, the networks used
in this work are very similar to what has been found to work well in earlier
studies on machine learning for barrier prediction. References [5, 7, 8, 9,
10, 14, 11, 13, 15] use ANNs with only one hidden layer. Cascade networks are
used in refs. [6] and [12]. The earlier rigid lattice variants of the method
also use a similar descriptor: a vector of lattice sites encoded with integers
according to the occupation state of each site. One of the differences in the
earlier studies compared to this work is the size of the LAE—up to hundreds of
atoms instead of the 26 used here. Expanding the LAE could thus provide an
additional way to improve accuracy, although in ref. [16] genetic programming
was successfully applied to predict barriers using LAE with up to 2nn sites in
the limited case of vacancy-assisted diffusion.

Despite the apparent lower accuracy of the ANN in the full barrier set, it
produced physically reasonable results for nanotip flattening, nanoparticle
shape relaxation and simulations of nanowire instability. The attempt
frequency value $2.81\times 10^{14}$ s-1, obtained by fitting the flattening
time of the nanotip on the 110 surface to the flattening time in MD, falls
right within the range of the previously obtained values ($0.7\times 10^{14}$
[19] and $3.1\times 10^{14}$ s-1 [20]). The flattening times for the nanotips
on the 100 and the 111 surfaces are also much closer to the MD results than
with the KMC model proposed in ref. [19] (see table 3), although the
flattening mechanism is different than seen in MD. Due to limitations of the
rigid lattice KMC, the reorientation of the crystal structure, which is
commonly seen within the nanotip oriented in ¡100¿ direction [19, 44, 45], is
not possible within this approach.

As was mentioned in the results section, the thermodynamics of the Cu surface
system appear to be well captured by our model, even though it was not
explicitly trained to do so. Evidence towards this is given both by the
prediction errors of the energy differences (see fig. 11) and the relaxation
of the nanoparticle shapes (fig. 16). Only a small fraction of the jumps have
their energy differences predicted in the opposite order.

Regarding the instability of the 110 surface that was observed at $T>950$ K,
it is worth mentioning that the real Cu 110 surface is known to self-roughen
at high temperatures. The lower bound for the roughening temperature
$T_{\mathrm{R}}$ was experimentally estimated to be at 600 °C (approx. 900 K)
initially by Mochrie [46] and later by Zeppenfeld et al. [47]. Kern estimated
$T_{\mathrm{R}}$ to be 1070 K [48]. Häkkinen et al. found a clear roughening
transition to happen at 1000 K in molecular dynamics simulations [49].

In experiments, the roughening phenomenon is usually observed indirectly from
changes in the diffraction of e.g. X-rays or He-atoms from the surface. In
computational and theoretical models, the roughening is perceived as the
proliferation of atomic steps, as the formation energy of the steps disappears
upon reaching $T=T_{\mathrm{R}}$. In molecular dynamics simulations,
roughening precedes _premelting_ , where the whole surface becomes covered by
a thin liquid-like layer at a temperature below the bulk melting temperature
[50]. While roughening—formation of atomic steps—could in principle be
modelled in rigid lattice KMC, premelting certainly cannot. Further studies
are required as to whether the behaviour observed in this work is a good
description of the roughening phenomenon, or if the transition at 1000 K is
more of a coincidence. In any case, we advise caution when modelling the Cu
110 surface at elevated temperatures with this ANN KMC parameterization.

As for the results for the stability of individual and crossing thin ¡110¿
nanowires, the model suggests that fragmentation occurs first near the
junctions of two nanowires. Individual nanowires took approximately 250 ns to
break, while the first breaking occurred in approximately 8 ns in the systems
of two nanowires. This behaviour is similar to what was earlier shown for Au
wires in both simulations and experiments [28]: nanowire systems tend to break
at the junction points earlier and at lower temperatures than individual
nanowires. The Au and Cu systems are rather well comparable, since both metals
are fcc materials with similar relative surface energies on the low index
facets. Some experimental evidence for similar behaviour in Cu exists:
Mallikarjuna et al. observed increased electrical resistance in Cu nanowire
mesh if photonic welding was carried on excessively long [38]. They attributed
the loss of resistance to breakup of nanowire junctions, which they also saw
in FE-SEM images. Oh et al. also saw broken Cu nanowire junctions after using
very high current in Eddy current welding [43], although they propose that
breaking was due to oxidisation of the junction; in our model, oxygen is not
present. Despite this limitation, the similarity of our results compared with
these experiments and with previous results for Au wires, provides sufficient
validation of the proposed model.

##  5 Conclusions

We have developed an artificial neural network that can predict Cu surface
migration energy barriers with sufficient accuracy and implemented it into a
Kinetic Monte Carlo model. The parameterization is fully three-dimensional,
i.e. it is applicable to arbitrarily rough surfaces. The Kinetic Monte Carlo
model is able to accurately simulate the energy minimisation of Cu
nanoclusters, the thermal flattening of small nanotips as well as the
fragmentation of nanowires. Our model predicts the 110 surface to be unstable
at temperatures above 1000 K, very near to the known roughening temperature of
this Cu surface. The exact mechanism of the instability may be different from
actual roughening, so we advise some caution when using the model at high
temperatures.

## Acknowledgements

J. Kimari was supported by a CERN K-contract and Academy of Finland (Grant No.
313867). V. Jansson was supported by Academy of Finland (Grant No. 285382). F.
Djurabekova acknowledges gratefully the financial support of Academy of
Finland (Grant No. 269696) and MEPhI Academic Excellence Project (Contract No.
02.a03.21.0005). E. Baibuz was supported by a CERN K-contract and the doctoral
program MATRENA of the University of Helsinki. V. Zadin and S. Vigonski were
supported by Estonian Research Council grant PUT 1372. Computing resources
were provided by the Finnish IT Center for Science (CSC) and the Finnish Grid
and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-
infras-2016072533).

## Data availability

The raw data required to reproduce these findings are available to download
from http://urn.fi/urn:nbn:fi:att:cde3019b-d849-41d2-a04a-9f7eb2972f5a. The
processed data required to reproduce these findings are available to download
from http://urn.fi/urn:nbn:fi:att:a4120c3b-9535-405d-b768-4a972bce0b1b.

## References

## References

  * [1] X. Yi, A. E. Sand, D. R. Mason, M. A. Kirk, S. G. Roberts, K. Nordlund, S. L. Dudarev, Direct observation of size scaling and elastic interaction between nano-scale defects in collision cascades, EPL (Europhysics Letters) 110 (3) (2015) 36001.
  * [2] C. Antoine, F. Peauger, F. Le Pimpec, Erratum to: Electromigration occurences and its effects on metallic surfaces submitted to high electromagnetic field: A novel approach to breakdown in accelerators, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 670 (2012) 79–94.
  * [3] CLIC, M. J. Boland, U. Felzmann, P. J. Giansiracusa, T. G. Lucas, R. P. Rassool, C. Balazs, T. K. Charles, K. Afanaciev, I. Emeliantchik, et al., Updated baseline for a staged compact linear collider, CERN Yellow Reports 4 (0).  doi:http://dx.doi.org/10.5170/CERN-2016-004.
URL https://e-publishing.cern.ch/index.php/CYR/article/view/335

  * [4] F. Djurabekova, L. Malerba, C. Domain, C. Becquart, Stability and mobility of small vacancy and copper-vacancy clusters in bcc-fe: An atomistic kinetic monte carlo study, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 255 (1) (2007) 47–51.
  * [5] F. G. Djurabekova, R. Domingos, G. Cerchiara, N. Castin, E. Vincent, L. Malerba, Artificial intelligence applied to atomistic kinetic monte carlo simulations in fe–cu alloys, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 255 (1) (2007) 8–12.
  * [6] N. Castin, R. Pinheiro Domingos, L. Malerba, Use of computational intelligence for the prediction of vacancy migration energies in atomistic kinetic monte carlo simulations, International Journal of Computational Intelligence Systems 1 (4) (2008) 340–352.
  * [7] N. Castin, L. Malerba, Prediction of point-defect migration energy barriers in alloys using artificial intelligence for atomistic kinetic monte carlo applications, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 267 (18) (2009) 3148–3151.
  * [8] N. Castin, L. Malerba, G. Bonny, M. I. Pascuet, M. Hou, Modelling radiation-induced phase changes in binary fecu and ternary fecuni alloys using an artificial intelligence-based atomistic kinetic monte carlo approach, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 267 (18) (2009) 3002–3008.
  * [9] N. Castin, L. Malerba, Calculation of proper energy barriers for atomistic kinetic monte carlo simulations on rigid lattice with chemical and strain field long-range effects using artificial neural networks, The Journal of Chemical Physics 132 (7) (2010) 074507.
  * [10] N. Castin, M. I. Pascuet, L. Malerba, Modeling the first stages of cu precipitation in $\alpha$-fe using a hybrid atomistic kinetic monte carlo approach, The Journal of Chemical Physics 135 (6) (2011) 064502.
  * [11] N. Castin, M. I. Pascuet, L. Malerba, Mobility and stability of large vacancy and vacancy–copper clusters in iron: An atomistic kinetic monte carlo study, Journal of Nuclear Materials 429 (1-3) (2012) 315–324.
  * [12] N. Castin, J. R. Fernández, R. C. Pasianot, Predicting vacancy migration energies in lattice-free environments using artificial neural networks, Computational Materials Science 84 (2014) 217–225.
  * [13] N. Castin, L. Messina, C. Domain, R. Pasianot, P. Olsson, Improved atomistic monte carlo models based on ab-initio-trained neural networks: Application to fecu and fecr alloys, Physical Review B 95 (21) (2017) 214117.
  * [14] M. I. Pascuet, N. Castin, C. S. Becquart, L. Malerba, Stability and mobility of cu–vacancy clusters in fe–cu alloys: A computational study based on the use of artificial neural networks for energy barrier calculations, Journal of nuclear materials 412 (1) (2011) 106–115.
  * [15] L. Messina, N. Castin, C. Domain, P. Olsson, Introducing ab initio based neural networks for transition-rate prediction in kinetic monte carlo simulations, Physical Review B 95 (6) (2017) 064112.
  * [16] K. Sastry, D. D. Johnson, D. E. Goldberg, P. Bellon, Genetic programming for multitimescale modeling, Physical Review B 72 (8) (2005) 085438.
  * [17] S. Verma, T. Rehman, A. Chatterjee, A cluster expansion model for rate constants of surface diffusion processes on ag, al, cu, ni, pd and pt (100) surfaces, Surface Science 613 (2013) 114–125.
  * [18] J. Kimari, V. Jansson, S. Vigonski, E. Baibuz, R. Domingos, V. Zadin, F. Djurabekova, Data sets and trained neural networks for cu migration barriers, submitted for publication (2020).
  * [19] V. Jansson, E. Baibuz, F. Djurabekova, Long-term stability of cu surface nanotips, Nanotechnology 27 (26) (2016) 265708.
  * [20] E. Baibuz, S. Vigonski, J. Lahtinen, J. Zhao, V. Jansson, V. Zadin, F. Djurabekova, Migration barriers for surface diffusion on a rigid lattice: challenges and solutions, Computational Materials Science 146 (2018) 287–302.
  * [21] G. Mills, H. Jónsson, Quantum and thermal effects in h2 dissociative adsorption: evaluation of free energy barriers in multidimensional quantum systems, Physical review letters 72 (7) (1994) 1124.
  * [22] G. Mills, H. Jónsson, G. K. Schenter, Reversible work transition state theory: application to dissociative adsorption of hydrogen, Surface Science 324 (2) (1995) 305–337.
  * [23] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of computational physics 117 (1) (1995) 1–19.
  * [24] M. S. Stave, D. E. Sanders, T. J. Raeker, A. E. DePristo, Corrected effective medium method. v. simplifications for molecular dynamics and monte carlo simulations, The Journal of Chemical Physics 93 (6) (1990) 4413–4426.
  * [25] S. B. Sinnott, M. S. Stave, T. J. Raeker, A. E. DePristo, Corrected effective-medium study of metal-surface relaxation, Physical Review B 44 (16) (1991) 8927.
  * [26] A. Stukowski, Visualization and analysis of atomistic simulation data with ovito–the open visualization tool, Modelling and Simulation in Materials Science and Engineering 18 (1) (2009) 015012.
  * [27] J. Zhao, E. Baibuz, J. Vernieres, P. Grammatikopoulos, V. Jansson, M. Nagel, S. Steinhauer, M. Sowwan, A. Kuronen, K. Nordlund, et al., Formation mechanism of fe nanocubes by magnetron sputtering inert gas condensation, ACS nano.
  * [28] S. Vigonski, V. Jansson, S. Vlassov, B. Polyakov, E. Baibuz, S. Oras, A. Aabloo, F. Djurabekova, V. Zadin, Au nanowire junction breakup through surface atom diffusion, Nanotechnology 29 (1) (2017) 015704.
  * [29] V. Jansson, E. Baibuz, A. Kyritsakis, S. Vigonski, V. Zadin, S. Parviainen, A. Aabloo, F. Djurabekova, Growth mechanism for nanotips in high electric fieldsSubmitted for publication.  arXiv:https://arxiv.org/abs/1909.05825.
  * [30] V. Jansson, A. Kyritsakis, S. Vigonski, E. Baibuz, V. Zadin, A. Aabloo, F. Djurabekova, Tungsten migration energy barriers for surface diffusion: a parameterization for kmc simulations, Modelling and Simulation in Materials Science and Engineering 28 (3) (2020) 035011.
  * [31] O. Trushin, A. Karim, A. Kara, T. S. Rahman, Self-learning kinetic monte carlo method: application to cu (111), Physical Review B 72 (11) (2005) 115401.
  * [32] A. Latz, L. Brendel, D. E. Wolf, A three-dimensional self-learning kinetic monte carlo model: application to ag (111), Journal of Physics: Condensed Matter 24 (48) (2012) 485005.
  * [33] S. Nissen, Implementation of a fast artificial neural network library (fann), Tech. rep., Department of Computer Science University of Copenhagen (DIKU) (2003).
URL http://fann.sf.net

  * [34] E. Jones, T. Oliphant, P. Peterson, et al., SciPy: Open source scientific tools for Python (2001–).
URL http://www.scipy.org/

  * [35] C. Igel, M. Hüsken, Improving the rprop learning algorithm, in: Proceedings of the second international ICSC symposium on neural computation (NC 2000), Vol. 2000, Citeseer, 2000, pp. 115–121.
  * [36] S. E. Fahlman, C. Lebiere, The cascade-correlation learning architecture, in: Advances in neural information processing systems, 1990, pp. 524–532.
  * [37] A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng, K. W. Jacobsen, The atomic simulation environment—a python library for working with atoms, Journal of Physics: Condensed Matter 29 (27) (2017) 273002.
URL http://stacks.iop.org/0953-8984/29/i=27/a=273002

  * [38] K. Mallikarjuna, H.-J. Hwang, W.-H. Chung, H.-S. Kim, Photonic welding of ultra-long copper nanowire network for flexible transparent electrodes using white flash light sintering, RSC advances 6 (6) (2016) 4770–4779.
  * [39] E. Baibuz, S. Vigonski, J. Lahtinen, J. Zhao, V. Jansson, V. Zadin, F. Djurabekova, Data sets of migration barriers for atomistic Kinetic Monte Carlo simulations of Cu self-diffusion via first nearest neighbour atomic jumps, Data in Brief 17 (2018) 739–743.
  * [40] F. Granberg, S. Parviainen, F. Djurabekova, K. Nordlund, Investigation of the thermal stability of cu nanowires using atomistic simulations, Journal of Applied Physics 115 (21) (2014) 213518.
  * [41] M. Toimil Molares, A. Balogh, T. Cornelius, R. Neumann, C. Trautmann, Fragmentation of nanowires driven by rayleigh instability, Applied physics letters 85 (22) (2004) 5337–5339.
  * [42] H. Li, J. M. Biser, J. T. Perkins, S. Dutta, R. P. Vinci, H. M. Chan, Thermal stability of cu nanowires on a sapphire substrate, Journal of Applied Physics 103 (2) (2008) 024315.
  * [43] J. S. Oh, J. S. Oh, T. H. Kim, G. Y. Yeom, Efficient metallic nanowire welding using the eddy current method, Nanotechnology 30 (6) (2018) 065708.
  * [44] H. S. Park, K. Gall, J. A. Zimmerman, Shape memory and pseudoelasticity in metal nanowires, Physical Review Letters 95 (25) (2005) 255504.
  * [45] W. Liang, M. Zhou, F. Ke, Shape memory effect in cu nanowires, Nano Letters 5 (10) (2005) 2039–2043.
  * [46] S. Mochrie, Thermal roughening of the copper (110) surface: An x-ray diffraction experiment, Physical review letters 59 (3) (1987) 304.
  * [47] P. Zeppenfeld, K. Kern, R. David, G. Comsa, No thermal roughening on cu (110) up to 900 k, Physical review letters 62 (1) (1989) 63.
  * [48] K. Kern, Thermal dynamics of (110) fcc metal surfaces, in: Surface X-Ray and Neutron Scattering, Springer, 1992, pp. 69–71.
  * [49] H. Häkkinen, J. Merikoski, M. Manninen, J. Timonen, K. Kaski, Roughening of the cu (110) surface, Physical review letters 70 (16) (1993) 2451.
  * [50] H. Häkkinen, M. Manninen, Computer simulation of disordering and premelting of low-index faces of copper, Physical review B 46 (3) (1992) 1725.


# Benchmarking Amplitude Estimation on a Superconducting Quantum Computer

Salvatore Certo<EMAIL_ADDRESS>Deloitte Consulting, LLP  Anh Dung Pham
Deloitte Consulting, LLP  Daniel Beaulieu  Deloitte Consulting, LLP

###### Abstract

Amplitude Estimation (AE) is a critical subroutine in many quantum algorithms,
allowing for a quadratic speedup in various applications like those involving
estimating statistics of various functions as in financial Monte Carlo
simulations. Much work has gone into devising methods to efficiently estimate
the amplitude of a quantum state without expensive operations like the Quantum
Fourier Transform (QFT), which is especially prohibitive given the constraints
of current NISQ devices. Newer methods have reduced the number of operations
required on a quantum computer and are the most promising near-term
implementations of the AE subroutine. While it remains to be seen the exact
circuit requirements for a quantum advantage in applications relying on AE, it
is necessary to continue to benchmark the algorithm’s performance on current
quantum computers and the circuit costs associated with such subroutines.

Given these considerations, we expand on results from previous experiments in
using Maximum Likelihood Estimation (MLE) to approximate the amplitude of a
quantum state and provide empirical upper bounds on the current feasible
circuit depths for AE on a superconducting quantum computer. Our results show
that MLE using optimally-compiled circuits can currently outperform naive
sampling for up to 3 Grover Iterations with a circuit depth of 131, which is
higher than reported in other experimental results. This functional benchmark
is one of many that will be continually monitored against current quantum
hardware to measure the necessary progress towards quantum advantage.

##  I Introduction

Many algorithms have been created which can provide a quadratic speedup of
certain calculations montanaro , including those in finance like Monte Carlo
simulations for derivative pricing and Value at Risk (VaR). These algorithms
involve applying functions to one or more random variables and extracting the
expected value of the function which is the quantity of interest. For example,
finding the fair price for a European Call Option includes evaluating the
payoff function based on the strike price and the distribution of the expected
spot prices for the asset at maturity. In a quantum setting, this example
would involve several steps, including loading the spot price distribution
into a quantum state, computing the payoff function, and finally extracting
the expected value from the amplitude of a quantum state egger2019credit W19
S19 O19 .

While these steps leverage unique properties of quantum mechanics, it is only
in the last step of extracting the expected value where the quadratic speedup
is obtained. The retrieval of the expected value of these computations rely on
Quantum Amplitude Estimation, which converges as $O(1/M)$ compared to the
classical convergence rate of $O(1/\sqrt{M})$montanaro . Therefore there is a
motivation to perform Amplitude Estimation (AE) in an efficient and accurate
way, and to benchmark its success on current hardware. Formally, the goal of
AE is to provide an approximation of $a$ prepared by an operator $\mathcal{A}$
such as in

$\mathcal{A}\mathinner{|{0_{n}}\rangle}\mathinner{|{0}\rangle}=\sqrt{1-a}\mathinner{|{\psi_{0}}\rangle}_{n}\mathinner{|{0}\rangle}+\sqrt{a}\mathinner{|{\psi_{1}}\rangle}_{n}\mathinner{|{1}\rangle}$
(1)

In other words, there is some set of quantum operations $\mathcal{A}$ (e.g.
the payoff function of a call option) that maps a quantity of interest (e.g.
the expected value of the option) into a quantum state
$\mathinner{|{\psi_{1}}\rangle}\mathinner{|{1}\rangle}$ and the goal is to
find the amplitude of that state, $\sqrt{a}$. The canonical method tapp  was
to use Quantum Phase Estimation (QPE) to provide an estimate for the
amplitude. This approach relied on using a controlled version of the grover
iterate $\mathcal{Q}$, where $\mathcal{Q}$ is defined as

$\mathcal{Q}=\mathcal{A}\mathcal{S}_{0}\mathcal{A^{\dagger}}\mathcal{S}_{\psi}$

Here
$\mathcal{S}_{\psi}=\mathcal{I}-2\mathinner{|{\psi_{1}}\rangle}_{n}\mathinner{\langle{\psi_{1}}|}_{n}\bigotimes\mathinner{|{1}\rangle}\mathinner{\langle{1}|}$
and
$\mathcal{S}_{0}=\mathcal{I}-2\mathinner{|{0}\rangle}_{n+1}\mathinner{\langle{0}|}_{n+1}$.
In other words, $\mathcal{S}_{\psi}$ adds a negative sign to the
$\mathinner{|{\psi_{1}}\rangle}$ states and $\mathcal{S}_{0}$ reflects around
the zero state.

By then performing Quantum Fourier Transform (QFT) on the control register of
size $m$ we would get an estimate for $a$.

Figure 1: Canonical algorithm for Amplitude Estimation by Quantum Phase
Estimation.

This method has features that make it unsuitable for current Noisy
Intermediate Scale Quantum devices. The qubit topography of many quantum
systems make controlled gates expensive and can quickly increase the depth of
the circuit. Similarly, QFT requires expensive operations that are generally
not practically implementable in current devices.

There has been active research into other ways to conduct AE, some even
involving variational methods plekhanov2021variational . Most current methods
involve non-controlled applications of $\mathcal{Q}$ and some classical post
processing to create an estimate of $a$ grinko2019iterative S20 . These
methods rely on the fact that
$\mathcal{A}\mathinner{|{0_{n}}\rangle}=\cos{\theta}\mathinner{|{\psi_{0}}\rangle}+\sin{\theta}{\mathinner{|{\psi_{1}}\rangle}}$.
If we can get an accurate estimate for $\theta$, then $\sin^{2}{\theta}=a$.
Included in these methods are the Maximum Likelihood Amplitude Estimation
(MLAE) algorithm S20 , which applies successive applications of the Grover
Iterate $\mathcal{Q}$ and finds the $\theta$ that is most likely given the
results. A benefit of MLAE is that it allows us to control the number of
applications as well as the number of circuit evaluations for the
$\mathcal{Q}$ oracle.

Given a number $m_{k}$ of repetitions of $\mathcal{Q}$, the number of circuit
evaluations $N_{k}$, and the number of times $h_{k}$ that
$\mathinner{|{\psi_{1}}\rangle}\mathinner{|{1}\rangle}$ was measured, the
probability of measuring
$\mathinner{|{\psi_{1}}\rangle}\mathinner{|{1}\rangle}$ is
$P(m_{k};a)=\sin^{2}{((2m_{k}+1)}\theta_{a})$. The likelihood function is
therefore

$\mathcal{L}(h;a)=\prod_{k=0}^{M}[P(m_{k};a)]^{h_{k}}[1-P(m_{k};a)]^{N_{k}-h_{k}}$
(2)

The number of repetitions of $\mathcal{Q}$ as well as the number of circuit
evaluations per each repetition define the sequencing schedule for the
algorithm. Different schedules have been proposed, including the Linearly
Incremental Sequence (LIS), the Exponential Incremental Sequence (EIS), and
the Power Law schedule. The experiments we conducted used the simple LIS
schedule.

##  II Motivation for Experiment

Experimental results were recently published lowdepth  that showed evidence of
amplitude estimation outperforming naive sampling and low error rates on a
state-of-the-art trapped ion quantum computer. It is known that the different
types of quantum computers have advantages and disadvantages, and while some
work has been published on benchmarking AE on superconducting hardware
rao2020quantum , there has been little published on actual circuit depths and
their effects on amplitude estimation. While quantum engineers use high level
languages to construct circuits, they are first compiled to solve the many
constraints a particular quantum computer has, including qubit topography and
actual basis gate implementations. This can have dramatic effects on the
actual circuit depths, number and types of gates used, thereby altering the
performance of the algorithm than what may be expected based on the non-
compiled circuit. Therefore it is critical to continue to evaluate the status
of current computers as well as the effect of circuit depth on their
performance, on which we report in this paper.

##  III Amplitude Estimation Circuits used in MLE Algorithm

One can construct any arbitrary circuit to use for amplitude estimation, the
only requirement being that there exists a well-defined quantum state from
which we can benchmark the results. For our purposes, we wanted to design an
experiment that is reproduce-able, scalable, and has advantageous properties
related to the algorithm we are benchmarking. The circuits used in lowdepth
share some of these properties. They originally were proposed to be used as a
unary data loader, a method to load unit vectors and matrices into a quantum
system. As it pertains to amplitude estimation and shown in lowdepth , it also
has nice properties for both the oracle $\mathcal{S}_{\psi}$ and diffuser
$\mathcal{A}\mathcal{S}_{0}\mathcal{A^{\dagger}}$ of the Grover Iterate,
requiring single Z gates for both $\mathcal{S}_{\psi}$ and $\mathcal{S}_{0}$.

The experiments previously mentioned used unary data loaders to load two unit
vectors and compute their inner product in the amplitude of a quantum state
that is to be estimated. It relied on the $\mathcal{RBS}$ gate which is a two
qubit unitary matrix that is described as:

$\mathcal{RBS}(\theta)=\begin{pmatrix}1&0&0&0\\\
0&\cos{\theta}&\sin{\theta}&0\\\ 0&-\sin{\theta}&\cos{\theta}&0\\\
0&0&0&1\end{pmatrix}$ (3)

Based on the unit vectors, there is a method for finding $\theta$. While we
refer to lowdepth  and johri2020nearest  for more information on constructing
the specific gates, we will briefly provide an example of using them. Given
two unit vectors, $\mathinner{|{a}\rangle}$ and $\mathinner{|{b}\rangle}$, we
can find angles through the procedure in johri2020nearest  to construct the
circuit below that will compute the square of their inner product
$\braket{a}{b}$ in the amplitude of the first qubit.

Figure 2: Circuit for computing the inner product between two 2 dimensional
unit vectors.

For instance,

$\mathinner{|{a}\rangle}=\begin{bmatrix}-0.96184207&0.27360523\end{bmatrix}$

$\mathinner{|{b}\rangle}=\begin{bmatrix}0.96845237&0.24919873\end{bmatrix}$

$\braket{a}{b}^{2}=\mathinner{|{01}\rangle}=.745314789774813$

We used similar initial gate configurations in Qiskit as lowdepth , however
both the hardware we ran our circuits on as well as the compilers changed the
gate configurations. Still, shown in Figure 3 and Figure 4 are the initial
circuits for the 4 dimensional inner product circuits which represent
$\mathcal{A}$ and
$\mathcal{A}\mathcal{S}_{0}\mathcal{A^{\dagger}}\mathcal{S}_{\psi}\mathcal{A}$.

Figure 3: Inner product estimation circuit for $\mathcal{A}$ Figure 4: Circuit
for a single grover iteration of
$\mathcal{A}\mathcal{S}_{0}\mathcal{A^{\dagger}}\mathcal{S}_{\psi}\mathcal{A}$

In the 4 dimensional case, $\mathcal{A}$ will load the states
$a\mathinner{|{0001}\rangle}+b\mathinner{|{0010}\rangle}+c\mathinner{|{0100}\rangle}+d\mathinner{|{1000}\rangle}$
where $a,b,c,d$ are the amplitudes and $a^{2}+b^{2}+c^{2}+d^{2}=1$. We are
interested in the amplitude of $a$, the probability of the
$\mathinner{|{0001}\rangle}$ state which represents the squared inner product
of the 4 dimensional unit vectors.

While we tested a trivial example on 2 qubits, the main result of our work was
in estimating the inner product of 2 random 4 dimensional unit vectors. In
this case, we constructed $\mathcal{A}$ and
$\mathcal{A}\mathcal{S}_{0}\mathcal{A^{\dagger}}\mathcal{S}_{\psi}\mathcal{A}$
as shown above. The quantum computer we ran our circuits on took the high
level circuits and compiled them to the allowable basis gates. We then used
Qiskit and Pytket to optimize these compiled circuits as aggressively as
allowed.

To estimate the amplitude of the first qubit which represents the inner
product of the unit vectors, we added Grover Iterates from $(0...7)$ where $0$
is just $\mathcal{A}$. Once compiled based on the 3 different methodologies,
we ran the circuit 500 times and used the results as input to our Maximum
Likelihood Estimation algorithm. Therefore by using $m_{k}=7$ of repetitions
of $\mathcal{Q}$, we also used the previous results of $m_{k}=(0...6)$ for our
calculation.

We generated random unit vectors whose inner product is uniformly distributed
and instantiated our initial distribution for $\theta$ in the MLE algorithm to
be a uniform distribution between $[0,\frac{\pi}{2}]$. For each set of unit
vector pairs, we compared the MLE output and traditional sampling method to
the exact inner product, giving us our estimate of relative error.

##  IV IBM Nairobi Specifications

The experiments conducted were done using the IBM Nairobi superconducting
computer. It has 7 qubits with a quantum volume of 32. The exact
specifications are detailed in the specifications table.

Figure 5: Depth, gate count, and gate composition of various compilation
methods

##  V Comparison of Optimized Circuits

Figure 6: Comparison of $\mathcal{A}$ circuits.

The RBS gates outlined above have significant depths when compiled on IBM
hardware. In fact, they can reach a depth of over 500 for 7 Grover Iterations.
Knowing the need to reduce these depths and optimize the circuits as much as
possible, we explored two well known open source compilers provided by Qiskit
Qiskit  and Pytket pytket . There are many ways to optimally compile a quantum
circuit, none of which are trivial. Considerations need to be made for the
goal of the algorithm, the specific hardware, and the types of gates required
along with their typical error rates.

Both Qiskit and Pytket have a similar set of compilation passes, first
translating the high level circuit into the required basis gates of the
hardware, then attempting to optimize the circuit to reduce depth and noise.
As referenced in pytket , Pytket’s optimization compiler falls into two broad
categories: peephole optimizations and macroscopic optimizations. Peephole
optimizations are similar to those in classical settings, where a sliding
window looks for local patterns that can be substituted with equivalent sub-
circuits with lower gate counts or depth. Macroscopic optimization looks for
non-local structures that can be synthesized to increase parallelization or
fewer gates, potentially leading to another pass of peephole optimization. We
refer to pytket  for more detailed information.

In figure 6 we show a sample of circuits for certain 4 dimensional vectors
with no Grover iterations. Both Qiskit and Pytket substantially reduced the
two-qubit interactions, gate counts, and overall circuit depth. In figure 5 we
detail the number of gates and circuit depth per Grover iteration, as well as
the total (sum of all 7 sub-circuits) counts per gate type. We can see that
both Qiskit and Pytket reduced the number of gates and the circuit depth; the
circuit depth was reduced by almost half and was a substantial improvement
over the standard compilation procedure.

Furthermore, both found the minimum number of CNOT gates that could be used,
while the standard compiler used more than double. We leave it to further
experiments to benchmark the different ways to optimize quantum circuits; the
results we present are not direct comparisons between either approaches. If we
were to do so, we would choose a variety of metrics, algorithms, and circuit
types to quantitatively measure them. Here we are simply showing in this
specific case both are worthwhile to use versus the standard compiler.

##  VI Minimal Depth Amplitude Estimation Results

Here we will review the results of estimating the amplitude of random unit
vectors for both 2 and 4 dimensions. We used a linear sequencing schedule for
MLE, of 500 shots each up to 7 Grover Iterations. The experiment was repeated
50 times each, and the results averaged. For naive sampling, we ran the
circuits for up to 3000 shots in increments of 500.

For the 2 qubit experiments, both the ordinarily-compiled circuits as well as
the optimized circuits outperformed naive sampling. Interestingly, the larger
number of shots for naive sampling did not reduce the error, which is most
likely due to the fact that the inherent noise found on current hardware
prohibits better amplitude estimates. As shown in Figure 7, The lowest
recorded mean error of .0017 was after 7 Grover iterations, although after 3
Grover iterations the decrease in error became minimal.

While the 2 qubit experiment is trivially small, it demonstrated that MLE does
outperform simple averaging of samples and after even 1 Grover iteration the
error decreased by half. It also showed the impact of reducing circuit depth,
complexity, and optimizing it for a particular quantum computer.

Figure 7: 2 qubit AE results for ordinary transpiled and optimized depth
circuits. The circles represent the maximum number of Grover Iterations used
in MLE, with 500 shots per each iteration. Figure 8: 4 qubit AE results for
ordinary compiled and optimized depth circuits. The optimized circuits were
done using level 2 optimization provided by Pytket and level 3 optimization by
Qiskit. The best results achieved in our experiments were found after 3 Grover
iterations on the Pytket optimally-compiled circuits, achieving a mean error
of .0442.

The 4 qubit results are more significant as the experiment still showed MLE
outperforming naive sampling with much larger circuit depths and 2-qubit
gates. The circuit compiler provided by both Qiskit and Pytket reduced the
error rates for both simple averaging and MLAE, with Pytket providing a
shallower circuit with lower mean error. Even with an optimal circuit layout,
the depths reached a maximum of around 300 and the errors did deteriorate
after 3 Grover iterations caused by decoherence of the qubits and compounded
errors throughout the circuit. Both Qiskit and Pytket have different levels of
circuit compilers, ranging from (0,1,2,3) for Qiskit and (0,1,2) for Pytket.
Each of these compilation levels differ in their aggressiveness for optimizing
the ciruits and reducing the depth. We used the most aggressive compilers from
each and the results are shown in Figure 8.

It is worthy to note that while Pytket seemed more robust to noise in the
higher number of Grover Iterations, this could be the result of several
causes. We used 50 random vectors for each experiment, and so the exact
distribution of this small sample size could account for the discrepancy.
Qiskit had a lower sampling error, which we also believe could have been the
result of noise in the quantum system as well as the random vectors.
Nonetheless, both offer substantial improvement over non-optimized circuits
and show how significant an impact a well-compiled circuit performs over an
ordinarily-compiled one. Even the results from the naive averaging of the non-
optimized circuits had a much larger variance than the optimized circuits from
both Qiskit and Pytket, showing that the compilation provided by both
substantially reduces unwanted noise and stabilizes the results.

Compared to lowdepth , our optimized circuits had higher depths than the
experiments performed on a trapped ion computer, which reported two-qubit
gates up to ninety two and a maximum circuit depth of sixty two. The depths
reported here are much larger compared to previous experiments as we were able
to achieve a minimum mean error of around .04 with circuit depths of 131 for 3
Grover iterations. In fact, while the circuits constructed were minimal with
respect to IBM’s superconducting computers, they are some of the deepest
circuits reported thus far that still obtained a significant reduction in
error while performing MLAE. That said, recent work has determined that the
threshold for quantum advantage in certain Monte Carlo simulations requires a
t-depth of around 54 million threshold , and so while current hardware is out
of reach for such performance, we hope that functional benchmark tests such as
this experiment continue to evaluate the maturity of current systems in the
pursuit of such an advantage.

##  VII Conclusion

In this paper we have expanded on previous work lowdepth  in using circuits
that compute the inner product between unit vectors to benchmark the Maximum
Likelihood Amplitude Estimation algorithm. While we have shown that both the
circuit optimizers by Qiskit and Pytket work very well in reducing noise and
increasing the allowable number of Grover iterations we can perform in the
algorithm, we have more importantly provided an empirical upper bound on the
maximum allowable depth for AE on current hardware, a metric that will be will
closely monitored as the industry progresses.

##  VIII Disclaimer

About Deloitte: Deloitte refers to one or more of Deloitte Touche Tohmatsu
Limited (“DTTL”), its global network of member firms, and their related
entities (collectively, the “Deloitte organization”). DTTL (also referred to
as “Deloitte Global”) and each of its member firms and related entities are
legally separate and independent entities, which cannot obligate or bind each
other in respect of third parties. DTTL and each DTTL member firm and related
entity is liable only for its own acts and omissions, and not those of each
other. DTTL does not provide services to clients. Please see
www.deloitte.com/about to learn more.

Deloitte is a leading global provider of audit and assurance, consulting,
financial advisory, risk advisory, tax and related services. Our global
network of member firms and related entities in more than 150 countries and
territories (collectively, the “Deloitte organization”) serves four out of
five Fortune Global 500® companies. Learn how Deloitte’s approximately 330,000
people make an impact that matters at www.deloitte.com. This communication
contains general information only, and none of Deloitte Touche Tohmatsu
Limited (“DTTL”), its global network of member firms or their related entities
(collectively, the “Deloitte organization”) is, by means of this
communication, rendering professional advice or services. Before making any
decision or taking any action that may affect your finances or your business,
you should consult a qualified professional adviser. No representations,
warranties or undertakings (express or implied) are given as to the accuracy
or completeness of the information in this communication, and none of DTTL,
its member firms, related entities, employees or agents shall be liable or
responsible for any loss or damage whatsoever arising directly or indirectly
in connection with any person relying on this communication. Copyright © 2021.
For information contact Deloitte Global.

## References

  * [1] Qiskit: An open-source framework for quantum computing, 2021.
  * [2] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  Quantum Computation and Information, page 53–74, 2002.
  * [3] Shouvanik Chakrabarti, Rajiv Krishnakumar, Guglielmo Mazzola, Nikitas Stamatopoulos, Stefan Woerner, and William J. Zeng.  A threshold for quantum advantage in derivative pricing.  Quantum, 5:463, Jun 2021.
  * [4] Daniel J Egger, Ricardo García Gutiérrez, Jordie Cahue Mestre, and Stefan Woerner.  Credit risk analysis using quantum computers.  IEEE Transactions on Computers, 2020.
  * [5] Tudor Giurgica-Tiron, Sonika Johri, Iordanis Kerenidis, Jason Nguyen, Neal Pisenti, Anupam Prakash, Ksenia Sosnova, Ken Wright, and William Zeng.  Low depth amplitude estimation on a trapped ion quantum computer, 2021\.
  * [6] Dmitry Grinko, Julien Gacon, Christa Zoufal, and Stefan Woerner.  Iterative quantum amplitude estimation.  npj Quantum Information, 7(1):1–6, 2021.
  * [7] Sonika Johri, Shantanu Debnath, Avinash Mocherla, Alexandros Singh, Anupam Prakash, Jungsang Kim, and Iordanis Kerenidis.  Nearest centroid classification on a trapped ion quantum computer, 2020\.
  * [8] Ashley Montanaro.  Quantum speedup of monte carlo methods.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2181):20150301, Sep 2015.
  * [9] Román Orús, Samuel Mugel, and Enrique Lizaso.  Quantum computing for finance: overview and prospects.  Reviews in Physics, page 100028, 2019.
  * [10] Kirill Plekhanov, Matthias Rosenkranz, Mattia Fiorentini, and Michael Lubasch.  Variational quantum amplitude estimation, 2021.
  * [11] Pooja Rao, Kwangmin Yu, Hyunkyung Lim, Dasol Jin, and Deokkyu Choi.  Quantum amplitude estimation algorithms on ibm quantum devices, 2020.
  * [12] Seyon Sivarajah, Silas Dilkes, Alexander Cowtan, Will Simmons, Alec Edgington, and Ross Duncan.  t—ket⟩: a retargetable compiler for nisq devices.  Quantum Science and Technology, 6(1):014003, Nov 2020.
  * [13] Nikitas Stamatopoulos, Daniel J Egger, Yue Sun, Christa Zoufal, Raban Iten, Ning Shen, and Stefan Woerner.  Option pricing using quantum computers.  Quantum, 4:291, 2020.
  * [14] Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tomoki Tanaka, Tamiya Onodera, and Naoki Yamamoto.  Amplitude estimation without phase estimation.  Quantum Information Processing, 19(2):75, 2020.
  * [15] Stefan Woerner and Daniel J Egger.  Quantum risk analysis.  npj Quantum Information, 5(1):1–8, 2019.