# Binary Operations on Neuromorphic Hardware with Application to
Linear Algebraic Operations and Stochastic Equations

Oleksandr Iaroshenko  Andrew T. Sornborger<EMAIL_ADDRESS>1Information
Sciences, CCS-3, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

###### Abstract

Non-von Neumann computational hardware, based on neuron-inspired, non-linear
elements connected via linear, weighted synapses – so-called neuromorphic
systems – is a viable computational substrate. Since neuromorphic systems have
been shown to use less power than CPUs for many applications, they are of
potential use in autonomous systems such as robots, drones, and satellites,
for which power resources are at a premium. The power used by neuromorphic
systems is approximately proportional to the number of spiking events produced
by neurons on-chip. However, typical information encoding on these chips is in
the form of firing rates that unarily encode information. That is, the number
of spikes generated by a neuron is meant to be proportional to an encoded
value used in a computation or algorithm. Unary encoding is less efficient
(produces more spikes) than binary encoding. For this reason, here we present
neuromorphic computational mechanisms for implementing binary two’s complement
operations. We use the mechanisms to construct a neuromorphic, binary matrix
multiplication algorithm that may be used as a primitive for linear
differential equation integration, deep networks, and other standard
calculations. We also construct a random walk circuit and apply it in Brownian
motion simulations. We study how both algorithms scale in circuit size and
iteration time.

##  I Introduction

In the context of neural systems, where, due to the expense of synaptic
transmission, a spike costs more than the absence of a spike, sparse
communication is more efficient. This fact has been studied both in biological
neural systems Harris et al. (2012); Levy and Baxter (2002) as well as in
neuromorphic systems Davies .

However, biological neural computations as well as algorithms implemented on
neuromorphic chips typically use some form of unary coding, i.e. information
is encoded in firing rates. Though energy efficient at low spiking rates,
unary coding can be wasteful in terms of the number of spikes used for a given
numerical resolution. Each wasted spike contributes excess power used by the
neural circuit.

A previous attempt to use more efficient spike coding in a neuromorphic
machine learning algorithm used a mixed binary-unary encoding to make
multiplication more efficient Diehl et al. (2016). Here, we present three
implementations of circuits using binary operations: 1) a fully binary
implementation of multiplication, where a two’s complement number encoded in
spikes multiplies a two’s complement number encoded in synapses, 2) an
extension of this circuit to vector-matrix multiplication, and 3) a simpler,
and highly parallelizable, circuit for performing a random walk.

Our neural circuits use a pulse-gating framework to control and coordinate the
flow of spikes between a number of cores. They are made up of five neuron
types and standard synaptic connections are used for the implementation.

By encoding numbers in binary, we logarithmically reduce the number of spikes
needed to represent a binary-encoded number or vector. This reduction comes at
the cost of increasing the neurosynaptic footprint needed to represent the
neural circuit.

Matrix multiplication has a broad range of uses across the spectrum of
computational algorithms. The neural circuits that we develop below may be
used for various applications, including implementing kernels for machine
learning to more prosaic computational algorithms. Here, we consider 1)
simulations of linear dynamics, and 2) simulations of Brownian motion. The
first application is straightforward with the availability of the matrix-
vector multiplication circuit that we demonstrate. We implement the second
application in an easily parallelizable neural circuit.

Both of the applications that we consider have been studied in neuromorphic
circuits previously Severa et al. (2018); Parekh et al. (2017); Reeder et al.
(2019). The circuits we present here use two’s complement binary operations
and scale poly-logarithmically in the number of neurons as a function of
numerical resolution, whereas the modular spatial code using a residue
numerical system presented in Severa et al. (2018) offers a compromise between
unary and binary codings with intermediate scaling. The time scaling for our
circuits is also poly-logarithmic in the numerical resolution, whereas with
unary encoding, scaling would be linear in the numerical resolution (i.e.
exponentially more expensive than binary). The neural circuits that we present
here should prove useful in providing low-spiking rate, low-power algorithms
for linear algebraic applications in neuromorphic circuits.

##  II Methods

###  II.1 Neural dynamics

We implemented our two’s complement algorithms on Intel’s neuromorphic
research processor codenamed Loihi Davies et al. (2018). A single Loihi chip
consists of $128$ logical cores with each core having $1024$ compartments and
an upper bound of $2^{14}=16,384$ synapses per core. A simple neuron occupies
one compartment whereas more sophisticated neurons can be built from several
compartments. The neural dynamics is computed at discrete timesteps called
ticks.

For our binary multiplication circuit, we make use of three types of simple
neurons:

  * •

$\gamma_{0}$ – a TRANSFER neuron, which emits a spike one tick after it
receives an incoming spike from any of its inputs;

  * •

$\gamma_{1}$ – a GATING neuron, which emits a spike one tick after receiving
at least two incoming spikes at the same tick;

  * •

$\gamma_{2}$ – a COUNTER neuron, which emits a number of spikes equal to the
number of incoming spikes at time $t_{0}$, at a rate of one spike per tick
starting from time $t_{0}+1$.

Loihi’s neuronal compartments convert incoming spikes to a current which is
multiplied by a synaptic weight and accumulated as an internal voltage:

$\begin{cases}u_{i}(t)=u_{i}(t-1)\cdot(2^{12}-\delta_{i}^{(u)})\cdot
2^{-12}+2^{6}\cdot\sum_{j}w_{ij}s_{j}(t),\\\
v_{i}(t)=v_{i}(t-1)\cdot(2^{12}-\delta_{i}^{(v)})\cdot
2^{-12}+u_{i}(t)+u_{j}^{\textup{bias}}.\end{cases}$ (1)

Here $i$ and $j$ are the indices of the receiving and emitting neurons
(resp.), $u$ is the current in the receiving neuron, $v$ is the voltage in the
receiving neuron. $\delta_{i}^{(u,v)}\in[0,2^{12}]$ are current and voltage
decay constants respectively, $w_{ij}$ is the synaptic weight,
$s_{j}(t)\in\\{0,1\\}$ are upstream emitted spikes, and
$u_{j}^{\textup{bias}}$ is a constant bias current. The variable
$t\in\mathbb{N}$ represents time (in units of ticks).

We use the same synaptic weight throughout, $w_{ij}\equiv w=100$, and a single
spike current $u_{s}=2^{6}$ (by default in the Loihi API) as shown in equation
(1). When the voltage exceeds a threshold, $v_{\textup{th}}$, a spike is
emitted and the voltage is reset to zero. Both the incoming current and the
voltage can have their own exponential decay rules defined by the decay
constants $\delta_{i}^{(u,v)}$ and the system of equations (1). In all neurons
$\gamma_{j}$, $j=0,1,2$, the voltage decay is immediate
($\delta_{i}^{(v)}=2^{12}$) and resets to zero at every tick. The current
decay is immediate for $\gamma_{0}$ and $\gamma_{1}$ but is linear for the
spike counter, $\gamma_{2}$, which we explain below.

With this modeling framework, the TRANSFER and GATING neurons are
straightforward to model. The TRANSFER neuron has a voltage threshold just
below the product of spike current and synaptic weight,
$v_{\text{th}}=u_{s}w-1$. So, as soon as it receives an input spike, the
voltage jumps over threshold and emits a spike at the next tick. The GATING
neuron has a voltage threshold above the spike current-weight product,
$v_{\text{th}}=u_{s}w+1$, but below twice this value. Thus, it does not
respond to a single spike but emits a spike if it receives two (or more)
spikes. If a GATING neuron has only 2 dendrites, this behavior is equivalent
to a logical AND operation and hence we use it to gate the transfer of spiking
information dependent on whether a gating spike and an information spike
arrive simultaneously at the dendrites of the GATING neuron.

The COUNTER neuron, $\gamma_{2}$, uses current decay to drive the continuous
emission of spikes during a period of time equal to the current-weight product
of a set of input spikes to a set of dendrites. The current decay per tick
equals a single spike’s current. Therefore, the incoming current persists for
a time period equal to the current-weight product received at a given tick
then drops linearly to zero. This neuron works correctly (i.e. as a neuron to
count out a current-spike product as a set of spikes) if it does not receive
any spikes while its incoming current is decaying to zero. To model COUNTER
neuron behavior we set the current decay to zero ($\delta_{i}^{(u)}=0$) but
connect the neuron to itself with an inhibitory synapse. Therefore when it
emits a spike, it sends itself a negative current, equivalent to subtracting a
single spike’s current. In this way, the current decreases one spike current
per tick while emitting spikes until its current value drops to zero.

###  II.2 Two’s Complement Multiplication

In order to multiply two $M$-bit integers in two’s complement notation, we
first extend them to $N=2M$-bit numbers. By extend here, we mean duplicating
the highest bit to fill out the additional bits. For example, $1011\rightarrow
11111011$ or $0011\rightarrow 00000011$. Input to our binary multiplication
circuit is a vector of spikes (the multiplier), where existence of a spike
indicates a $1$ and non-existence of a spike indicates a $0$. The lowest bit
is input on an axon representing $2^{0}$ and the highest bit is input on an
axon representing $2^{N-1}$, and similarly with axons in between.

The vector of spikes propagates through a synaptic connectivity (i.e. a
matrix) that encodes a binary number (the multiplicand). This encoding is done
by shifting the binary representation (synapse of $1$ represents binary digit
of $1$ and synapse of $0$ represents binary digit of $0$) of the multiplicand
across the set of axonal synapses (see Fig. 2 for an encoding of the two’s
complement number $-1$). Shifts are performed with truncation at the highest
bit.

Our multiply subcircuit (core) uses COUNTER neurons. Therefore, when the
multiplier and multiplicand are multiplied, the multiply operation outputs a
sequence of at most $N$ sets of spikes representing a mixed binary-unary
representation.

The rest of the circuit is dedicated to sorting this mixed representation such
that pairs of numbers output by the multiply circuit may be carried and added,
then recombined, with subsequent carries and adds performed iteratively until
all carries are performed and only a binary representation remains.

###  II.3 Constituent Cores Making up the Full Circuit

Our full circuit is composed of six different types of core. The purpose of
these cores is 1) to turn the circuit on or off; 2) to control circuit timing
and to control how information is distributed across the circuit; 3) to
multiply input spikes representing the multiplier by synaptic weights
representing the multiplicand; 4) to spatially distribute the resulting
spikes; 5) to perform carry/add operations on pairs of spike registers; and 6)
to aggregate intermediate results and iterate the computation to completion.

These tasks are performed by on/off, clock, multiply, multiplexer, carry/add,
and resum cores (resp.), whose functions are outlined below (see Figs. 1-5).
We present example cores used for multiplying two $4$-bit numbers. The
generalization to an arbitrary number of bits is straightforward.

####  II.3.1 On/Off Core

The on/off core (Fig. 1) has two neurons responsible for turning the circuit
on or off. This is implemented using a connection to the clock core. The ON
neuron sends a spike to the first neuron in the clock core. It starts the
clock and the incoming spike starts circulating inside the clock core. A spike
from the OFF neuron inhibits all activity in the clock core and stops the
clock.

####  II.3.2 Clock Core

Figure 1: The On/Off and Clock Cores. Here, and in figures below, inputs
arrive from the top and outputs are emitted to the right. on and off inputs to
the clock core (left two columns) are indicated at top left. A cyclic circuit
(clock core) implementing an $N$-tick clocking cycle is the main set of axons
and neurons on the right. Recurrent connections implementing the cycle are
shown to the right. Output is additionally sent downstream to the multiplexer
circuit. (See text for description.)

The clock core (Fig. 1) has $N$ $\gamma_{0}$ neurons which indicate the state
of an $N$-tick clocking cycle. After the first neuron receives a spike, it
emits the spike to the multiplexer (below) and also sends the spike to the
second neuron. The second neuron emits the spike to the multiplexer and sends
it to the third neuron and so on. The last neuron sends the spike to the first
neuron after emitting it to the multiplexer. The emitted spikes communicate
the clock state to the multiplexer core, and serve as gating spikes to
regulate the multiplexer core’s function.

####  II.3.3 Multiply Core

The multiply core is designed to calculate the partial products of the
multiplicand, that are encoded in synaptic weights, by each bit of the
multiplier, which are encoded in spikes input to the multiply core. Basically,
this operation mimics the standard two’s compliment multiplication algorithm
as it would be performed manually, where the partial products of the
multiplier and multiplicand are calculated and then added together. In order
to add multiple binary partial products we use a combination of carry/add,
multiplexer, and resum cores. We use $\gamma_{2}$ neurons for the multiply
core, in which each neuron $i$ will emit $k_{i}$ spikes during $k_{i}$ ticks,
where $k_{i}$ is the summed current-weight input to neuron $i$.

Figure 2 shows the positive synapses that represent the number $-1$ in $8$-bit
two’s complement notation (i.e. $11111111$). For each subsequent spiking input
from the multiplier to the multiply core the binary representation of the
multiplicand is shifted to accommodate for the different positions of the bits
in the multiplier. For example, if the multiplier is binary $11111111$, then
the spikes come from all input channels. E.g. the first neuron receives a
single spike at time $t_{0}$ and emits a spike at time $t_{0}+1$. In turn the
second neuron receives two spikes at time $t_{0}$ and emits spikes during two
ticks from $t_{0}+1$ to $t_{0}+2$. In general if the multiplier and the
multiplicand are $a$ and $b$ respectively, then the spike from channel $i$
(which means that $a$ has $1$ at the $i$th position in two’s complement
notation) leads to the output that encodes $2^{i}\cdot b$ which is $b$,
shifted by $i$ positions.

Figure 2: The Multiply Core. The two’s complement representation of
$-1=11111111$ is shown encoded in the synapse in this example. This
representation is shifted and truncated across the set of inputs. All neurons
are COUNTER neurons in this core. See text for further details.

For this 4-bit implementation, the shifted multiplicand is truncated at
8-bits, since carries beyond this bit are unnecessary.

####  II.3.4 Multiplexer Core

The role of the multiplexer is to provide the correct input for the carry-add
cores from multiply and resum cores based on timed gating spikes from the
clock. It has $3N^{2}$ neurons that consist of $3N$ groups of $N$ neurons.
Each group fires only at a corresponding clock state. This is done in order to
synchronize the input from different carry-adds to the resum core. Since the
input is coming from the multiply and the resum cores and they use spike
counter neurons, which emit a sequence of spikes during a certain period of
time after receiving the input, the clock and synchronization are needed in
order for all outputs to combine together at a single tick.

The input for the multiplexer comes either from the multiply core or the resum
core but not both simultaneously since the multiply core shuts down after the
first $N$ ticks after receiving input and the resum core does not process any
spikes before that. Both the multiply and the resum cores can keep firing for
at most $N$ ticks after receiving a single input as this is the maximum number
of spikes that a single neuron can receive. Each of $N$ groups in the
multiplexer receives spiking input from the multiply and resum cores but
receives different gating spikes from the clock depending on the state in its
cycle. The gating spikes will only allow through-propagation of spikes that
came from the input core in the $j$th group, where $j\in 1,\dots,N$ is the
clock state.

Figure 3: The Multiplexer Core. Gating spikes from the clock (left register)
enable through-propagation of spikes arriving either from the multiply core or
from the resum core.

####  II.3.5 Carry/Add Core

The carry/add core receives two binary numbers simultaneously from the
multiplexer as an input (one of the inputs is delayed to achieve simultaneity,
see Fig. 6) and processes them using $\gamma_{2}$ spike counter neurons. Since
the multiplexer outputs $N$ numbers, we need $N/2$ carry-add cores to add them
in pairs. As output the carry/add core emits (at most) two sets of spikes
during two ticks after performing a single-bit carry-add operation. These
spikes are transferred to the resum core, which combines them all together for
further processing.

Figure 4: The Carry/Add Core. This core combines input from two summands and,
for two coincident spikes, performs a carry and adds the spikes. Note that
there is an assumption that if only one of the summands has a spike on a
particular input, it is summand 1. This assumption may be made since spike
sequences from the multiply and resum cores consist of a sequence of all ones,
followed by zeros. I.e. there are no zeros interspersed in the spiking
sequence.

This core operates under the assumption that two binary numbers $a$ and $b$
are coming as a sequential output of a spike counter neuron. This means that
if, for a certain position $i$, $a_{i}$ is zero, then $b_{i}$ is zero too. If
the $\gamma_{2}$ neuron did not emit a spike at time $t$, it will not emit it
at any further time unless it receives a new input.

The summation scheme is shown on Figure 4. If, for a certain position,
$b_{i}=1$ (then $a_{i}=1$, too), then the corresponding neuron receives both
positive and negative current values which cancel each other. If, in addition,
the previous bit had two ones ($b_{i-1}=1$), then that bit carries the value
of one to $i$th position and adds a positive current value there. Thus, only
if $a_{i-1}=1$, $b_{i-1}=1$, and $a_{i}=1$, but $b_{i}=0$, will neuron $i$
fire two spikes, otherwise it will fire zero or one spike.

####  II.3.6 Resum Core

The purpose of the resum core is to receive the spikes from all $N/2$ carry-
add cores at the same time and output the total number of spikes received by
each bit. This is modeled with $\gamma_{2}$ neurons. Since each carry-add
operation outputs not more than two spikes per neuron and there are $N/2$
carry-adds, the maximum number of spikes per neuron that the resum cores
receives and outputs is $N$. Here the difference from other cores that use
$\gamma_{2}$ neurons is that not all spikes come at the same tick: the input
is coming from carry-add cores which emit spikes during a period of two ticks.
However, each neuron still accumulates the total number of spikes and outputs
them correctly since the number of spikes coming to each neuron at the second
tick cannot exceed the number of spikes that came at the first tick due to the
nature of $\gamma_{2}$ output from the carry-add cores.

Figure 5: The Resum Core. Consisting of COUNTER neurons, this core adds input
spikes and outputs a new mixed unary-binary representation that is sent to the
multiplexer for a new carry-add iteration.

To keep resumming numbers and carrying the overflow bits, the resum core sends
the output back into the multiplexer which again sends pairs of numbers into
carry-add cores at appropriate clock ticks.

###  II.4 The Full Binary Matrix Multiplication Circuit

Figure 6: The Full Circuit. All cores are shown configured for two’s
complement multiplication. Note the delays after output from the multiplexer.
This coordinates the multiplexer output such that spikes from pairs of
registers arrive simultaneously at the carry-add cores. See text for further
details.

A binary multiplication circuit may be constructed by combining the cores
described above as shown in Fig. 6. This circuit may be extended to
accommodate multiple inputs representing the multiplication of synaptically
encoded matrices with spike-encoded input vectors (see Binary Vector-Matrix
Multiplication, below).

The maximum number of clock cycles required for all bits to carry over is
$N-1$ with each cycle having $N$ ticks. The total number of ticks required for
the multiplication of two binary numbers is $T=\log X(\log X-1)+4$ where $X$
is the resolution of the two’s complement numbers. The total number of neurons
used in this circuit is $3\log X+1.5(\log X)^{2}$ and the number of synapses
is $2\log X+2(\log X)^{2}+4.5(\log X)^{3}$.

###  II.5 Example Output

Figure 7: An example of the two’s complement binary multiplication of
$-1\times-1$. Red spikes from the multiply core and blue from the resum core
are shown. Successive carry-add operations transform the initial mixed unary-
binary encoding into a pure two’s complement representation (in this case the
number $00000001$.

In Fig. 7 we illustrate how the multiply core computes the partial sums that
are later summed using carry-add and resum cores for the product $-1\times-1$.
This example is chosen because it has the most nonzero bits in two’s
complement representation which is $-1\equiv 11111111$.

The vertical axis represents $8$ bits in multiply and resum cores and the
horizontal axis represents time. The multiply core emits spikes during $8$
ticks which are shown in red. These sequences of spikes are added in pairs
using multiplexer and carry-add cores and sent to the resum core. The work of
the resum core is shown as blue spikes, where the vertical lines indicate the
resum cycles. The resum cycle is shifted compared to the clock since it takes
$2$ ticks for spikes to go from the clock to resum core.

All the bits beyond the first one are summed up until eventually the summation
in the $8$th bit leads to an overflow and disappearance of spikes. The result
is shown at the beginning of the last resum cycle which represents $2^{0}=1$.

###  II.6 Binary Vector-Matrix Multiplication

Figure 8: Multiply Core for Inner Product. This example circuit shows the
two’s complement encoding of the vector $(-1,-1,-1,-1)$. To compute an inner
product, the elements of the other vector in the inner product are input in
the registers indicated as element $i$.

The two’s complement binary multiplication circuit described above can
straightforwardly be generalized to multiply $P\times Q$ matrices encoded in
synapses with $Q\times 1$ element vectors encoded in spikes.

To see this, consider the example multiply circuit shown in Fig. 8. Here, each
of the elements of a $4\times 1$ input vector are first encoded in spikes.
They are then input to a multiply core that has been extended to compute the
inner product of one row of the matrix with the input vector.

This results in a similar spiking output to that of the above circuit. The
main difference is that many more spikes can potentially be output due to the
multiple terms in the inner product. Nonetheless, with the same cores
described above, the full computation can be performed for all rows of the
matrix of interest, resulting in a binary representation of the vector-matrix
computation output.

A naïve implementation of the above concept results in a circuit with a total
number of ticks of $T=Q\log X(Q\log X-1)+4$, total neurons of $3Q\log
X+1.5(Q\log X)^{2}$, and total synapses of $2Q\log X+2(Q\log X)^{2}+4.5(Q\log
X)^{3}$. However, by incorporating a second set of clock, multiplexer, carry-
add, and resum cores, the time may be reduced at the expense of using more
neurons and synapses.

To do this, the first clock performs a first carry-add cycle, but then
triggers the second clock, which has a cycle length half that of the first
clock. Spikes output from the first resum core are propagated to the second
multiplexer, which only takes inputs for half the time of the first clock (as
the initial carry-adds decrease the number of spikes considerably) and outputs
half the number of carry-add pairs as the first multiplexer. This
significantly decreases the time needed for the full multiplication procedure
to complete.

More complex carry-add cycle management could further reduce the time needed
for a matrix-vector multiplication procedure, but we have not investigated
further refinements in detail.

With this matrix multiplication circuit in hand, one can then integrate
discretized partial differential equations, deep nets, or other linear
algebraic operations in the usual way. For instance, the diffusion equation
may be discretized to a set of ordinary differential equations. The time
evolution is then computed iteratively via matrix multiplication.

###  II.7 Random Walk Circuit

In this section, we present a circuit which computes a 1-dimensional random
walk that starts at a location $x_{0}$ and chooses to step $\delta_{\pm}=\pm
1$ with probability $p_{\pm}=0.5$. We again use two’s complement encoding for
N-bit numbers. For independent random walkers, the full multiplication circuit
is unnecessary but carry-add operations are still necessary. This allows for a
significant reduction in circuit size. We introduce new types of neurons that
we use in this circuit:

  * •

$\gamma_{3}$ – a CONSTANT neuron, which emits spikes every $M$th tick;

  * •

$\gamma_{4}$ – a RANDOM neuron, which emits spikes every tick with probability
$p=0.5$.

The behavior of the $\gamma_{3}$ neuron is straightforward and is modeled with
the help of bias current and refractory delay. The period $M$ is equal to the
length of a single summation cycle which is $N+1$. In order to model
randomness we use Loihi’s capability to add random noise to the voltage
threshold. This is implemented by setting noise to take only two values (0 and
2) and bias current to be equal to the voltage threshold. Thus the total
voltage exceeds the threshold with probability $0.5$ as needed.

The full circuit consists of two single neurons, $\gamma_{3}$ (constant) and
$\gamma_{4}$ (random), and three cores: a perturbation core of $\gamma_{1}$
neurons, a location core of $\gamma_{0}$ neurons, and a carry core of
$\gamma_{1}$ neurons. The circuit needs $N+1$ ticks after supplying the
perturbation value $\pm 1$ to the location core to perform the summation
operation. If an overflow happens in the first digit, it can need up to $N$
more operations to carry the overflow to more significant digits.

The diagrams in Figs. 9, 10, and 11 show the circuit elements. The synaptic
weights in the perturbation core combine the spikes coming from neurons
$\gamma_{3}$ and $\gamma_{4}$ in order to output the binary representation of
either $+1$ or $-1$ depending on the output of $\gamma_{4}$ (Figure 9). The
large red circle in the first neuron represents a double weight which is
needed to pass the larger threshold of $\gamma_{1}$.

Figure 9: The Perturbation Core. Here a spike is input at every iteration
(constant) and a spike is randomly input at every iteration. In concert with
the given synapses, an output of either $+1$ or $-1$ is output. The output
drives the downstream location and carry circuits in the random walker.

The location and carry cores receive the same input (starting position,
perturbation, carry, and location) but process it with different neurons:
$\gamma_{0}$ for the location core and $\gamma_{1}$ for the carry core
(Figures 10 and 11). The circuit works in such a way that at each tick these
two cores receive spikes from not more than 2 inputs (either starting position
and perturbation, location and perturbation, or location and carry). The sum
of location and carry represents the current location and it can take up to
$N$ ticks to carry the overflow digits after receiving the perturbation
spikes. After $N$ ticks the carry stops producing any spikes and the location
core holds the complete current location result.

The number of ticks required for this circuit to make $M$ steps in terms of
resolution is $T=2+M(\log X+1)$. The number of neurons used for a single agent
is $2+3\log X$ and the number of synapses is $2\log X+8(\log X)^{2}$.

Figure 10: The Location Core. A starting location (used only once for
initialization), perturbation from the perturbation core, and carry from the
carry core are input. The carry is output and also fed back to the location
register. Figure 11: The Carry Core. This core has the same structure as the
location core, but with location and carry registers reversed.

The full circuit is depicted on Figure 12.

Figure 12: The Random Walk Circuit. This simple circuit combines the
perturbation, location, and carry cores to perform a random walk. See text for
details.

##  III Results

Fig. 13 shows the random walk of a single agent in $8$-bit space. Fig. 14
shows diffusion by several agents in $8$-bit space. Fig. 15 shows the
distribution of the locations of $500$ agents at times $t=0,50,130$ in $5$-bit
space. Here the space was taken small to illustrate the effect of the periodic
boundary conditions which eventually lead from a Gaussian to uniform
distribution.

Figure 13: The random walk of a single agent in $8$-bit resolution as a
function of time. Figure 14: The random walks of several agents in $8$-bit
resolution as a function of time. Figure 15: The distributions of agents’
locations as the distribution evolves from a delta function, to a normal
distribution, to an asymptotically flat distribution.

In Fig. 16, we show spike propagation as a function of time in perturbation,
location, and carry cores superimposed. Here, three repeated perturbations of
$-1$ increment the location from an initial value of $0$, successively to
$-1$, then $-2$ (two’s complement). The action of the location and carry cores
is evident as the spikes are successively reduced to final two’s complement
numbers at each iteration.

Figure 16: Illustration of spikes in a random walk scheme that starts at
$x_{0}=1$ and makes three steps, each time receiving a perturbation of $-1$.
Spikes from all three cores are shown superposed.

The performance of the random walk circuit on Loihi is shown in Figure 17 for
different resolutions and numbers of agents. All runs contain $1000$ steps and
we expect the performance to be independent of the number of agents since
Loihi propagates all agents in parallel. We see that the time per tick is
almost independent of the resolution of the number which results in time per
random walk step being proportional to $\log X$ where $X$ is the resolution.

However, the plot shows an initial linear growth of time with the number of
agents. This indicates the counterintuitive result that only asymptotically is
good parallelization achieved. Without access to the full scheduling and
distribution of the circuit on chip it is difficult to pinpoint the reason for
the initial growth in iteration time versus number of agents.

Figure 17: The performance of the random walk (time per random walk step)
algorithm implemented on Loihi as a function of the number of agents
distributed on the chip. Note that, at a certain point, the number of $16$-bit
walkers exceeded the size of our Kapoho Bay system, so we were only able to
assess up to $2000$ $16$-bit walkers.

##  IV Discussion

We have constructed a neuromorphic circuit that performs a matrix-vector
product in two’s complement representation. The timing of the computation is
controlled by a clock, and a set of other cores is used to carry out the
reduction of the representation to a final, two’s complement encoded binary
integer.

We have shown how this circuit can be extended for application to matrix
multiplication and hence used for various applications, such as the
integration of linear equations.

Additionally, we introduced a simpler circuit, still using two’s complement
numbers to perform a random walk and shown that this circuit parallelizes well
and asymptotes to a constant time per step for large numbers of agents.

Our work here demonstrates that, at the cost of a larger neuronal footprint,
binary coded multiplication may be performed on circuits with sizes and run
times that are poly-logarithmic in the resolution of the numbers being
multiplied. This results in a significant savings in multiplication time and
spikes generated relative to unary coding, where number resolution is
proportional to spike number. This can result in significant reductions in
power consumed due to the relationship between spike number and power
consumption.

###### Acknowledgements.

Acknowledgements. We thank Alpha Renner for helpful conversations. We thank
Intel’s Loihi development team and the Intel Neuromorphic Research Community
for their support. This work was supported by the LANL ASC Beyond Moore’s Law
project (O.I.) and by the US Department of Energy National Nuclear Security
Administration’s Office of Defense Nuclear Nonproliferation Research &
Development (DNN R&D) at Los Alamos National Laboratory under contract
89233218CNA000001 (A.T.S.). LANL approval designation: LA-UR-21-22286

## References

  * Harris et al. [2012] Julia J Harris, Renaud Jolivet, and David Attwell.  Synaptic energy use and supply.  _Neuron_ , 75(5):762–777, 2012.
  * Levy and Baxter [2002] William B Levy and Robert A Baxter.  Energy-efficient neuronal computation via quantal synaptic failures.  _Journal of Neuroscience_ , 22(11):4746–4755, 2002.
  * [3] Mike Davies.  Exploring neuromorphic computing for AI: Why spikes? (part one).  _Intel, Artificial Intelligence, Monthly Newsletter_.
  * Diehl et al. [2016] Peter U. Diehl, Bruno U. Pedroni, Andrew Cassidy, Paul Merolla, Emre Neftci, and Guido Zarella.  TrueHappiness: Neuromorphic emotion recognition on TrueNorth.  _International Joint Conference on Neural Networks (IJCNN)_ , pages 1–10, 2016.
  * Severa et al. [2018] William Severa, Rich Lehoucq, Ojas Parekh, and James B Aimone.  Spiking neural algorithms for Markov process random walk.  In _2018 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–8. IEEE, 2018.
  * Parekh et al. [2017] Ojas Parekh, Cynthia A Phillips, Conrad D James, and James B Aimone.  Constant depth and subcubic size threshold circuits for matrix multiplication.(2017), 2017.
  * Reeder et al. [2019] Leah Evelyn Reeder, James Bradley Aimone, and William Mark Severa.  The future of computing: Integrating scientific computation on neuromorphic systems.  Technical report, Sandia National Lab.(SNL-NM), Albuquerque, NM (United States), 2019.
  * Davies et al. [2018] Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham N. Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios D. Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, Yuyun Liao, Chit-Kwan Lin, Andrew Lines, Ruokun Liu, Deepak Mathaikutty, Steven McCoy, Arnab Paul, Jonathan Tse, Guruguhanathan Venkataramanan, Yi-Hsin Weng, Andreas Wild, Yoonseok Yang, and Hong Wang.  Loihi: A neuromorphic manycore processor with on-chip learning.  _IEEE Micro_ , 38(1):82–99, 2018.  doi: 10.1109/MM.2018.112130359.  URL https://doi.org/10.1109/MM.2018.112130359.


# Initialization of quantum simulators by sympathetic cooling

Meghana Raghunandan<EMAIL_ADDRESS>Institut für
Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, 30167
Hannover, Germany  Fabian Wolf  QUEST Institut, Physikalisch-Technische
Bundesanstalt, Bundesallee 100, 38116 Braunschweig, Germany  Christian
Ospelkaus  QUEST Institut, Physikalisch-Technische Bundesanstalt, Bundesallee
100, 38116 Braunschweig, Germany  Institut für Quantenoptik, Leibniz
Universität Hannover, Welfengarten 1, 30167 Hannover, Germany  Piet O. Schmidt
QUEST Institut, Physikalisch-Technische Bundesanstalt, Bundesallee 100, 38116
Braunschweig, Germany  Institut für Quantenoptik, Leibniz Universität
Hannover, Welfengarten 1, 30167 Hannover, Germany  Hendrik Weimer  Institut
für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, 30167
Hannover, Germany

###### Abstract

Simulating computationally intractable many-body problems on a quantum
simulator holds great potential to deliver insights into physical, chemical,
and biological systems. While the implementation of Hamiltonian dynamics
within a quantum simulator has already been demonstrated in many experiments,
the problem of initialization of quantum simulators to a suitable quantum
state has hitherto remained mostly unsolved. Here, we show that already a
single dissipatively driven auxiliary particle can efficiently prepare the
quantum simulator in a low-energy state of largely arbitrary Hamiltonians. We
demonstrate the scalability of our approach and show that it is robust against
unwanted sources of decoherence. While our initialization protocol is largely
independent of the physical realization of the simulation device, we provide
an implementation example for a trapped ion quantum simulator.

##  I Introduction

Quantum simulation is an emergent technology that can potentially solve
important open problems related to high-temperature superconductivity,
interacting quantum field theories, or many-body localization Georgescu2014 .
While a series of experiments demonstrated the successful implementation of
Hamiltonian dynamics within a quantum simulator Greiner2002 ; Schneider2008 ;
Jordens2008 ; Friedenauer2008 ; Lanyon2011 ; Aidelsburger2013 ; Miyake2013 ;
Alvarez2015 ; Mazurenko2017 ; Bernien2017 ; Zhang2017 ; Guardado-Sanchez2018 ;
Lienhard2018 , these works had the simulator initialized in an easily
accessible state such as a product state. Consequently, adiabatic evolution
from an initial Hamiltonian whose ground state can be prepared to the final
Hamiltonian of interest has been used. However, this approach becomes
challenging across quantum phase transitions, especially if the transition is
of first order.

Our strategy to overcome this problem builds on the recent advances in using
dissipative quantum systems to engineer interesting many-body states as the
attractor states of such an open quantum many-body system Diehl2008 ;
Verstraete2009 ; Weimer2010 ; Barreiro2011 ; Carr2013a ; Rao2013 ; Lin2013 ;
Shankar2013 ; Cormick2013 ; Morigi2015 ; Neto2017 . In the past, these
dissipative state engineering schemes have been limited to ground states of
stabilizer or frustration-free Hamiltonians Verstraete2009 ; Weimer2010 ;
Weimer2011 ; Roghani2018 , whose ground state can be found by performing local
optimizations alone. Unfortunately, almost all many-body Hamiltonians of
interest lie outside this class, requiring generalization of the dissipative
state preparation procedure.

Figure 1: Sympathetic cooling of a quantum simulator. (a) A system of $N$
spins performing the quantum simulation is interacting with an additional bath
spin that is dissipatively driven. (b) Sketch of the energy level structure
showing resonant energy transport between the system and the bath, after which
the bath spin is dissipatively pumped into its ground state. (c) Level scheme
for the implementation with trapped 40Ca+ ions.

Here, we present a previously unexplored paradigm for the dissipative
initialization of a quantum simulator. We consider a coupling of the many-body
system performing the quantum simulation to an auxiliary particle that is
dissipatively driven. Crucially, the energy splitting within the auxiliary
particle is chosen such that it becomes resonant with the many-body excitation
gap of the system of interest, i.e., the difference of the ground-state energy
and the energy of the first excited state. Under such a resonance condition,
the energy of the quantum simulator is efficiently transferred to the
auxiliary particle such that the former is cooled sympathetically Cormick2013
; Neto2017 . Although this setup is only resonant at a single energy, the
density of states increases exponentially with energy, resulting in the
lowest-lying excitations being the bottleneck for fast ground-state
preparation, see the Supplementary Materials for details. While the value of
the many-body excitation gap is usually unknown before performing the
simulation, we demonstrate that the gap can actually be determined from the
quantum simulation data in a spectroscopic measurement. Hence, the dissipative
initialization process provides important information about the many-body
system of interest at the same time. Notably, we show that the cooling by a
single auxiliary particle is efficient, and it is especially robust against
unwanted noise processes occurring in the quantum simulator.

(a) |  |  (b) |
---|---|---|---
Figure 2: Sympathetic cooling of the transverse field Ising model in the
ferromagnetic phase ($J/g=5$, $N=5$, $f_{x,y,z}=\\{1,1.1,0.9\\}$). The speed
of the cooling dynamics as well as the final energy of the system depend on
the system-bath coupling $g_{sb}$ for $\gamma/g=1.9$ (a) and the dissipation
rate $\gamma$ for $g_{sb}/g=1.15$ (b). The ground state energy is indicated by
the dashed line. The insets show that the ground state can be prepared with
greater than 90% fidelity.

To be explicit, we consider different paradigmatic one-dimensional (1D) spin
$1/2$ many-body systems coupled to a single dissipatively driven auxiliary
bath spin (see Fig. 1). This setup can be readily generalized to bosonic or
fermionic many-body systems with a larger local Hilbert space, to settings
incorporating several bath particles, and to higher spatial dimensions. In the
following, we assume a 1D chain of $N$ spins governed by the Hamiltonian
$H_{sys}$. One boundary spin of the system is coupled to the auxiliary bath
spin via an interaction Hamiltonian of the form
$H_{int}=g_{sb}\sum\limits_{x,y,z}f_{i}\sigma_{i}^{(N)}\sigma_{i}^{(b)}$,
where $g_{sb}$ is the strength of the system-bath interaction and the
$\sigma_{i}$ refer to Pauli matrices. The exact values of the dimensionless
parameters $f_{i}$ are not particularly important. In the models studied here,
we find that it is either favorable to choose them roughly equal or have one
dominant contribution. In addition, to avoid any symmetries in the interaction
preventing the cooling of certain degrees of freedom, it is beneficial to
assign slightly different values to them.

The Hamiltonian of the bath spin $H_{bath}$ is given by
$H_{b}=(\Delta/2)\sigma_{z}^{(b)}$. The dissipation channel acting on the bath
spins performs dissipative spin flips from the up spin state to the down spin
state occurring with a rate $\gamma$. Then, the total dynamics is described by
a quantum master equation in Lindblad form

$\frac{d}{dt}\rho=-\frac{i}{\hbar}[H,\rho]+\gamma\left(\sigma_{-}^{(b)}\rho\sigma_{+}^{(b)}-\frac{1}{2}\left\\{\sigma_{+}^{(b)}\sigma_{-}^{(b)},\rho\right\\}\right),$
(1)

where $H=H_{sys}+H_{bath}+H_{int}$ is the total Hamiltonian of the $N+1$ spin
system Breuer2002 .

We would like to stress that such a setup imposes only modest requirements for
an experimental implementation, which works equally well for both analog and
digital quantum simulators. In particular, we note that our setup does not
require control over individual particles of the quantum simulator. In our
case, it is sufficient to merely be able to control the bath particle
independently of the rest of the system. In addition, the dissipative dynamics
can be induced by measuring the spin state of the bath spin followed by a spin
flip conditional on measuring the spin in the up state.

In the methods section, we give a detailed implementation guide for a trapped
ion quantum simulator.

##  II Results

###  II.1 Ising chain in a transverse field

As the first paradigmatic model, we consider the Ising model in a transverse
field, given by the Hamiltonian

$H_{sys}=g\sum\limits_{i=1}^{N}\sigma_{z}^{(i)}-J\sum\limits_{i=1}^{N-1}\sigma_{x}^{(i)}\otimes\sigma_{x}^{(i+1)},$
(2)

where $g$ is the strength of the transverse field, and $J$ is the coupling
constant for the Ising interaction. As the Pauli matrices do not commute with
each other, it is impossible to minimize the interaction terms and the
magnetic field term at the same time, meaning that already this simple model
lies outside of the class of frustration-free Hamiltonians. The transverse
field Ising model is known to undergo a quantum phase transition at $g=J$ from
a paramagnetic phase ($g>J$) to a ferromagnetic phase ($g<J$) Sachdev1999 . In
the following, we will set the energy splitting of the bath spin $\Delta$ to
be identical to the many-body gap $\Delta E=E_{1}-E_{0}$ of the transverse
field Ising model, where $E_{0}$ ($E_{1}$) is the energy of the ground state
(first excited state). In the ferromagnetic phase, the ground state becomes
doubly degenerate for large system sizes. Because we are not interested in
cooling into a particular ground state, $E_{1}$ refers to the first excited
state above the ground state manifold. Below, we will demonstrate that
choosing the bath spin splitting as $\Delta=\Delta E$ leads to optimal
cooling, and we will show how to extract the (a priori unknown) energy gap
$\Delta E$ from the quantum simulation results.

Let us now analyze the cooling performance of the setup by tracking the system
energy $\langle H_{sys}\rangle$ of the transverse field Ising model in wave-
function Monte Carlo simulations of $N=5$ spins, initially in the
experimentally accessible state of all spins pointing up. Figure 2 shows that
the energy of the system decreases rapidly and finally approaches a value that
is close to the numerically calculated ground-state energy. The cooling
performance depends on the choice of the system-bath coupling $g_{sb}$ and the
dissipation rate $\gamma$. In the following, we assume that the time available
for the cooling remains fixed. Then, if $g_{sb}$ is too small, the cooling
dynamics is very slow. On the other hand, if $g_{sb}$ is too large, then the
system and the bath spin will become strongly entangled, and the cooling
performance is reduced. Similarly, if $\gamma$ is too small, then the cooling
is slowed down in the same way, while a too large value of $\gamma$ will lead
to a quantum Zeno suppression of the energy transfer required for the cooling
process. Hence, there should be an optimal choice for $g_{sb}$ and $\gamma$,
which leads to a minimum in energy within the available time.

To find this optimal choice, we use a model-independent quantity to measure
the cooling performance. For this, we calculate the fidelity of the state of
the system with respect to the ground-state manifold of the transverse field
Ising model. The fidelity $f$ is given by

$f=\langle\Pi_{g}\rangle=\text{Tr}\left\\{\rho(t)\Pi_{g}\right\\},$ (3)

where $\Pi_{g}=\sum\limits_{i}|\psi_{0}^{i}\rangle\langle\psi_{0}^{i}|$ is the
sum of the projectors onto the ground states Nielsen2000 . As the inset of
Fig. 2a and 2b shows, the ground state can be prepared with more than 90%
fidelity for the optimal choice of $g_{sb}=1.15\,g$ and $\gamma=1.9\,g$.

We can also relate the fidelity $f$ to the system energy $\langle
H_{sys}\rangle$. For this, we introduce a dimensionless excitation energy
$\epsilon$, measured in units of the many-body gap $\Delta E$, i.e.

$\epsilon=\frac{\langle H_{sys}\rangle-E_{0}}{\Delta E}.$ (4)

In the low-energy limit $\epsilon\ll 1$ and assuming that the excitation
energy is mostly concentrated in low-energy excitations, $\epsilon$ is related
to the fidelity according to $\epsilon=1-f$.

We have also checked that our cooling procedure works independently of the
choice of $J/g$, i.e., both in the ferromagnetic phase and in the paramagnet,
as well as independently of the initial state (see the Supplementary Materials
for details.) Even in the critical regime ($J/g\sim 1$), where the many-body
gap is closing, we observe a similar cooling performance. To substantiate this
point, and also to demonstrate that our cooling protocol is not limited to a
particular model, we turn to the especially challenging case of a critical
Heisenberg chain in the following section.

###  II.2 Antiferromagnetic Heisenberg model

As a second paradigmatic quantum many-body model, we investigate the
antiferromagnetic Heisenberg chain, given by the system Hamiltonian

$H_{sys}=J\sum\limits_{i=1}^{N-1}\sum\limits_{j=x,y,z}\sigma_{j}^{(i)}\otimes\sigma_{j}^{(i+1)}.$
(5)

This model exhibits an $SU(2)$ symmetry and serves as the critical point of a
Kosterlitz-Thouless transition when the strength of the $\sigma_{z}\sigma_{z}$
interaction is varied Schollwock2004 . As the many-body gap vanishes in the
thermodynamic limit, this model represents a particularly challenging case for
our cooling protocol. In addition, the ground state at the critical point is
highly entangled Latorre2004 ; hence, we also test the capability of our
cooling protocol to prepare entangled quantum many-body states.

(a) |  |  (b) |
---|---|---|---
Figure 3: Sympathetic cooling of the antiferromagnetic Heisenberg model
($N=4$, $g_{sb}/J=0.2$, $\gamma/J=0.6$, $f_{x,y,z}=\\{0.4,2.3,0.3\\}$). (a)
The efficiency of the cooling procedure depends on the choice of the bath spin
splitting $\Delta$. (b) The optimal cooling leading to the lowest system
energy $\langle H_{sys}\rangle$ corresponds to setting $\Delta$ to the many-
body gap $\Delta E$ (vertical dashed line). The same minimum is observed when
measuring the energy $E_{dis}$ that is being dissipated during the cooling
process. The ground state energy is indicated by the horizontal dashed line.

The antiferromagnetic Heisenberg model adds one minor complication concerning
its ground state preparation compared with the Ising model. Because of an
approximate symmetry conserving certain spin-wave excitations, the ground
state cooling performance is limited when the system-bath coupling is
restricted to the last spin of the chain. We resolve this issue by an
additional system-bath coupling of strength $g_{sb}/2$ to the second last spin
of the chain. Fig. 3 shows the cooling performance in terms of the system
energy $\langle H_{sys}\rangle$, with an initial state of spins pointing up
and down alternately, as a function of the splitting of the bath spin
$\Delta$. As in the case of the transverse field Ising model, $\langle
H_{sys}\rangle$ decreases rapidly and reaches a final value that is close to
the ground state energy $E_{0}$. In addition, the cooling is optimal when
$\Delta$ is chosen to be identical to the many-body gap $\Delta E$ ($f=0.9$).
Hence, experimentally measuring $H_{sys}$ as a function of $\Delta$ allows one
to obtain the value of the many-body gap $\Delta E$, which in itself is an
important quantity to understand a quantum many-body system. We also find the
final state to be highly entangled (see the Supplementary Materials for
details).

However, on many quantum simulation architectures, it might be difficult to
experimentally measure the system energy $H_{sys}$, as this will typically
require one to perform tomography on all the operators that appear in the
system Hamiltonian. Further challenges arise in architectures where not all
coupling constants in the Hamiltonian can be perfectly controlled, leading to
additional uncertainties in the estimated value of $\Delta E$.

Fortunately, it is possible to obtain $\Delta E$ by measuring only the bath
spin. The key idea is to measure the energy $E_{dis}$ that is dissipated
during the cooling dynamics. Crucially, this energy is related to the number
of quantum jumps $N_{jump}$ by the relation $E_{dis}=N_{jump}\Delta$, as a
quantum jump will lower the energy of the bath spin by $\Delta$. We note that
there are two different ways to obtain $N_{jump}$. First, one can directly
count the number of quantum jumps, e.g., by counting the number of emitted
photons, if the dissipative flip of the bath spin is realized by a spontaneous
emission event. In many setups, however, collecting each emitted photon with
high probability might be too challenging. However, as a second method, one
can also obtain $N_{jump}$ via the integrated probability to find the bath
spin in the up state according to

$N_{jump}=\gamma\int\limits_{0}^{t_{p}}\text{Tr}\left\\{\sigma_{+}^{(b)}\sigma_{-}^{(b)}\rho(t)\right\\}\,dt$
(6)

where $t_{p}$ is the total preparation time. As shown in Fig. 3, the minimum
of $E_{dis}$ is almost identical to the minimum in $H_{sys}$, corresponding to
the case where the splitting of the bath spin $\Delta$ is identical to the
many-body gap $\Delta E$. We note that if the system-bath coupling $g_{sb}$ or
the dissipation rate $\gamma$ is chosen too large, then the difference between
the minima in $\langle H_{sys}\rangle$ and $E_{dis}$ becomes appreciably
larger. We also observe that $E_{dis}$ is slightly larger in magnitude than
the system energy; this can be attributed to the fact that even in the limit
of large times, a finite probability for quantum jumps remains as the ground
state of the system Hamiltonian is not a perfect dark state of the quantum
master equation Lemeshko2013a  due to the finite system-bath coupling
$g_{sb}$. These additional jumps can also happen for non optimal values of
$\Delta$, leading to a broadening of the dissipated energy $E_{dis}$ in Fig.
3b compared with the system energy $\langle H_{sys}\rangle$.

###  II.3 Efficiency of the cooling protocol

Figure 4: Scalability of the cooling protocol. The preparation time $t_{p}$ to
reach a final dimensionless energy of $\epsilon=0.2$ grows linearly on a log-
log scale, i.e., $t_{p}\sim N^{\alpha}$. The solid line is a fit to the data
according to $\alpha={\color[rgb]{0,0,0}3.1\pm 0.1}$.

For any quantum state preparation protocol, it is crucial to determine how its
properties behave when the size of the system is increased. A protocol is
called efficient when the resources required (i.e., the preparation time) grow
at most polynomially with the system size. To determine the scaling with
system size in an unbiased way, we compute the preparation time $t_{p}$ that
is required to cool the system down to a fixed dimensionless energy
$\epsilon$, while the system bath coupling $g_{sb}$ and the dissipation rate
$\gamma$ are chosen such that the cooling is optimal. Within our numerical
simulations, we use a standard nonlinear optimization scheme (see the Methods
section for details). In an actual quantum simulator, one can use a hybrid
algorithm in which the energies measured on the quantum device are fed back
into the classical optimization algorithm Kokail2019 .

Figure 4 shows the scaling behavior of $t_{p}$ for the transverse field Ising
model. Although the system is cooled across the phase transition into the
ferromagnet, the preparation time grows only polynomially with the system
size. The same scaling is also observed for the antiferromagnetic Heisenberg
model (see the Supplementary Materials for details). This scaling behavior
underlines that our cooling procedure is already scalable when using only a
single bath spin. As the number of particles is often a scarce resource in a
quantum simulator, the required minimal overhead for the initialization allows
us to use almost all of the particles for the actual quantum simulation.

###  II.4 Performance under decoherence

So far, the only source of decoherence in our considerations stems from the
dissipative flips of the bath spin. However, in most quantum simulation
architectures, there will also be unwanted decoherence processes in the system
performing the quantum simulation. Therefore, it is crucial to determine the
consequences of this additional decoherence on the performance of our cooling
protocol.

As an additional source of decoherence, we consider $\sigma_{z}$ spin flips in
the quantum simulation of the transverse field Ising model, applied with a
rate $\kappa$ to all $N$ spins of the quantum simulator. In the ferromagnetic
phase, such a spin flip will create two neighboring domain-wall excitations,
i.e., when applied to the ground state, the dimensionless energy will
approximately increase to $\varepsilon\approx 2$. This type of decoherence
represents a worst case scenario of all local decoherence processes. Hence, we
expect that this scenario is quite generic and that our findings should also
apply to other many-body models.

To analyze the consequences of these additional decoherence channels, we
consider the quantity $\kappa t_{p}$, which is essentially the probability of
any spin to undergo a decoherence event during the preparation time. Then,
tracking how the energy $\epsilon$ behaves as a function of $\kappa t_{p}$
allows us to assess the robustness of our cooling protocol under additional
decoherence.

Figure 5 shows the system energy for different decoherence rates, from which
the behavior of $\epsilon$ is calculated. Crucially, we find that the system
contains one excitation, $\epsilon\approx 1$ at a value of $\kappa
t_{p}\approx 2$. This means that the system picks up one excitation when on
average all the spins have undergone a decoherence event. This is in stark
contrast to the scaling observed in adiabatic state preparation protocols,
where the error probability is typically given by the probability that a
single spin undergoes a decoherence event, i.e., proportional to $N\kappa
t_{p}$ Weimer2012 . This improved robustness against decoherence can be
attributed to the fact that our state preparation protocol itself is
dissipative and therefore can self-correct decoherence events.

Figure 5: Cooling performance in the presence of decoherence in the quantum
simulator for the transverse field Ising chain ($J/g=5$, $N=4$). The inset
shows the dimensionless energy $\epsilon$ as a function of the product $\kappa
t_{p}$, where $t_{p}$ was taken from the dynamics without decoherence
corresponding to a ground state preparation fidelity of $f=0.9$ (dashed line).

###  II.5 Experimental realization

The proposed initialization protocol can be implemented in a trapped ion
system with state-of-the-art technology, e.g. by confining a 1D ion string in
a linear Paul trap. Here, we propose an implementation with 40Ca+-ions in a
setup similar to the one described in Jurcevic2014 . The spin states are
encoded in the optical qubit, $|\\!\downarrow\rangle=|S_{1/2},m=+1/2\rangle$
and $|\\!\uparrow\rangle=|D_{5/2},m=+5/2\rangle$ (see Fig. 1c) with an energy
splitting of $\hbar\omega_{0}$, coherently manipulated by radial laser beams.
E.g. the rightmost ion serves as the bath spin (index $b$), while its laser-
induced coupling to the neighboring ion (index $s$) implements the system-bath
coupling. The bath ion can be isolated from the system interaction by shelving
the population to an auxiliary state
$|\text{aux}\rangle_{b}=|D_{5/2},m=-5/2\rangle_{b}$ with a laser beam
addressing only the bath ion. An experimental realization requires the
implementation of the system and system-bath Hamiltonians. For simplicity, we
suggest to implement $H_{sys}$ and $H_{sb}$ in an interleaved fashion by
trotterizing the total interaction Lloyd1996 ; Lanyon2011 .

In trapped ion systems, $H_{sys}$ for the transverse field Ising model
Friedenauer2008  has been realized with up to 53 qubits Zhang2017 . For this
purpose, a global bichromatic laser beam with frequency $\omega_{0}\pm\delta$
implements a gate operation by coupling to all radial modes. If $\delta$ is
larger than the center-of-mass mode frequency, then the resulting spin-spin
coupling coefficient shows a power law scaling $J_{i,j}\propto
1/|i-j|^{\alpha}$ Kim2009 , where $\alpha$ can be varied between $0$ and $3$
by changing the radial confinement. Implementation of the Heisenberg model is
possible by interleaving the spin-spin coupling gates with single-qubit
rotations performing a basis change from $\sigma_{x}$ to $\sigma_{y}$ and
$\sigma_{z}$.

We propose to implement $H_{sb}$ with a separate laser that provides single
ion addressing for the bath spin and the neighboring system spin. A Mølmer-
Sørensen gate Benhelm2008 ; Roos2008  on the radial motional modes bridges two
different energy gaps, $\omega_{s}$ and $\omega_{0}$, similar to a two-species
gate Tan2015 , and provides a $\sigma^{(N)}_{x}\sigma^{(b)}_{x}$-type coupling
of the spins. For the bath spin, the laser frequencies will be
$\omega_{0}\pm\delta$ and for the system spin will be $\omega_{s}\pm\delta$
with $\omega_{s}=\Delta E/\hbar$ for optimal cooling. Tuning the latter
frequency corresponds to searching for the resonance condition described in
the main text. Again, $\sigma^{(N)}_{x}\sigma^{(b)}_{x}$-gates interleaved
with single qubit rotations on both ions implement
$\sigma^{(N)}_{x}\sigma^{(b)}_{x}$, $\sigma^{(N)}_{y}\sigma^{(b)}_{y}$, and
$\sigma^{(N)}_{z}\sigma^{(b)}_{z}$. The coupling between the bath spin and the
second last spin for the Heisenberg model can be realized by extending the
addressing laser to the second last spin such that the power law of the
system-bath interaction has an exponent $\alpha_{sb}\sim 1$. This comes at the
expense of an additional interaction between the last and the second-last
spins of the system, which is significantly weaker than $J$ and can therefore
be neglected. This additional coupling may also be canceled using an
additional addressing laser.

Assuming $\Delta E$ is already known, repumping from $|\\!\uparrow\rangle_{b}$
to P3/2 and a subsequent spontaneous decay to $|\\!\downarrow\rangle_{b}$ on
the bath ion can be used to provide a channel for dissipation. The strength of
dissipation, $\gamma$, within the trotterized scheme can be adjusted by the
repumping laser intensity, i.e. the repumping probability during each Trotter
cycle. For determination of $\Delta E$ by recording $N_{jump}$, every
scattered photon during the repump process has to be detected. This is
accomplished by an electron shelving scheme in which the population in
$|\\!\downarrow\rangle_{b}$ is hidden in state $|\text{aux}\rangle_{b}$ and a
potentially scattered photon bringing the bath ion from
$|\\!\uparrow\rangle_{b}$ to $|\\!\downarrow\rangle_{b}$ is detected by
measuring fluorescence on the $|\\!\downarrow\rangle_{b}$ (S1/2) to P1/2
transition. To avoid a perturbation of the system spins, the detection laser
has to be tightly focused onto the bath ion.

Figure 6: Cooling performance of an Ising-like chain of $5+1$ ions of
$t_{p}=80\hbar/g=24$ ms. The blue line shows the dynamics in the decoherence-
free case resulting in a fidelity of $f=0.92$, while the orange line indicates
the dynamics under a collective decoherence mechamism with rate
$\kappa_{c}=3.3$ Hz, resulting in $f=0.89$. The dashed line indicates the
ground state energy of the system.

To be more specific, we assume 5 40Ca+ ions in a linear chain with single ion
axial and radial trapping frequencies of $\omega_{z}=2\pi\times 0.15$ MHz and
$\omega_{r}=2\pi\times 0.5$ MHz, respectively Jurcevic2014 . With a resonant
Rabi frequency of $2\pi\times 125$ kHz for all ions, $J_{i,j}$ ranges between
$2\pi\times 2.7$ kHz and $2\pi\times 1.2$ kHz for a detuning of
$\delta-\omega_{r}\approx 2\pi\times 10.5$ kHz, while the largest Lamb-Dicke
parameter is given by $\eta_{\textrm{max}}=0.128$. For these parameters, the
spacing between the bath spin and the nearest system spin of around
$14~{}\upmu\text{m}$ is sufficiently large to provide a factor of $10^{-7}$
suppression of the scattering rate for the electron shelving detection on the
neighboring ion for a beam focused to $2.6~{}\upmu\text{m}$ on the bath ion.

In such a setup, the dominant decoherence mechanism is arising from global
magnetic field fluctuations. This process can be expressed in terms of a jump
operator of the form $c^{\prime}=\sqrt{\kappa_{c}}\sum_{i}\sigma_{z}^{(i)}$,
assuming a decoherence rate of $\kappa_{c}=3.3$ Hz Ruster2016 . Fig. 6 shows
the cooling of such a system to an optimized ground state fidelity of $f=0.92$
in the decoherence-free case, while the presence of decoherence leads to a
fidelity of $f=0.89$.

An alternative to single ion addressing is to employ another isotope for the
bath ion, such as 44Ca+. The large isotope shifts of 850 MHz on the S1/2-P1/2
transition and 5.3 GHz on the qubit transition Solaro2018 ; Gebert2015  will
significantly relax the focussing requirements at the expense of having to
achieve an appropriately ordered ion crystal Splatt2009 .

##  III Discussion

Here, we demonstrated how adding a dissipatively driven auxiliary particle can
sympathetically cool a quantum simulator into low-energy states. Our approach
is efficient even when using only a single bath spin, and it exhibits strong
robustness against unwanted decoherence occurring in the quantum simulator.
Future directions include investigating the scaling behavior when optimally
varying the coupling constants of the bath in time and when adding multiple
bath spins. In the latter case, it will also be of interest to choose
different splittings of the bath spins, allowing engineering of tailored bath
spectral functions for the quantum simulator.

##  IV Materials and Methods

###  IV.1 Numerical simulations

All numerical simulations were performed using a wave-function Monte Carlo
approach provided by the QuTiP library Johansson2013 , extended to a massively
parallelized version Raghunandan2018 . Results were obtained by averaging over
1,000 Monte Carlo trajectories. We note that we are interested in the long
time limit of a weakly dissipative system, i.e., a regime where tensor network
algorithms are breaking down Kshetrimayum2017 . Numerical optimization of the
coupling constants $g_{sb}$ and $\gamma$ was carried out using a Nelder-Mead
algorithm. We typically obtain convergence within approximately 50 runs of the
simulation, which does not significantly depend on the size of the system.

##  V Acknowledgements

This work was funded by the Volkswagen Foundation and the DFG within SFB 1227
(DQ-mat, projects A01, A04, and B05).

##  VI Author contributions

M.R. and H.W. designed the cooling protocol. M.R. performed the numerical
simulations. F.W., C.O., and P.O.S. designed the experimental implementation
proposal. All authors contributed to the writing of the manuscript.

##  VII Competing interests

The authors declare that they have no competing interests.

##  VIII Data and materials availability

All data needed to evaluate the conclusions in the paper are present in the
paper and/or the Supplementary Materials. Additional data related to this
paper may be requested from the authors.

## References

  * (1) Georgescu, I. M., Ashhab, S. & Nori, F.  Quantum simulation.  _Rev. Mod. Phys._ 86, 153–185 (2014).
  * (2) Greiner, M., Mandel, O., Esslinger, T. W., Hänsch, T. & Bloch, I.  Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms.  _Nature_ 415, 39 (2002).
  * (3) Schneider, U. _et al._ Metallic and Insulating Phases of Repulsively Interacting Fermions in a 3D Optical Lattice.  _Science_ 322, 1520–1525 (2008).
  * (4) Jördens, R., Strohmaier, N., Günter, K., Moritz, H. & Esslinger, T.  A Mott insulator of fermionic atoms in an optical lattice.  _Nature_ 455, 204–207 (2008).
  * (5) Friedenauer, A., Schmitz, H., Glueckert, J. T., Porras, D. & Schaetz, T.  Simulating a quantum magnet with trapped ions.  _Nature Physics_ 4, 757–761 (2008).
  * (6) Lanyon, B. P. _et al._ Universal Digital Quantum Simulation with Trapped Ions.  _Science_ 334, 57– (2011).
  * (7) Aidelsburger, M. _et al._ Realization of the Hofstadter Hamiltonian with Ultracold Atoms in Optical Lattices.  _Phys. Rev. Lett._ 111, 185301 (2013).
  * (8) Miyake, H., Siviloglou, G. A., Kennedy, C. J., Burton, W. C. & Ketterle, W.  Realizing the Harper Hamiltonian with Laser-Assisted Tunneling in Optical Lattices.  _Phys. Rev. Lett._ 111, 185302 (2013).
  * (9) Álvarez, G. A., Suter, D. & Kaiser, R.  Localization-delocalization transition in the dynamics of dipolar-coupled nuclear spins.  _Science_ 349, 846–848 (2015).
  * (10) Mazurenko, A. _et al._ A cold-atom Fermi-Hubbard antiferromagnet.  _Nature (London)_ 545, 462–466 (2017).
  * (11) Bernien, H. _et al._ Probing many-body dynamics on a 51-atom quantum simulator.  _Nature_ 551, 579 (2017).
  * (12) Zhang, J. _et al._ Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator.  _Nature (London)_ 551, 601–604 (2017).
  * (13) Guardado-Sanchez, E. _et al._ Probing the Quench Dynamics of Antiferromagnetic Correlations in a 2D Quantum Ising Spin System.  _Phys. Rev. X_ 8, 021069 (2018).
  * (14) Lienhard, V. _et al._ Observing the Space- and Time-Dependent Growth of Correlations in Dynamically Tuned Synthetic Ising Models with Antiferromagnetic Interactions.  _Phys. Rev. X_ 8, 021070 (2018).
  * (15) Diehl, S. _et al._ Quantum states and phases in driven open quantum systems with cold atoms.  _Nature Phys._ 4, 878–883 (2008).
  * (16) Verstraete, F., Wolf, M. M. & Ignacio Cirac, J.  Quantum computation and quantum-state engineering driven by dissipation.  _Nature Phys._ 5, 633–636 (2009).
  * (17) Weimer, H., Müller, M., Lesanovsky, I., Zoller, P. & Büchler, H. P.  A Rydberg quantum simulator.  _Nature Phys._ 6, 382–388 (2010).
  * (18) Barreiro, J. T. _et al._ An open-system quantum simulator with trapped ions.  _Nature_ 470, 486 (2011).
  * (19) Carr, A. W. & Saffman, M.  Preparation of Entangled and Antiferromagnetic States by Dissipative Rydberg Pumping.  _Phys. Rev. Lett._ 111, 033607 (2013).
  * (20) Rao, D. D. B. & Mølmer, K.  Dark Entangled Steady States of Interacting Rydberg Atoms.  _Phys. Rev. Lett._ 111, 033606 (2013).
  * (21) Lin, Y. _et al._ Dissipative production of a maximally entangled steady state of two quantum bits.  _Nature_ 504, 415–418 (2013).
  * (22) Shankar, S. _et al._ Autonomously stabilized entanglement between two superconducting quantum bits.  _Nature_ 504, 419–422 (2013).
  * (23) Cormick, C., Bermudez, A., Huelga, S. F. & Plenio, M. B.  Dissipative ground-state preparation of a spin chain by a structured environment.  _New Journal of Physics_ 15, 073027 (2013).
  * (24) Morigi, G. _et al._ Dissipative Quantum Control of a Spin Chain.  _Phys. Rev. Lett._ 115, 200502 (2015).
  * (25) de Moraes Neto, G. D., Teizen, V. F., Montenegro, V. & Vernek, E.  Steady many-body entanglements in dissipative systems.  _Phys. Rev. A_ 96, 062313 (2017).
  * (26) Weimer, H., Müller, M., Büchler, H. P. & Lesanovsky, I.  Digital quantum simulation with Rydberg atoms.  _Quant. Inf. Proc._ 10, 885–906 (2011).
  * (27) Roghani, M. & Weimer, H.  Dissipative preparation of entangled many-body states with Rydberg atoms.  _Quantum Sci. Technol._ 3, 035002 (2018).
  * (28) Breuer, H.-P. & Petruccione, F.  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2002).
  * (29) Sachdev, S.  _Quantum Phase Transitions_ (Cambridge University Press, Cambridge, 1999).
  * (30) Nielsen, M. A. & Chuang, I. L.  _Quantum computation and quantum information_ (Cambridge University Press, Cambridge, 2000).
  * (31) Schollwöck, U., Richter, J., Farnell, D. J. & Bishop, R. F. (Editors).  _Quantum Magnetism_.  Lecture Notes in Physics, Vol. 645 (Springer, Berlin, 2004).
  * (32) Latorre, J. I., Rico, E. & Vidal, G.  Ground State Entanglement in Quantum Spin Chains.  _Quant. Inf. Comput._ 4, 48–92 (2004).
  * (33) Lemeshko, M. & Weimer, H.  Dissipative binding of atoms by non-conservative forces.  _Nature Commun._ 4, 2230 (2013).
  * (34) Kokail, C. _et al._ Self-verifying variational quantum simulation of lattice models.  _Nature_ 569, 355–360 (2019).
  * (35) Weimer, H., Yao, N. Y., Laumann, C. R. & Lukin, M. D.  Long-Range Quantum Gates using Dipolar Crystals.  _Phys. Rev. Lett._ 108, 100501 (2012).
  * (36) Jurcevic, P. _et al._ Quasiparticle engineering and entanglement propagation in a quantum many-body system.  _Nature_ 511, 202–205 (2014).
  * (37) Lloyd, S.  Universal Quantum Simulators.  _Science_ 273, 1073–1078 (1996).
  * (38) Kim, K. _et al._ Entanglement and Tunable Spin-Spin Couplings between Trapped Ions Using Multiple Transverse Modes.  _Physical Review Letters_ 103 (2009).
  * (39) Benhelm, J., Kirchmair, G., Roos, C. F. & Blatt, R.  Towards Fault-Tolerant Quantum Computing with Trapped Ions.  _Nat Phys_ 4, 463–466 (2008).
  * (40) Roos, C. F.  Ion Trap Quantum Gates with Amplitude-Modulated Laser Beams.  _New J. Phys._ 10, 013002 (2008).
  * (41) Tan, T. R. _et al._ Multi-element logic gates for trapped-ion qubits.  _Nature_ 528, 380–383 (2015).
  * (42) Ruster, T. _et al._ A long-lived Zeeman trapped-ion qubit.  _Applied Physics B_ 122, 254 (2016).
  * (43) Solaro, C., Meyer, S., Fisher, K., DePalatis, M. V. & Drewsen, M.  Direct Frequency-Comb-Driven Raman Transitions in the Terahertz Range.  _Phys. Rev. Lett._ 120, 253601 (2018).
  * (44) Gebert, F. _et al._ Precision Isotope Shift Measurements in Calcium Ions Using Quantum Logic Detection Schemes.  _Phys. Rev. Lett._ 115, 053003 (2015).
  * (45) Splatt, F. _et al._ Deterministic Reordering of 40Ca+ Ions in a Linear Segmented Paul Trap.  _New J. Phys._ 11, 103008 (2009).
  * (46) Johansson, J., Nation, P. & Nori, F.  QuTiP 2: A Python framework for the dynamics of open quantum systems.  _Comp. Phys. Comm._ 184, 1234–1240 (2013).
  * (47) Raghunandan, M., Wrachtrup, J. & Weimer, H.  High-Density Quantum Sensing with Dissipative First Order Transitions.  _Phys. Rev. Lett._ 120, 150501 (2018).
  * (48) Kshetrimayum, A., Weimer, H. & Orús, R.  A simple tensor network algorithm for two-dimensional steady states.  _Nature Commun._ 8, 1291 (2017).
  * (49) Kubo, R.  The Spin-Wave Theory of Antiferromagnetics.  _Phys. Rev._ 87, 568–580 (1952).
  * (50) Vidal, G. & Werner, R. F.  Computable measure of entanglement.  _Phys. Rev. A_ 65, 032314 (2002).
  * (51) Eisert, J., Cramer, M. & Plenio, M. B.  Colloquium : Area laws for the entanglement entropy.  _Rev. Mod. Phys._ 82, 277–306 (2010).
  * (52) Gross, D., Flammia, S. T. & Eisert, J.  Most Quantum States Are Too Entangled To Be Useful As Computational Resources.  _Phys. Rev. Lett._ 102, 190501 (2009).