# Gradient-based Neuromorphic Learning on Dynamical RRAM Arrays

Peng Zhou  Jason K. Eshraghian  Dong-Uk Choi  Wei D. Lu  and Sung-Mo Kang  P.
Zhou, D. Choi, and S. M. Kang are with the Department of Electrical and
Computer Engineering, UC, Santa Cruz, CA, USA.J. K. Eshraghian and W. D. Lu
are with the Department of Electrical Engineering and Computer Science,
University of Michigan, Ann Arbor, MI, USA.

###### Abstract

We present MEMprop, the adoption of gradient-based learning to train fully
memristive spiking neural networks (MSNNs). Our approach harnesses intrinsic
device dynamics to trigger naturally arising voltage spikes. These spikes
emitted by memristive dynamics are analog in nature, and thus fully
differentiable, which eliminates the need for surrogate gradient methods that
are prevalent in the spiking neural network (SNN) literature.

Memristive neural networks typically either integrate memristors as synapses
that map offline-trained networks, or otherwise rely on associative learning
mechanisms to train networks of memristive neurons. We instead apply the
backpropagation through time (BPTT) training algorithm directly on analog
SPICE models of memristive neurons and synapses. Our implementation is fully
memristive, in that synaptic weights and spiking neurons are both integrated
on resistive RAM (RRAM) arrays without the need for additional circuits to
implement spiking dynamics, e.g., analog-to-digital converters (ADCs) or
thresholded comparators. As a result, higher-order electrophysical effects are
fully exploited to use the state-driven dynamics of memristive neurons at run
time. By moving towards non-approximate gradient-based learning, we obtain
highly competitive accuracy amongst previously reported lightweight dense
fully MSNNs on several benchmarks.

###### Index Terms:

Neuromorphic computing, memristor, spiking neural network, supervised
learning, backpropagation.

##  I Introduction

Figure 1: An overview of our MEMprop approach. The MSNN architecture and the
resulting computational graph consisting of memristive dynamics.

Spiking Neural Networks (SNNs) impose neuroscience-inspired constraints to
modern deep learning algorithms, and have accordingly demonstrated significant
improvements in runtime efficiency. By moving from full precision and fixed
precision activations of artificial neuron models over to temporally-encoded
data representations captured by spiking neurons, neuromorphic hardware has
shown significant savings in energy consumption and latency [1, 2, 3, 4, 5,
6].

The broad success of error backpropagation to train deep learning models has
ushered in a plethora of related training algorithms adapted for SNNs, most of
which are guided by surrogate gradient descent to overcome the non-
differentiability of discrete spikes [7, 8]. This proliferation of SNN usage
is complemented by the development of modular deep learning programming
packages that have optimized autodifferentiation for CUDA acceleration [9, 10,
11, 12, 13].

In parallel to these advances in training SNNs, the past decade has seen huge
strides in brain-inspired devices, circuits, and architectures that integrate
neuronal dynamics to improve the hardware integration of SNNs and its
constituent parts. Memristors and resistive RAM (RRAM) make up an immense part
of such exploratory research in SNN implementation as they are a natural
bridge between SNN algorithms and accelerators [14, 15]. They have been widely
employed as both synapses and as spiking neurons.

At the ionic level, memristive synapses have been integrated into systems that
naturally implement the spike-timing-dependent-plasticity (STDP) update rule
using higher-order device dynamics [16, 17, 18]. An alternative use of ion-
driven dynamics is when implementing the memristor as a neuron, where
nonlinear conductance evolution gives rise to abrupt switching that can be
used to emit sudden voltage spikes. This approach is typically coupled with
capacitive integration, and has been referred to as a ‘neuristor’ [19, 20, 21,
22], and a ‘memristive integrate-and-fire’ (MIF) neuron [23, 24, 3].
Similarly, membrane leakage in biological neurons can be implemented using
resistive dissipation as in neuristors, or via volatile ionic drifting
dynamics observed in single devices [25] and also in nanowire networks [26,
27, 28].

Moving up to the architectural level, RRAM has been shown as a promising
candidate for compute-in-memory (CIM) architectures due to their ability to
parallelize matrix-vector multiplication independently of time complexity when
integrated as large-scale, modular arrays [29, 30, 31, 32, 33]. Rather than
neurons, memristive synapses map neural network weights to device
conductances. In general, RRAM CIM architectures are intended to be trained
offline with weights mapped on-chip for inference and deployment. As such,
RRAM synapses should be stationary and only used for weight read-out. Higher-
order dynamical behaviors of memristors are abstracted away, and treated as
non-idealities.

An additional challenge with RRAM-based CIM is the cost of communicating
analog current signals along lengthy bit-lines and conversion into the digital
domain. These issues have spurred the use of binary activations in the form of
spike-based CIM accelerators which have shown to alleviate the burdens of
mixed-signal computation, by removing the need for large analog-to-digital
(ADC) data converters [34].

The majority of deep learning acceleration using memristors can be classified
into one of the above cases: memristive neurons, memristive synapses that
learn via associative learning, and CIM accelerators. A small set of designs
have combined memristive neurons and memristive synapses together [35, 36].
Doing so is a commendable feat, as the natural switching dynamics of
memristive systems are relied on to achieve data-driven tasks. The cost of
allowing hardware behave naturally is that a designer can no longer rely on
synchronous, clock-driven processing, and is subject to fault injections that
are a result of nonlinear ionic dynamics. Letting the natural dynamics of
memristive hardware to ‘teach itself’ exacerbates the challenges of training
MSNNs. This challenge has limited the demonstration of fully memristive SNNs
to unsupervised learning tasks that have shown to solve simple, low-
dimensional pattern recognition via local learning rules (typically STDP) and
associative learning. These include the classification of several characters
and numbers (such as a subset of the MNIST dataset).

In this paper, we push beyond the classical limits of fully memristive SNNs,
and demonstrate competitive performance on real-world datasets that go beyond
static pattern recognition using a MSNN that includes both memristive synapses
and memristive neurons. We achieve this by combining higher-order memristive
dynamics with the gradient optimization process. While fault injections during
weight update have been accounted for in the past, this is the first time
nonlinear memristive dynamics are included within the computational graph used
in the gradient calculation process, enabling more precise fine-tuning of
network parameters.

The action potentials generated by memristive neurons are a result of
nonlinear ion-driven dynamics, which are functionally equivalent to a chain of
nonlinear operations in a computational graph. As each spike is a result of
naturally arising physics-driven phenomena rather than discrete switching in
digitally-composed SNNs, all functions are fully differentiable and the
gradient is well-defined. Therefore, there is no need to resort to surrogate
gradient approaches that have become ubiquitous when training SNNs. The
voltage action potential is scaled by memristive synapses, and the resultant
current drives downstream layers of memristive neurons.

An overview of our approach is shown in Fig. 1 and a summary of our
contributions are provided below:

  * •

We present a fully memristive SNN utilizing both memristive neurons and
synapses, and is shown to achieve state-of-the-art accuracy for lightweight
dense network architectures on real-world datasets. These datasets are more
complex than what has been demonstrated in the past on MSNNs while using
naturally arising physics-driven phenomena in RRAM;

  * •

We propose MEMprop, the application of error backpropagation directly to SPICE
circuit models of memristive neurons such that higher-order device dynamics
are fully utilized in the learning process. Memristive dynamics are broken
down into a series of composable, differentiable functions and used during
gradient descent;

  * •

SNNs are shown to be trainable without the need for gradient approximations
around hard thresholds, such as with surrogate gradient descent, and

  * •

The need for ADCs is amortized by relying entirely upon analog current
injection at the input, and analog action potentials at the output.

In doing so, the rich dynamics of nonlinear devices can be fully leveraged in
larger systems, in a way that moves far beyond applying fault injection and
variation into the training process. We achieve state-of-the-art accuracy for
lightweight dense network architectures on several static and neuromorphic
datasets, which pushes the fringe of what previous MSNNs have shown to be
capable of.

##  II Background

There is an incredibly broad span of work that integrates memristors with
brain-inspired architectures, from low-level analog action potential emulation
[23], to discrete spiking dynamics [34], up to non-spiking CIM processors
[29]. We focus our background on prior work that use nonlinear dynamics in
memristive neurons together with memristive synapses, with associated
demonstration of synaptic optimization to achieve a data-driven outcome.

###  II-A Fully Memristive SNNs

Full memristive SNNs refer to arrays that utilize nonlinear switching dynamics
of memristors to trigger action potentials, and are coupled with memristive
weights that are used as neural network parameters. The 8$\times$8 crossbar
array presented in [35] successfully integrates a fully memristive SNN,
including memristive synapses and neurons. The synaptic array is trained using
unsupervised STDP to classify four letters in a 24-pixel grid. While the task
achieved is considerably simple, the fully memristive experimental
demonstration sets the stage to build upon new training methods.

The work in [37] uses half-wave rectification interposed between crossbar
arrays to process the ReLU activation in the analog domain. While not fully
memristive, nor a ‘spiking’ network, it demonstrates a pivotal example of
successively passing analog activations between RRAM crossbars without
intermediate data conversion, in a manner akin to how our analog action
potentials are transmitted between layers. The training process relies on
gradient-based optimization where device non-idealities are injected during
the forward pass. In doing so, a test set accuracy of 93.63$\%$ is obtained on
the MNIST dataset.

In Refs. [38] and [39], convolutional SNNs with memristors are used, where
both used pretrained on non-spiking networks that are mapped or converted into
the spiking domain. Both networks obtained competitive accuracy on the MNIST
dataset, though did not demonstrate performance on more complex, real-world
data. This may potentially be due to the large differences between the
networks that were trained and the MSNN that was implemented. In Ref. [40], a
dense MSNN is adopted using a similar approach to what is used here, and as
such, has minimal hardware requirements at run-time. The training process
abstracts away the switching dynamics of the memristive neuron into a firing
rate, and may be the reason why a relatively low accuracy of 83.2% was
achieved on the MNIST dataset. A more detailed comparison is illustrated in
Section V, subsection B.

Almost all of these works offer compelling demonstrations using in-house
fabricated arrays, either with standalone crossbars or back-end-of-the-line
(BEOL) integrated arrays with foundry-made chips. In contrast, we have aimed
to make our work as device-agnostic as presently possible by using
commercially available Knowm memristors and its corresponding model [41].

Although Knowm memristors are not known for their reliability, their
metastable switching dynamics are accounted for within the gradient
calculation step, as it forms part of the computational graph. In contrast,
Kiani et al. use the memristors in the forward-pass only, as their devices are
not intended to be reprogrammed during inference, and thus do not require
switching to generate spiking dynamics. Gradients can therefore be
deterministically calculated partially off-chip [37], whereas our approach
harnesses memristive dynamics in the forward-pass computation in our network.
As such, our MSNN approach can leverage the benefits of spike-based
processing, such as sparse processing and lower data collision rates.

###  II-B Memristive Learning Frameworks

To ease and emulate the training process of memristive networks, a variety of
valuable frameworks have been developed that each address various niches.
These include MemTorch [42], NeuroSim [43], and the IBM Analog Hardware
Acceleration Kit [44], which implement non-spiking networks that adopt mixed-
signal bit-line charge/current accumulation/summation processing. In these
simulators, memristive dynamics are accounted for during weight updates, and
are otherwise fixed during inference. To complement these tools, NeuroPack
[45] specifically targets the simulation of spiking networks, where memristive
dynamics are also factored in during the weight update process, and fixed
during inference. Spiking dynamics are triggered by pulse-based input
voltages.

The closest relation to our proposed work comes from Demirag et al. who offer
a software implementation of a real-time recurrent learning variant [46, 47]
using phase change device models to train a MSNN [48].

Much like these simulators, our approach MEMprop integrates SPICE-level
memristor models into the training process. But distinct from these simulators
and training methods, we utilize arrays of memristive integrate-and-fire (MIF)
neuristors that trigger self-induced action potentials during the forward-pass
computation. To emit spikes, the state of a memristor must evolve and
ultimately switch to provide sudden discharge pathways which triggers spikes.
Therefore, the dynamical characteristics of a memristor will alter the
gradient itself, whereas all prior approaches account for dynamics in the
weight update step when reprogramming memristive synapses.

Our approach applies the backpropagation algorithm directly to SPICE-level
neuristor models that emit analog action potentials as a result of charge-
accumulation and metastable threshold switching. That is to say, the
memristive dynamics in SPICE models are not fixed during the forward-pass, but
rather, they are dynamically extracted during the forward-pass to fire
neuristor-induced spikes. The BPTT algorithm is applied directly to SPICE
models, thus integrating memristive neurons as part of the gradient
calculation step, rather than isolating their dynamics to synaptic weight
updates as in prior implementations. This offers a dynamical and totally new
way to train fully memristive networks that uses state-based dynamics as part
of the gradient calculation step.

In terms of hardware implementation, the conventional use of RRAM in circuits
often requires significant overhead to convert analog currents into digital
voltages and consumes large amount of power [29]. In many instances, the
power/area demands of the ADCs and and Digital-to-Analog Converters (DAC) far
exceeds the overhead brought on by RRAM, offsetting the advantages of
memristors. In contrast, our spike-based approach eliminates the need for ADCs
and DACs such that the cost of peripheral circuits are substantially reduced.

###  II-C Error Backpropagation Through SPICE models

Our proposed approach enables us to scale up the complexity of learnable tasks
in fully MSNNs by directly applying gradient descent to the nonlinear state
evolution of memristive neurons and synapses. Both neurons and synapses in
biological neural networks are modeled using memristors. The MIF neuron model
is designed to achieve distinct depolarization, hyperpolarization, and
repolarization voltage phases with a minimal set of circuit elements.
Memristive synapses act as interconnects between layers of neurons.

To train the fully memristive network, we propose MEMprop, an application of
error backpropagation to large-scale networks derived from SPICE circuit
models. This enables dynamical, time-varying memristive neurons to learn and
thus achieve much higher accuracy on data-driven tasks than has been
previously reported with MSNNs. By relying on the analog spiking
characteristics that naturally occur in the MIF neuron model, the non-
differentiability of spike-based activations are completely avoided. This
means MEMprop does not rely on surrogate gradient techniques that are commonly
used to train SNNs, which calculates biased gradient estimators to circumvent
the dead neuron problem [8]. To promote broad accessibility of our methods, we
use SPICE models of commercially available, low-cost memristors and
demonstrate the efficacy of MEMprop in a supervised deep learning framework.

##  III Methods

###  III-A Memristive Integrate-and-Fire Model

Table I: MIF model differential equations vs numerical integration. Variable | Continuous Time Derivative | Discrete Time Solution
---|---|---
$v$ | $\frac{{\mathrm{d}}v}{\mathrm{d}t}=\frac{I-G_{1}(v-E_{\rm rest})-G_{2}(v-E_{\rm reset})}{C}$ | $v^{t+1}=\frac{I_{t}-G_{1}^{t}(v^{t}-E_{\rm rest})-G_{2}^{t}(v^{t}-E_{\rm reset})}{C}+v^{t}$
$G_{1}$ | $G_{1}=\frac{x_{1}}{R_{\rm on1}}+\frac{1-x_{1}}{R_{\rm off1}}$ | $G_{1}^{t+1}=\frac{x_{1}^{t+1}}{R_{\rm on1}}+\frac{1-x_{1}^{t+1}}{R_{\rm off1}}$
$G_{2}$ | $G_{2}=\frac{x_{1}}{R_{\rm on2}}+\frac{1-x_{2}}{R_{\rm off2}}$ | $G_{2}^{t+1}=\frac{x_{2}^{t+1}}{R_{\rm on2}}+\frac{1-x_{2}^{t+1}}{R_{\rm off2}}$
$x_{1}$ | $\frac{{\mathrm{d}x_{1}}}{\mathrm{d}t}=\frac{1}{\tau_{1}}(\frac{1-x_{1}}{1+\exp(\frac{v_{\rm on1}-(v-E_{\rm rest})}{V_{\rm th}\cdot k_{\rm v}})}-\frac{x_{1}}{1+\exp(\frac{(v-E_{\rm rest})-v_{\rm off1}}{V_{\rm th}\cdot k_{\rm v}})})$ | $x_{1}^{t+1}=\frac{1}{\tau_{1}}\Big{(}\frac{1-x_{1}^{t}}{1+\exp(\frac{v_{\rm on1}-(v^{t}-E_{\rm rest})}{V_{\rm th}\cdot k_{\rm v}})}-\frac{x_{1}^{t}}{1+\exp(\frac{(v^{t}-E_{\rm rest})-v_{\rm off1}}{V_{\rm th}\cdot k_{\rm v}})}\Big{)}+x_{1}^{t}$
$x_{2}$ | $\frac{{\mathrm{d}x_{2}}}{\mathrm{d}t}=\frac{1}{\tau_{2}}(\frac{1-x_{2}}{1+\exp(\frac{v_{\rm on2}-(v-E_{\rm reset})}{V_{\rm th}\cdot k_{\rm v}})}-\frac{x_{2}}{1+\exp(\frac{(v-E_{\rm reset})-v_{\rm off2}}{V_{\rm th}\cdot k_{\rm v}})})$ | $x_{2}^{t+1}=\frac{1}{\tau_{2}}\Big{(}\frac{1-x_{2}^{t}}{1+\exp(\frac{v_{\rm on2}-(v^{t}-E_{\rm reset})}{V_{\rm th}\cdot k_{\rm v}})}-\frac{x_{2}^{t}}{1+\exp(\frac{(v^{t}-E_{\rm reset})-v_{\rm off2}}{V_{\rm th}\cdot k_{\rm v}})}\Big{)}+x_{2}^{t}$

The neuron model adopted in our MSNN is the MIF neuron model depicted in Fig.
2(a) [23]. Qualitatively, given a positive current injection $I(t)$, the
membrane potential $v(t)$ will charge up from the neuron’s resting potential
$E_{\rm rest}$. Once a sufficiently large electric field builds up across the
memristor $M_{2}$, it switches on, effectively shorting the output to $E_{\rm
reset}$ which charges back up to equilibrium, $E_{\rm rest}$. The evolution of
the membrane potential is shown in Fig. 2(b), and our prior work in Ref. [23]
provides an experimental demonstration using a pair of Knowm memristors.

The memristor model used is based on the generalized metastable switch (MSS)
model [41], which has been used to accurately describe a large range of
possible devices. In this model, a MSS is an idealized element that switches
with a given probability between two states as a function of its voltage and
temperature. A memristor is modeled by a collection of MSSs, which determines
the state-time dynamics that lead to non-volatile characteristics. $x$
characterizes the internal state as a variable normalized between 0 and 1, as
determined by the switching states of all MSSs.

Formally, the governing dynamics of the MIF neuron are characterized by the
system of differential equations in Table I. The membrane potential $v$ is
dependent on $G_{1}$ and $G_{2}$, the device conductances of $M1$ and $M2$,
the MIF capacitance $C$, and $E_{\rm rest}$ and $E_{\rm reset}$, which are
voltage biases in the MIF circuit.

(a)

(b)

Figure 2: (a) A memristive integrate-and-fire (MIF) neuron consists of two
memristors $M_{1}$ and $M_{2}$, connected to DC voltage sources $E_{\rm rest}$
and $E_{\rm reset}$, in parallel with a capacitor $C$. The MIF neuron is
provably minimal in generating a membrane potential traversing from a rest
voltage level to a threshold voltage level, then to a reset voltage level, and
then back to the rest potential when a current pulse is applied. (b) An action
potential waveform that MIF model resembles showing rest, threshold and reset
potentials.

$G_{1}$ and $G_{2}$ are dependent on $x_{1}$ and $x_{2}$, a pair of variables
governing the internal state of each device, where $R_{\rm on}$ and $R_{\rm
off}$ are the on/off resistances.

The rates of change of $x_{1}$ and $x_{2}$ are dependent on the state
evolution time constants $\tau_{1}$ and $\tau_{2}$, the thermal voltage
$V_{\rm th}$ at room temperature, and the volatility constant $k_{\rm
v}\in[0,1]$ in the range of 0 to 1. As $k\rightarrow 0$, the more retentive
the device is.

###  III-B Neural Network Layout

When accounting for the hardware implementation of a fully memristive neural
network, the voltage response of the MIF neuron must drive the memristive
synapses, and is correspondingly weighted by the synaptic conductances. This
is can be fully integrated in a crossbar according to Eq. (1) below:

$\mathbf{I}=\mathbf{G}\times\mathbf{V}$ (1)

$\displaystyle\begin{bmatrix}I_{1}\\\ I_{2}\\\ \vdots\\\ I_{n}\end{bmatrix}$
$\displaystyle=\left[\begin{array}[]{cccc}w_{00}&w_{01}&\cdots&w_{0m}\\\
w_{10}&w_{11}&\cdots&w_{1m}\\\ \vdots&\vdots&\ddots&\vdots\\\
w_{n0}&w_{n1}&\cdots&w_{nm}\\\ \end{array}\right]\begin{bmatrix}v_{1}\\\
v_{2}\\\ \vdots\\\ v_{m}\end{bmatrix}$

$I_{n}=\sum_{i=1}^{m}v_{i}\times w_{ni}$ (2)

where $\mathbf{I}$ is the output current vector, $\mathbf{V}$ is the input
voltage vector, and $\mathbf{G}$ is the conductance matrix of the crossbar.
The output current vector is generated via bit-line current summation as shown
in Eq. (2), and directly drives the input of the next MIF neuron layer.
Resistive loading may attenuate the output current in subsequent stages, but
this can be accounted for using a scaling factor, or otherwise buffered [49,
50, 23].

Fig. 3 shows a small-scale schematic of a fully memristive neural network with
five MIF neurons and a $5\times 3$ memristive crossbar, which receives five
fan-in input currents $I_{\rm in0}$ \- $I_{\rm in4}$ and generates three fan-
out output currents. $E_{\rm rest}$ and $E_{\rm reset}$ are the voltage
sources in the MIF circuit. $C_{\rm BL}$ is the bit-line parasitic capacitance
generated by the metal line, or otherwise by the drain-bulk capacitance of
select transistors (not shown), and is used as the membrane capacitance in the
MIF circuit.

Figure 3: Schematic of a fully memristive neural network. $C_{\rm BL}$
represents the parasitic bit-line capacitance.

###  III-C Forward-mode Solution of MIF Model

In order to train a network of MIF neurons and synapses using gradient
descent, the differential equations representing the MIF circuit dynamics
(middle column, Table I) are recast into discrete-time form (left column,
Table I). In doing so, the memristive dynamics can be captured in a
computational graph that evolves over time, much like a recurrent neural
network (RNN).

In practice, SPICE simulators use a variety of differential equation solvers,
such as the backward Euler method, and the 4th-order Runge-Kutta method (RK4).
For compatibility with the BPTT algorithm, we solve the differential equations
using the forwrd Euler method, which provides an explicit representation of
the next time step using present-time dynamics. The rich dynamics of the MIF
neuronal network are now accounted for in the MSNN, unrolled in time such that
gradient descent can be used to optimize the memristive synapses as a function
of the MIF evolution.

Many neural coding studies represent spike trains as a summation of time-
shifted Dirac delta pulses $\sum_{n}{\delta(t-t_{j}^{n})}$. As such a model of
spikes are an idealization, we use the spike train to modulate a time-
continuous, alpha-shaped input current $I$ modeled by the equation in Table II
to be compatible with real, physical systems. In this series of equations, $a$
is an internal state variable, $\tau_{\rm syn}$ is a time constant that
determines the shape of the alpha current, $W_{i}$ is the synaptic weight
between a pre-synaptic neuron $i$ and its associated post-synaptic neuron. It
is commonly regarded that such alpha-shaped waveforms correspond to the
response from biological neurons in the sensory periphery that respond with
graded potentials [51].

Table II: Alpha current differential equations vs numerical integration. Variable | Continuous Time Derivative | Discrete Time Solution
---|---|---
$I$ | $\tau_{\rm syn}\frac{{\mathrm{d}I}}{\mathrm{d}t}=a-I$ | $I^{t+1}=\frac{a^{t}-I^{t}}{\tau_{\rm syn}}+I^{t}$
$a$ | $\tau_{\rm syn}\frac{{\mathrm{d}a}}{\mathrm{d}t}=-a+W_{i}\cdot\sum_{n}{\delta(t-t_{i}^{n})}$ | $a^{t+1}=-\frac{a^{t}+W_{i}\cdot\sum_{n}{\delta(t-t_{i}^{n})}}{\tau_{\rm syn}}+a^{t}$

###  III-D MIF Single Neuron Simulation

To verify the accuracy of the forward Euler method, we conduct a single neuron
simulation using Python shown in Fig. 4 across 1,000 time steps, which
provides ample time for the membrane potential $v$ to traverse from spiking
$V_{\rm th}$, to the reset potential $E_{\rm reset}$, back to the resting
potential $E_{\rm rest}$. Each time step is of duration 0.01 ms over a total
duration of 10 ms.

The parameters used in the simulation are listed in Table III. The resting
potential, charge integration, thresholding, and reset dynamics are closely
reproduced in a SPICE simulator that uses the RK4 solver, which verifies the
solution generated by adopting the forward Euler method. Deep learning is
known to be tolerant to fault injections [52, 53], and so the tradeoff between
the numerical integration accuracy and the training complexity should be
considered. We adopted Forward Euler method which is able to balance the two.
Technically, our training method is expected to generalize to other numerical
integration methods as well, but with added computational complexity.

(a)

(b)

(c)

Figure 4: Simulation results of MIF neuron solved using forward Euler numerical integration. The results matches the quantitative dynamics of solved by a SPICE simulator. Difference between two methods is marked by grey area. (a) Alpha synaptic current dynamics. (b) Internal states $x_{1}$, $x_{2}$. (c) Voltage response $v$. Table III: MIF circuit parameters Parameter | Value | Parameter | Value
---|---|---|---
$v_{\rm off_{1}}$, $v_{\rm off_{2}}$ | 5 mV |  $v_{\rm on_{1}}$, $v_{\rm on_{2}}$ | 110 mV
$R_{\rm off_{1}}$, $R_{\rm off_{2}}$ | 0.1 M$\Omega$ |  $R_{\rm on_{1}}$, $R_{\rm on_{2}}$ | 1 k$\Omega$
$\tau_{1}$, $\tau_{2}$ | 1 ms | $\tau_{\rm syn}$ | 0.64 ms
$E_{\rm rest}$ | 0 mV | $E_{\rm reset}$ | 50 mV
$V_{\rm th}$ | 25 mV | $k_{\rm v}$ | 0.6
$C$ | 100 pF |  |

###  III-E BPTT in MSNNs

The discrete-time solution in Table I is illustrated as a directed, acyclic
graph in Fig. 1, where time flows from left to right. The MIF circuit
parameters are color-coded to show how each electrical characteristic impacts
the others at the next time step.

To train a network, a loss function is calculated using the membrane potential
$v$ of the output layer at each step. Note that the adjoint method in [54] is
not adopted here, as an intermediate state is required to calculate loss and
guide training for each time step. We are concerned with the spiking output at
all time steps, and not just the final state of the system which makes BPTT a
more optimal choice. In our example network provided in Fig. 1 with 10 output
neurons, the predicted MIF neuron is expected to spike most frequently by
aiming to increase the membrane potential across time steps, while the
incorrect target should be suppressed. As the membrane dynamics are
continuous, a fully analog MSNN can be trained without surrogate gradients, as
has become ubiquitous in deep SNNs trained via error backpropagation [8, 55].

The BPTT algorithm iteratively applies the chain rule from the output back to
the leaf nodes, $w$, to determine the update direction to optimize the
network. Prior demonstrations treat $w$ as the device conductance. In our
case, we also do this, but additionally generate spiking behaviors using
naturally occuring MIF dynamics rather than applying a hard threshold to the
membrane potential. Not only is this a more biorealistic representation, but
it can be fully integrated using RRAM crossbars.

##  IV Experimental Results

###  IV-A Network Architecture

To validate the use of MIF neurons in a deep learning framework, we used a
lightweight 3-layer dense SNN with 100 hidden MIF neurons (bottom-left of Fig.
1). Each MSNN is simulated for a duration of 10 ms over a span of 1,000
discrete time steps.

The input layer consists of a number of input features (784 for the static
datasets; 2,048 for a downsampled neuromorphic dataset). An alpha current
generator with an amplitude weighted by the input pixel intensity (Table II),
weighted by memristive synapses, inducing 100 MIF neurons to fire into another
set of memristive synapses, terminated by the output layer of MIF neurons. To
account for circuit loading effects, the current injected to each MIF layer is
attenuated by a scaling factor defined in the hyperparameters.

###  IV-B Datasets

Three datasets of increasing difficulty are used to assess the MSNN: MNIST,
FashionMNIST, and the DVS128 Gesture datasets. Despite being considered a
‘solved’ problem for quite some time now, it was often the case that the MNIST
dataset was the most challenging problem that could be solved by previously
reported MSNNs. In most cases, a subset or simplified alternative would be
used to demonstrate pattern recognition. Here, we demonstrate for the first
time that RRAM arrays that rely on internal dynamics are still capable of
processing more challenging problems and real-world datasets, namely,
FashionMNIST and neuromorphic data.

The MNIST [56] and the slightly more challenging FashionMNIST datasets [57]
are used to test the performance of the MSNN on temporally static data. Both
datasets consist of 60,000 28$\times$28 greyscale images in the training set,
and 10,000 images in the test set, with 10 output classes of handwritten
digits (MNIST) and clothing items/accessories (FashionMNIST/FMNIST). During
the training process, the alpha current inject is applied at every 100 steps
to promote sparse network activity.

For a more challenging test case, the DVS128 Gesture dataset [58] is used a
neuromorphic baseline to test the MSNN’s performance on event-driven data
filmed with an event-based camera [59]. Each sample only processes sufficient
changes in luminance, and consists of 11 different output classes of hand
gestures, such as clapping, arm rotation, and air guitar. Each sample is
downsampled to a resolution of $32\times 32\times 2$, where the channel depth
of 2 accounts for on and off spikes (positive and negative luminance changes).
The training data is integrated to fit within 100 discrete time steps, and the
testing data is integrated within 360 time steps, to account for GPU memory
constraints.

###  IV-C Training Process

The forward Euler solution of each MIF neuron was defined in PyTorch v1.10.1
in Python 3.8, where the autodifferentiation framework was used to keep track
of gradients of all forward-mode computations on-the-fly.

The loss used is the negative log-likelihood of the membrane potential $v(t)$
(Fig. 1) of the output layer. The network objective is to increase voltage
across the MIF neuron associated to the correct class. The total loss across
time-steps is summed prior to backpropagating the error through the SPICE
memristor model using the Adam optimizer [60]. The MSNNs are each trained for
50 epochs.

###  IV-D Results

The final test set accuracy is measured across 50 epochs over 5 trials for
each dataset with early stopping applied, and the average result is shown in
Fig. 5 with error bars.

The average accuracy on the MNIST dataset reaches 93.08%, which is very high
among fully MSNN considering almost no required peripheral circuitry for
hardware implementation and far exceeds other MSNN baselines that have been
previously reported. Although it remains somewhat below non-memristive
baselines which can obtain greater than 97% on the MNIST dataset using a
similar architecture, such models are typically computed using fixed- and
full-precision arithmetic rather than depending on naturally arising
memristive dynamics.

Where our approach tends to shine is on real-world data that goes beyond the
simplicity of handwritten digit recognition. Our fully MSNN performance holds
for the FashionMNIST dataset, where we obtained an average of 84.77%, and also
for the DVS128 Gesture dataset, with 82.63%. This provides the first
successful demonstration of training a MSNN on a neuromorphic dataset by
relying on naturally occurring device dynamics for spike emission.

Figure 5: Test set accuracies with error bars over 5 trials.

##  V Discussion and Conclusion

###  V-A Area, Power, and Latency

With our fully MSNN approach, one of the main advantages is that we avoid the
need for ADCs, and subsequently, we do not need to multiplex or serialize
data. To quantify the potential improvement in terms of power, area, and
latency, we assume a dense 3-layer architecture of 784-100-10 neurons. Such a
network requires 110 MIF neurons and 79.4k memristive synapses. For each
synaptic weight, a pair of devices are required to implement current
subtraction to represent both positive and negative weights. Given the array
structure in Fig. 3, a total of 0.16M RRAM cells are required, which needs 10
RRAM tiles of 128$\times$128 in size.

To provide an area estimation at a 65-nm technology node, the size of each
RRAM tile with 1T1R cells is 2.77×10-3mm2, and a complete array will occupy
2.77×10-2mm2 [61]. If a bit-line current summation approach was adopted for a
dense ANN with 8-bit networks, a current-mode SAR ADC each occupies 3×10-3mm2,
where 4 ADCs are shared across the column wires for each crossbar. In such a
case, the improvement of our approach without converting and serializing
activations by adopting a spike-based approach is a factor of 5.33$\times$.

Estimating the power consumption of asynchronous, spike-based workloads
requires profiling the spiking activity in the network for a given dataset.
Using the MNIST dataset, we found that an average of 2% of neurons are spiking
at any given time in the network. This may be modulated to promote varying
degrees of excitation in the network, for example, by applying a spike-
dependent regularization term to the objective function. The average
resistance in our network is close to the midpoint between the on/off
resistances: $R_{\rm off}$ = 100 k$\Omega$, $R_{\rm on}$ = 1000 $\Omega$,
average $R_{\rm ave}$ = 50.5 k$\Omega$. The voltage across each device when
emitting a spike follows a sharp peak of $v_{\rm on}$ = 110 mV, a refractory
period where the potential drops down to $v_{\rm off}$ = 5 mV, and the average
value assuming the spike is linearized for the window of interest is $v_{\rm
ave}$ = 57.5 mV can give an estimate of the average power

$P_{\rm ave}=\frac{v_{\rm ave}^{2}}{R_{\rm ave}}\times N_{\rm cells}\times
2\%,$ (3)

which is then numerically evaluated for a single 128$\times$128 tile to be
21.45$\mu$W, and 0.21mW when accounting for all 10 tiles. When accounting for
data conversion overhead where each 8-bit ADC consumes 2×10-4W, the total
power consumption increases to 8.21 mW. Our approach can offer an improvement
of 38.30$\times$.

We estimate latency by driving ADCs at an operating frequency of 40MHz. With 4
ADCs converting 64 column currents (assuming 2 bit-lines per activation), 32
data-lines are grouped to share a driving DAC which generates an 8-bit
serialized input operating across 8 distinct cycles, evaluating to input data
latency of 6.4$\mu$s. At the output, 32 bit-lines each share an ADC where 16
column currents must be converted, increasing the latency by 0.4$\mu$s.
Further assuming that data conversion dominates latency, this means that it
takes 6.8$\mu$s per vector matrix multiplication (VMM). A 1,000 time step
simulation, as used in our network, will extend this to a latency of 6.8 ms.
Additionally, each input is fed as an alpha current signal with a measurable
latency around 0.64 ms, giving the total latency for the conventional approach
of 7.44 ms, as opposed to MEMprop which only requires the additional switching
time of the memristors (on the order of 100s of nanoseconds added to the alpha
current latency). In our approach, currents and voltages are generated without
the need for conversion or serialization, and so there is no added ADC
latency. The alpha current injection and response already account for RC delay
through the memory array.

A summary of the above mentioned comparison is shown in Table IV.

Table IV: Power, area, and latency improvement in our MSNN Aspect | Our method | Mixed-Signal | Improvement
---|---|---|---
Area | 2.77 mm2 | 14.77 mm2 |  5.33×
Power | 0.21 mW |  8.21 mW | 38.30×
Latency | 0.64 ms |  7.44 ms | 11.63x

###  V-B Comparison

A comparative analysis against other memristive networks is provided in Table
V. Where metrics are not provided in the original sources (e.g., chip area,
power consumption, and inference latency), we have made an estimate where
sufficient data has been provided to extrapolate these values. The closest
MSNNs to our work are briefly described below.

In Ref. [38], a convolutional neural network with analog neurons and digital
spiking at the input layer is adopted. Each analog neuron consists of four
transmission gates, seven operational amplifiers, a comparator, and a
memristor. While this is substantially more complex than the MIF neuron, the
large overhead per neuron is offset by using time-division multiplexing access
(TDMA) to treat a single physical neuron as multiple algorithmic neurons.

Table V: Comparison among fully memristive spiking neural networks

Method | Architecture | Neuron No. | Refractory | Neuron feature | Training Method | MNIST | Complexity of Additional
---|---|---|---|---|---|---|---
of elements | Control Logic
[38] | 1Mem-C-M-C- | 4(5)+33T | ✗ | Mixed-signal | LIF rate fit |  | high
M-1024F-10F | TDMA | 97.1%
[39] | 784-C-M-C-M-10F | 5+2T | ✗ | Digital | Stochastic switch fit |  | high
ANN converts SNN | $\sim$96%
[40] | 784-100F-10F | 3 | ✗ | Analog | LIF shape fit |  | low
surrogate gradient | 83.2%
Our work | 784-100F-10F | 3(5) | ✓ | Analog | MIF direct train |  | low
MEMprop | 93.2%

  * 1

T: transistor. 2 Mem: memristor. 3 C: convolutional layer. 4 M: max pooling
layer. 5 F: fully connected layer.

Table VI: Test set accuracies for the MNIST, FashionMNIST, and DVS128 Gesture
Datasets.

Dataset | Batch | LR | Best | Avg. ($n=5$) | $\sigma$
---|---|---|---|---|---
MNIST | 128 | 1e-4 | 93.18 | 93.08 | 0.07
FMNIST | 128 | 1e-4 | 84.79 | 84.77 | 0.13
DVS-128 | 16 | 1e-4 | 83.21 | 82.63 | 0.95

  * 1

LR: learning rate $\eta$. 2 $\sigma$: standard deviation.

Wijesinghe et al. also adopt a convolutional SNN using stochastic switching to
implement probabilistic firing, which is a natural fit for devices that switch
based on random processes [39]. A non-spiking ANN is first trained and
subsequently converted to an SNN. This approach is known to perform optimally
for static datasets on reasonably deep SNNs, but it sets an upper-limit of
performance such that the SNN accuracy will not surpass the accuracy of the
ANN. Furthermore, ANN to SNN conversion is yet to show success on neuromorphic
datasets [55]. The total area of the design is reported to be 3 mm2, and the
latency for a single spike is 0.21 $\mu$s.

In Ref. [40], a 3 element neuron is constructed, consisting of a resistor, a
memristor, and a capacitor, which is closest in spirit to our own approach
here. The neuron is parameterized to fit to a leaky integrate-and-fire neuron
model, which relies on hard-thresholded spike generation which is a non-
differentiable function. The problem is circumvented by utilizing surrogate
gradient descent, and results in a test set accuracy of 83.2%. This
comparatively lower accuracy highlights the challenges of adopting neuristor-
like dynamics, and the difficulty of mapping conventional training methods,
both supervised (surrogate gradient descent) and unsupervised (STDP), to
dynamical RRAM arrays. We expect the lack of a direct correspondence between
thresholded leaky integrate-and-fire neurons and memristive neuron spiking is
what resulted in the accuracy degradation, though it sets the stage for
developing new training methods that account for dynamical, continuous-time
switching into the computational graph of the network. In absence of power
consumption metrics, we apply the same assumptions as in our MSNN
architecture, namely, that 2% of the network is active at any given time, and
estimate the total power consumption is 0.63 W.

No surrogate gradients or associated approximations are required due to the
continous-time nature of the spiking behavior in our experiments. The absence
of ADCs/DACs in conventional full-precision/fixed-precision bit-line current
summation approaches reduces the computational overhead. This is due to the
use of MEMprop. More details of our experiments are listed in Table VI.

###  V-C Concluding Remarks

To the best of our knowledge, this is the first work that has mapped low-
level, analog SPICE dynamics directly into an acylic computational graph that
is compatible with the backpropagation algorithm. Directly mapping the
behavior of dynamical RRAM arrays into an autodifferentiation framework has
enabled us to achieve promising performance on real-world datasets on fully
memristive SNNs, beyond simple pattern recognition and handwritten digit
classification. By harnessing a blend of the intrinsic behaviors of memristors
and the sparse network activity of SNNs, the power consumption and latency of
our approach surpasses that of all other similar methods without compromising
on accuracy.

We propose MEMprop with a lightweight dense neural network to validate this
idea of applying backpropagation directly to SPICE circuit model. In the
future, a deeper convolutional MSNN can be explored to further improve
accuracy and to solve more complex tasks. In such a case, the benefits of
MEMprop in terms of area, power, latency should still hold. Our approach can
account for fault injections by incorporating device-to-device variability at
the SPICE model level, which are then ported into gradient-based updates. This
is important in the commercial memristors available, as learning heterogeneous
device dynamics can improve the representational capacity of MSNNs. MEMprop
should be able to offer a method by which these variations are accounted for
during the training process.

While variation is tolerable, neuristor-based firing relies on device
switching which is a limiting factor in RRAM cells that exhibit limited cycle
endurance. This motivates the development of devices that have a wide range of
behaviors, but are resilient to consistent switching. Nonetheless, we
successfully trained our network to be sparsely activated, which reduces the
frequency of switching and relaxes this stringent demand on large device
endurance.

## References

  * [1] M. R. Azghadi, C. Lammie, J. K. Eshraghian, M. Payvand, E. Donati, B. Linares-Barranco, and G. Indiveri, “Hardware implementation of deep network accelerators towards healthcare and biomedical applications,” _IEEE Transactions on Biomedical Circuits and Systems_ , vol. 14, no. 6, pp. 1138–1159, 2020.
  * [2] X. Wang, X. Lin, and X. Dang, “Supervised learning in spiking neural networks: A review of algorithms and evaluations,” _Neural Networks_ , vol. 125, pp. 258–280, 2020.
  * [3] P. Zhou, D.-U. Choi, J. K. Eshraghian, and S.-M. Kang, “A fully memristive spiking neural network with unsupervised learning,” _arXiv preprint arXiv:2203.01416_ , 2022.
  * [4] P. Zhou and S. Abbaszadeh, “Towards real-time machine learning for anomaly detection,” in _2020 IEEE Nuclear Science Symposium and Medical Imaging Conference (NSS/MIC)_. IEEE, Octorber 2020, pp. 1–3.
  * [5] A. Tavanaei, M. Ghodrati, S. R. Kheradpisheh, T. Masquelier, and A. Maida, “Deep learning in spiking neural networks,” _Neural Networks_ , vol. 111, pp. 47–63, 2019.
  * [6] T. Wunderlich, A. F. Kungl, E. Müller, A. Hartel, Y. Stradmann, S. A. Aamir, A. Grübl, A. Heimbrecht, K. Schreiber, D. Stöckel _et al._ , “Demonstrating advantages of neuromorphic computation: a pilot study,” _Frontiers in Neuroscience_ , vol. 13, p. 260, 2019.
  * [7] P. J. Werbos, “Backpropagation through time: what it does and how to do it,” _Proceedings of the IEEE_ , vol. 78, no. 10, pp. 1550–1560, 1990.
  * [8] E. O. Neftci, H. Mostafa, and F. Zenke, “Surrogate gradient learning in spiking neural networks: Bringing the power of gradient-based optimization to spiking neural networks,” _IEEE Signal Processing Magazine_ , vol. 36, no. 6, pp. 51–63, 2019.
  * [9] J. K. Eshraghian, “snntorch,” 2021. [Online]. Available: https://snntorch.readthedocs.io
  * [10] S. Sheik and M. Sorbaro, “sinabs,” 2013. [Online]. Available: https://sinabs.ai
  * [11] D. Muir, F. Bauer, and P. Weidel, “rockpool,” 2019. [Online]. Available: https://doi.org/10.5281/zenodo.4045345
  * [12] H. Hazan, D. J. Saunders, H. Khan, D. Patel, D. T. Sanghavi, H. T. Siegelmann, and R. Kozma, “Bindsnet: A machine learning-oriented spiking neural networks library in python,” _Frontiers in Neuroinformatics_ , p. 89, 2018.
  * [13] C. Pehle and J. E. Pedersen, “Norse - A deep learning library for spiking neural networks,” Jan. 2021, documentation: https://norse.ai/docs/. [Online]. Available: https://doi.org/10.5281/zenodo.4422025
  * [14] L. O. Chua, “Memristor-the missing circuit element,” _IEEE Transactions on Circuit Theory_ , vol. 18, no. 5, pp. 507–519, 1971.
  * [15] L. O. Chua and S. M. Kang, “Memristive devices and systems,” _Proceedings of the IEEE_ , vol. 64, no. 2, pp. 209–223, 1976.
  * [16] S. H. Jo, T. Chang, I. Ebong, B. B. Bhadviya, P. Mazumder, and W. Lu, “Nanoscale memristor device as synapse in neuromorphic systems,” _Nano Letters_ , vol. 10, no. 4, pp. 1297–1301, 2010.
  * [17] T. Serrano-Gotarredona, T. Masquelier, T. Prodromakis, G. Indiveri, and B. Linares-Barranco, “Stdp and stdp variations with memristors for spiking neuromorphic learning systems,” _Frontiers in Neuroscience_ , vol. 7, p. 2, 2013.
  * [18] C.-Y. Lin, J. Chen, P.-H. Chen, T.-C. Chang, Y. Wu, J. K. Eshraghian, J. Moon, S. Yoo, Y.-H. Wang, W.-C. Chen _et al._ , “Adaptive synaptic memory via lithium ion modulation in rram devices,” _Small_ , vol. 16, no. 42, p. 2003964, 2020.
  * [19] M. D. Pickett, G. Medeiros-Ribeiro, and R. S. Williams, “A scalable neuristor built with mott memristors,” _Nature Materials_ , vol. 12, no. 2, pp. 114–117, 2013.
  * [20] L. Bao, J. Zhu, Z. Yu, R. Jia, Q. Cai, Z. Wang, L. Xu, Y. Wu, Y. Yang, Y. Cai _et al._ , “Dual-gated mos2 neuristor for neuromorphic computing,” _ACS Applied Materials & Interfaces_, vol. 11, no. 44, pp. 41 482–41 489, 2019.
  * [21] J. Del Valle, P. Salev, Y. Kalcheim, and I. K. Schuller, “A caloritronics-based mott neuristor,” _Scientific Reports_ , vol. 10, no. 1, pp. 1–10, 2020.
  * [22] H. Lim, V. Kornijcuk, J. Y. Seok, S. K. Kim, I. Kim, C. S. Hwang, and D. S. Jeong, “Reliability of neuronal information conveyed by unreliable neuristor-based leaky integrate-and-fire neurons: a model study,” _Scientific Reports_ , vol. 5, no. 1, pp. 1–15, 2015.
  * [23] S. M. Kang, D. Choi, J. K. Eshraghian, P. Zhou, J. Kim, B.-S. Kong, X. Zhu, A. S. Demirkol, A. Ascoli, R. Tetzlaff, and L. O. Chua, “How to build a memristive integrate-and-fire model for spiking neuronal signal generation,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 2021.
  * [24] S. Hao, X. Ji, S. Zhong, K. Y. Pang, K. G. Lim, T. C. Chong, and R. Zhao, “A monolayer leaky integrate-and-fire neuron for 2d memristive neuromorphic networks,” _Advanced Electronic Materials_ , vol. 6, no. 4, p. 1901335, 2020\.
  * [25] X. Zhu, Q. Wang, and W. D. Lu, “Memristor networks for real-time neural activity analysis,” _Nature Communications_ , vol. 11, no. 1, pp. 1–9, 2020\.
  * [26] A. Loeffler, R. Zhu, J. Hochstetter, A. Diaz-Alvarez, T. Nakayama, J. M. Shine, and Z. Kuncic, “Modularity and multitasking in neuro-memristive reservoir networks,” _Neuromorphic Computing and Engineering_ , vol. 1, no. 1, p. 014003, 2021.
  * [27] R. Zhu, A. Loeffler, J. Hochstetter, A. Diaz-Alvarez, T. Nakayama, A. Stieg, J. Gimzewski, J. Lizier, and Z. Kuncic, “Mnist classification using neuromorphic nanowire networks,” in _International Conference on Neuromorphic Systems 2021_ , 2021, pp. 1–4.
  * [28] J. Hochstetter, R. Zhu, A. Loeffler, A. Diaz-Alvarez, T. Nakayama, and Z. Kuncic, “Avalanches and edge-of-chaos learning in neuromorphic nanowire networks,” _Nature Communications_ , vol. 12, no. 1, pp. 1–13, 2021.
  * [29] F. Cai, J. M. Correll, S. H. Lee, Y. Lim, V. Bothra, Z. Zhang, M. P. Flynn, and W. D. Lu, “A fully integrated reprogrammable memristor–cmos system for efficient multiply–accumulate operations,” _Nature Electronics_ , vol. 2, no. 7, pp. 290–299, 2019.
  * [30] J. K. Eshraghian, K. Cho, and S. M. Kang, “A 3-d reconfigurable rram crossbar inference engine,” in _2021 IEEE International Symposium on Circuits and Systems (ISCAS)_. IEEE, 2021, pp. 1–5.
  * [31] Y. Zhang, Z. Wang, J. Zhu, Y. Yang, M. Rao, W. Song, Y. Zhuo, X. Zhang, M. Cui, L. Shen _et al._ , “Brain-inspired computing with memristors: Challenges in devices, circuits, and systems,” _Applied Physics Reviews_ , vol. 7, no. 1, p. 011308, 2020.
  * [32] Y. Li, “Analog computing using 1t1r crossbar arrays,” 2018.
  * [33] J. K. Eshraghian, S.-M. Kang, S. Baek, G. Orchard, H. H.-C. Iu, and W. Lei, “Analog weights in reram dnn accelerators,” in _2019 IEEE International Conference on Artificial Intelligence Circuits and Systems (AICAS)_. IEEE, 2019, pp. 267–271.
  * [34] J. K. Eshraghian, X. Wang, and W. D. Lu, “Memristor-based binarized spiking neural networks: Challenges and applications.” _IEEE Nanotechnology Magazine_ , 2022.
  * [35] Z. Wang, S. Joshi, S. Savel’ev, W. Song, R. Midya, Y. Li, M. Rao, P. Yan, S. Asapu, Y. Zhuo _et al._ , “Fully memristive neural networks for pattern classification with unsupervised learning,” _Nature Electronics_ , vol. 1, no. 2, pp. 137–145, 2018.
  * [36] Z. Tang, Y. Chen, S. Ye, R. Hu, H. Wang, J. He, Q. Huang, and S. Chang, “Fully memristive spiking-neuron learning framework and its applications on pattern recognition and edge detection,” _Neurocomputing_ , vol. 403, pp. 80–87, 2020.
  * [37] F. Kiani, J. Yin, Z. Wang, J. J. Yang, and Q. Xia, “A fully hardware-based memristive multilayer neural network,” _Science Advances_ , vol. 7, no. 48, p. eabj4801, 2021.
  * [38] J. Wang, S. Hu, X. Zhan, Q. Yu, Z. Liu, T. P. Chen, Y. Yin, S. Hosaka, and Y. Liu, “Handwritten-digit recognition by hybrid convolutional neural network based on hfo2 memristive spiking-neuron,” _Scientific Reports_ , vol. 8, no. 1, pp. 1–7, 2018.
  * [39] P. Wijesinghe, A. Ankit, A. Sengupta, and K. Roy, “An all-memristor deep spiking neural computing system: A step toward realizing the low-power stochastic brain,” _IEEE Transactions on Emerging Topics in Computational Intelligence_ , vol. 2, no. 5, pp. 345–358, 2018.
  * [40] Q. Duan, Z. Jing, X. Zou, Y. Wang, K. Yang, T. Zhang, S. Wu, R. Huang, and Y. Yang, “Spiking neurons with spatiotemporal dynamics and gain modulation for monolithically integrated memristive neural networks,” _Nature Communications_ , vol. 11, no. 1, pp. 1–13, 2020.
  * [41] T. W. Molter and M. A. Nugent, “The generalized metastable switch memristor model,” in _CNNA 2016; 15th International Workshop on Cellular Nanoscale Networks and Their Applications_. VDE, 2016, pp. 1–2.
  * [42] C. Lammie, W. Xiang, B. Linares-Barranco, and M. R. Azghadi, “Memtorch: An open-source simulation framework for memristive deep learning systems,” _Neurocomputing_ , vol. 485, pp. 124–133, 2022.
  * [43] P.-Y. Chen, X. Peng, and S. Yu, “Neurosim: A circuit-level macro model for benchmarking neuro-inspired architectures in online learning,” _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ , vol. 37, no. 12, pp. 3067–3080, 2018.
  * [44] M. J. Rasch, D. Moreda, T. Gokmen, M. Le Gallo, F. Carta, C. Goldberg, K. El Maghraoui, A. Sebastian, and V. Narayanan, “A flexible and fast pytorch toolkit for simulating training and inference on analog crossbar arrays,” in _2021 IEEE 3rd International Conference on Artificial Intelligence Circuits and Systems (AICAS)_. IEEE, 2021, pp. 1–4.
  * [45] J. Huang, S. Stathopoulos, A. Serb, and T. Prodromakis, “Neuropack: An algorithm-level python-based simulator for memristor-empowered neuro-inspired computing,” _arXiv preprint arXiv:2201.03339_ , 2022.
  * [46] R. J. Williams and D. Zipser, “A learning algorithm for continually running fully recurrent neural networks,” _Neural Computation_ , vol. 1, no. 2, pp. 270–280, 1989.
  * [47] G. Bellec, F. Scherr, A. Subramoney, E. Hajek, D. Salaj, R. Legenstein, and W. Maass, “A solution to the learning dilemma for recurrent networks of spiking neurons,” _Nature Communications_ , vol. 11, no. 1, pp. 1–15, 2020\.
  * [48] Y. Demirag, C. Frenkel, M. Payvand, and G. Indiveri, “Online training of spiking recurrent neural networks with phase-change memory synapses,” _arXiv preprint arXiv:2108.01804_ , 2021.
  * [49] X.-Y. Wang, P.-F. Zhou, J. K. Eshraghian, C.-Y. Lin, H. H.-C. Iu, T.-C. Chang, and S.-M. Kang, “High-density memristor-cmos ternary logic family,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 68, no. 1, pp. 264 – 274, 2021.
  * [50] X.-Y. Wang, C.-T. Dong, P.-F. Zhou, S. K. Nandi, S. K. Nath, R. G. Elliman, H. H.-C. Iu, S.-M. Kang, and J. K. Eshraghian, “Low-variance memristor-based multi-level ternary combinational logic,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , 2022.
  * [51] J. K. Eshraghian, S. Baek, J.-H. Kim, N. Iannella, K. Cho, Y. S. Goo, H. H. Iu, S.-M. Kang, and K. Eshraghian, “Formulation and implementation of nonlinear integral equations to model neural dynamics within the vertebrate retina,” _International Journal of Neural Systems_ , vol. 28, no. 07, p. 1850004, 2018\.
  * [52] C. Lammie, J. K. Eshraghian, W. D. Lu, and M. R. Azghadi, “Memristive stochastic computing for deep learning parameter optimization,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 68, no. 5, pp. 1650–1654, 2021.
  * [53] L. Pentecost, A. Hankin, M. Donato, M. Hempstead, G.-Y. Wei, and D. Brooks, “Nvmexplorer: A framework for cross-stack comparisons of embedded non-volatile memories,” _arXiv preprint arXiv:2109.01188_ , 2021.
  * [54] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” _Advances in neural information processing systems_ , vol. 31, 2018.
  * [55] J. K. Eshraghian, M. Ward, E. Neftci, X. Wang, G. Lenz, G. Dwivedi, M. Bennamoun, D. S. Jeong, and W. D. Lu, “Training spiking neural networks using lessons from deep learning,” _arXiv preprint arXiv:2109.12894_ , 2021\.
  * [56] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, 1998.
  * [57] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms,” _arXiv preprint arXiv:1708.07747_ , 2017.
  * [58] A. Amir, B. Taba, D. Berg, T. Melano, J. McKinstry, C. Di Nolfo, T. Nayak, A. Andreopoulos, G. Garreau, M. Mendoza _et al._ , “A low power, fully event-based gesture recognition system,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 7243–7252.
  * [59] L. Patrick, C. Posch, and T. Delbruck, “A 128$\times$128 120 db 15$\mu$ s latency asynchronous temporal contrast vision sensor,” _IEEE Journal of Solid-State Circuits_ , vol. 43, pp. 566–576, 2008.
  * [60] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [61] Q. Wang, X. Wang, S. H. Lee, F.-H. Meng, and W. D. Lu, “A deep neural network accelerator based on tiled rram architecture,” in _2019 IEEE International Electron Devices Meeting (IEDM)_. IEEE, 2019, pp. 14–4.


# Bound orbits near black holes with scalar hair

I M Potashov  Ju V Tchemarina  and A N Tsirulev  Faculty of Mathematics, Tver
State University, 35 Sadovyi, Tver, Russia, 170002<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We consider spherically symmetric black holes with minimally coupled scalar
fields and concentrate our attention on asymptotically flat self-gravitating
configurations having the event horizons located at radii much smaller than
$2m$. We think of such configurations as rigorous mathematical models of the
gravitating objects, surrounded by dark matter, in the centres of normal
galaxies. It turns out that the radius of the event horizon of a scalar field
black hole always less than the Schwarzschild radius of vacuum black hole of
the same mass and can be arbitrary close to zero. In astronomical
observations, a key role in distinguishing between black holes, wormholes, and
naked singularities plays measuring parameters of bound quasiperiodic orbits,
in particular, the shape of an orbit and the angle of precession of its
pericentre. We consider a typical two-parameter family of compact scalar field
black holes and compute numerically the shapes of some bound orbits. We find
that a key feature of bound orbits around a compact black hole is that the
angle between closest pericentre points is either negative or, at least, less
than that for the Schwarzschild black hole of the same mass.

##  1 Introduction

One of the main problems of modern astrophysics is the identification of self-
gravitating objects in the centers of galaxies [1, 2, 3]. Namely, are they
black holes, wormholes, naked singularities, or maybe boson stars? In any of
these cases, the observations of bound orbits in the centers of galaxies play
a key role for this question [4, 5, 6, 7, 8, 9, 10]). In this paper, we study
some properties of bound orbits near black holes in the model of a static,
self-gravitating, spherically symmetric scalar field. The field is supposed to
be minimally coupled to gravity and the spacetime to be asymptotically flat.
This model is motivated by the modern astronomical observations of the centers
of galaxies. Strongly gravitating objects in galactic centers cannot be
regarded as being vacuum because they are surrounded by dark matter. In our
approach, we model dark matter by a self-interacting scalar field. It is very
probably that scalar fields exist in nature, but even in the opposite case, we
have a powerful degree of freedom to construct a phenomenological model of
dark matter. It is true, at least, for the spherically symmetric objects,
because the self-interaction potential can be chosen arbitrarily in the wide
class of physically admissible (particularly, positive at infinity [11])
potentials. In the other words, we can approximately model any spherically
symmetric distribution of matter in a self-consistent manner.

In fact, the form of the potential and the distribution of dark matter near
the centre should be found from the astronomical observations. Therefore, we
use the so-called inverse problem method [12, 13, 14, 15, 16] for static,
spherically symmetric, scalar field configurations. This method allows us to
deal with, in some sense, all the physically admissible potentials at the same
time. Our aim is to study some distinctive features of the geodesic motion of
massive test particles near the horizons of scalar field black holes. We
restrict our attention to compact black holes, which have the radius of the
event horizon much smaller than the corresponding Schwarzschild radius and,
consequently, show most distinctive features of the shape of bound orbits. In
Section 2 we describe the mathematical background of the inverse problem
method. Section 3 is devoted to an analytical formulation of a simple two-
parameter family of compact black holes. In Section 4 we compare the shapes of
bounded orbits around a typical compact black hole with those around the
Schwarzschild black hole of the same mass.

In this paper we adopt the metric signature $(+\,-\,-\,-)$ and use the
geometrical system of units with $G=1,\,c=1$.

##  2 Spherically symmetric scalar field black holes

Self-gravitating minimally coupled scalar fields are described by the action

$\Sigma=\frac{1}{8\pi}\int\\!\left(-\frac{1}{2}S-\langle
d\phi,d\phi\rangle-2V(\phi)\right)\sqrt[]{|g|}\,d^{\,4}x\,,$ (1)

where $S$ is the scalar curvature, $\phi$ the scalar field, $V(\phi)$ the
self-interaction potential, and the angle brackets denote the scalar product
with respect to the metric. We will write the metric of a static spherically
symmetric spacetime in the Schwarzschild-like coordinates (that is, $r$ is the
area coordinate) in the form

$ds^{2}=Adt^{2}-\frac{\,dr^{2}}{f}-r^{2}(d\theta^{2}+\sin^{2}\\!\theta\,d\varphi^{2}),$
(2)

where the metric functions $A$ and $f$, as well as the field $\phi$, depend
only on the radial coordinate. Any static spherically symmetric solution of
the Einstein-Klein-Gordon equations for the action (1) with an arbitrary
physically admissible potential $V(\phi)$ obeys the quadrature formulae [16]

$A(r)=2r^{2}\\!\\!\int\limits_{\\!r}^{\,\,\infty}\frac{\,\xi-3m}{\,r^{4}}\,\mathrm{e}^{F}dr\,,\quad
f(r)=\mathrm{e}^{-2F}A\,,$ (3)
$F(r)=-\\!\int\limits_{\\!r}^{\,\,\infty}\\!{\phi^{\prime}}^{2}rdr\,,\quad\xi(r)=r+\int\limits_{\\!r}^{\,\,\infty}\\!\\!\left(1-\mathrm{e}^{F}\right)\\!dr\,,$
(4)
$\widetilde{V}(r)=\frac{1}{2r^{2}}\\!\left(1-3f+r^{2}{\phi^{\prime}}^{2}\\!f+2\,\mathrm{e}^{-F}\,\frac{\,\xi-3m}{r}\right),$
(5)

where a prime denotes differentiation with respect to $r$,
$\widetilde{V}(r)=V(\phi(r))$, and the positive parameter $m$ is the
Schwarzschild mass. The ’inverse problem method’ consists of specifying a
monotonic function $\phi(r)$ and finding successively the functions
$\mathrm{e}^{F(r)}$, $\xi(r)$, $A(r)$, $f(r)$, $\widetilde{V}(r)$, and the
potential $V(\phi)=\widetilde{V}(r(\phi))$. The quadratures (3) – (5) give us,
in some sense, a general solution of the Einstein-Klein-Gordon equations.

We assume that $\phi(r)$ is of class $C^{2}\\!\left[(0,\infty)\right]$ and has
the asymptotic behaviour

$\phi=O\\!\left(r^{-1/2-\textstyle{\alpha}}\right)\\!,\;\;\,r\rightarrow\infty\;\;\;(\alpha>0).$
(6)

The functions $F(r)$, $\mathrm{e}^{F(r)}$, and $\xi(r)$ obey the obvious
relations

$F^{\prime}>0,\;\;\;\xi>r,\;\;\;0<\,\xi^{\prime}=\mathrm{e}^{F}\leqslant\,1,\;\;\;\xi^{\prime\prime}\\!=\left(\mathrm{e}^{F}\right)^{\prime}\geqslant
0\;\;\;\mbox{for all}\;\;r>0,$ (7)
$\mathrm{e}^{F}\\!=1+o\big{(}r^{-1}\big{)},\quad\xi\\!=r+o(1),\;\;\;r\rightarrow\infty,$
(8)
$\mathrm{e}^{F}=\mathrm{e}^{F(0)}+o(r),\;\;\;\xi\\!=\xi(0)+\xi^{\prime}(0)r+o(r),\;\;r\rightarrow
0,$ (9)

where

$\mathrm{e}^{F(0)}\equiv\lim\limits_{r\rightarrow
0}\mathrm{e}^{F(r)},\;\quad\xi(0)>0,\quad\;\;\xi^{\prime}(0)=\mathrm{e}^{F(0)}\geqslant
0.$ (10)

Note also that $\mathrm{e}^{F(0)}>0$ if and only if $\phi(r)$ grows slower
than $\ln{r}$ as $r\rightarrow 0$. Substituting the expansions (8) and (9)
into the quadrature (3) and using the conditions (7) and (10), we find that
$A(r)$ has the asymptotic behaviours

$A(r)=1-\frac{3m}{r}\,+\,o(1/r),\;\;r\rightarrow\infty,\qquad
A(r)=\frac{2}{3}\>\\!\frac{\xi(0)-3m}{r}\,\mathrm{e}^{F(0)}\,+\,O(1),\;\;r\rightarrow
0.$ (11)

It means, first, that $m$ is the Schwarzschild mass and, second, that this
metric function defines a black hole if and only if111It defines a naked
singularity if $3m<\xi(0)$, while if $3m=\xi(0)$, we obtain either a regular
solution or a naked singularity.

$3m>\xi(0).$ (12)

Furthermore, it follows directly from the quadrature (3) and the second
inequality in (7) that for a given scalar field the event horizon radius
$r_{\\!\scriptscriptstyle h}$ of a black hole with mass $m>\xi(0)/3$ is less
then the corresponding Schwarzschild radius, that is,
$r_{\\!\scriptscriptstyle{}h}<2m$.

In what follows, we will consider compact black holes for which
$3m\rightarrow\xi(0)+0$ and, consequently,
$r_{\\!\scriptscriptstyle{}h}\ll{}2m$ (see Fig 1).

##  3 A family of compact black holes

In this section we consider a typical two-parameter family of compact black
holes; the assertion about typicality is based on many numerical and
analytical experiments with other suitable examples of spherically symmetric
scalar fields. In the purely analytical treatment of the quadratures (3) –
(5), it is more convenient to start with some specially chosen function
$\mathrm{e}^{F(r)}$, because the functions $\mathrm{e}^{F(r)}$, $\xi(r)$, and
$\phi$ uniquely determine each other.

The strictly monotonic piecewise analytic function

$\mathrm{e}^{F(r)}=1-\alpha+\frac{\,2\alpha}{7}\\!\;r^{3},\;\;0\leqslant{}r\leqslant{}1,\quad\mbox{and}\quad\mathrm{e}^{F(r)}=1-\frac{\,2\alpha}{r^{3}}+\frac{\,9\alpha}{7r^{4}},\;\;1\leqslant{}r<\infty$
(13)

has continuous derivatives up to second order at the point $r=1$. By direct
integration in (4) and then in (3), we obtain

$\xi=\frac{\,3\alpha}{2}+(1-\alpha)\\!\;r+\frac{\,\alpha\,}{14}\\!\;r^{4},\;\;0\leqslant{}r\leqslant{}1,\quad\mbox{and}\quad\xi=r+\frac{\,\alpha\,}{r^{2}}-\frac{\,3\alpha}{7r^{3}}\,,\;\;1\leqslant{}r<\infty,$
(14) $A=-\frac{(1-\alpha)(2m-\alpha)}{r}+(1\\!-\\!\alpha)^{2}\,+\\\
\left(\frac{6}{7}\\!\;(2m-\alpha)\ln{}r+\frac{8}{5}-\frac{88}{105}\\!\;\alpha-\frac{54}{49}m\\!\right)\\!\alpha{}r^{2}-\frac{5}{7}\\!\;\alpha(1-\alpha)r^{3}-\frac{\alpha^{2}}{98}\\!\;r^{6},\;\;\;\;\;0\leqslant{}r\leqslant{}1,\qquad\;\;$
(15) $\quad
A=1\,-\,\frac{2m}{r}\,-\,\frac{2\alpha}{5r^{3}}\,+\,\left(\\!m+\frac{1}{7}\right)\\!\frac{2\alpha}{r^{4}}\,-\,\frac{54\alpha{}m}{49r^{5}}\,-\,\frac{\alpha^{2}}{2r^{6}}\,+\,\frac{10\alpha^{2}}{21r^{7}}\,-\,\frac{27\alpha^{2}}{245r^{8}}\,,\;\;1\\!\leqslant{}\\!r\\!<\\!\infty,\;\;$
(16)

where $\alpha\\!\in\\![0,1)$ is the parameter of ’intensity’ of the scalar
field. The field can be obtained by numerical solving of the problem
$\phi^{\prime}=\sqrt{F^{\prime}/r\,},\,\,\phi(\infty)=0$; at infinity
$\phi\sim\sqrt{8\alpha/3}\,r^{-3/2}$. Note that the length scale in (13) –
(16) is fixed by choosing the matching point at $r=1$. In studying bound
orbits, we are interested only in the metric functions $A$ and $f$ and can,
therefore, use an arbitrary unit of length. On the other hand, the solution
(3) – (5) is invariant under the scale transformations

$r\rightarrow r/a,\quad m\rightarrow m/a,\quad V\rightarrow
a^{2}V,\quad\;\;a>0.$

It means that by applying $a=m$ in this transformation, one can choose, as it
is usually done in general relativity, the mass of a black hole as the current
unit of length.

We see that the radius of the horizon $r_{\\!\scriptscriptstyle
h}\\!\rightarrow\\!0$ as $3m\rightarrow\xi(0)\\!+\\!0$ ($\xi(0)=3\alpha/2$)
or, equivalently, as $m\rightarrow\alpha/2\\!+\\!0$. For compact black holes,
one can estimate directly from (16) that $r_{\\!\scriptscriptstyle
h}\approx(2m-\alpha)/(1-\alpha)$. In what follows, we choose $\alpha=0.6$. For
$m=m_{\\!\scriptscriptstyle B\\!H}=\alpha/2+0.02=0.32$, we obtain a compact
black hole with $r_{\\!\scriptscriptstyle h}\approx 0.01$. The corresponding
metric function $A_{\scriptscriptstyle B\\!H}$, as well as the Schwarzschild
solution $A_{\scriptscriptstyle Sch}$ of the same mass and the naked
singularity solution $A_{\scriptscriptstyle N\\!S}$ with $\alpha=0.6$ and
$m=0.29$, are shown in Fig. 1. In the region $r>3r_{\\!\scriptscriptstyle
Sch}=1.92$ (outside the Schwarzshcild innermost stable circular orbits), where
we also have $\mathrm{e}^{F(r)}\\!\\!\>\big{|}_{\alpha=0.6}\approx 1$, all the
three spacetime geometries practically coincide.

Figure 1:  The functions $\xi(r)$ (left) and $A(r)$ (right). Parameters:
$r_{\\!*}=0.15$, $r_{\\!\scriptscriptstyle h}=0.01$, $\alpha=0.6$,
$m_{\\!\scriptscriptstyle BH}=0.32$, $m_{\\!\scriptscriptstyle Sch}=0.32$,
$m_{\\!\scriptscriptstyle NS}=0.29$.

##  4 Bound orbits near the event horizons of compact black holes

Figure 2: Top-left panel: $\alpha=0.6$, $m=\alpha/2+0.02=0.32$,
$\Delta\varphi=-\pi$, $r_{\\!\scriptscriptstyle min}=0.31$,
$r_{\\!\scriptscriptstyle max}=1.62$, $J=0.6$, $E^{2}=0.689$. In the other
plots, $\alpha=0.6$, $r_{\\!\scriptscriptstyle max}\approx 8.04$. Top-right
panel: $m=\alpha/2+10^{-4}=0.3001$, $r_{\\!\scriptscriptstyle h}=5\cdot
10^{-4}$, $\Delta\varphi=-2\pi/3$, $J=0.262$, $E^{2}=0.926$. Bottom-left
panel: $m=\alpha/2+10^{-5}=0.30001$, $r_{\\!\scriptscriptstyle h}=5\cdot
10^{-5}$, $\Delta\varphi=-3.05$, $J=0.02$, $E^{2}=0.925$. Bottom-right panel:
the Schwarzschild black hole, $m=0.3001$, $\Delta\varphi=+3.01$, $J=1.2$,
$E^{2}=0.946$.

In a spherically symmetric spacetime, the Lagrangian for geodesics is constant
along each geodesic and does not depend explicitly on the coordinates $t$ and
$\varphi$. We write the corresponding integrals of motion as

$\frac{dt}{ds}=\frac{E}{A}\,,\qquad\frac{d\varphi}{ds}=\frac{J}{\,r^{2}}\,,\qquad\left(\frac{dr}{ds}\right)^{\\!2}=\mathrm{e}^{-2F}\\!\left(E^{2}-{V\vphantom{\underline{A}}}_{\\!\\!eff}\right),\quad{V\vphantom{\underline{A}}}_{\\!\\!eff}=A\\!\left(1+\frac{\,J^{2}}{\,r^{2}}\right),$
(17)

where $E$ is the specific energy and $J$ the specific angular momentum of a
massive free particle. The effective potential of a black hole spacetime
vanishes at the horizon, tends to unity as $r\rightarrow\infty$, and has, for
sufficiently large $J$, at least one minimum and one maximum outside the
horizon. We are primarily interested in the shape of a bound orbit and the
angle of precession $\Delta\varphi$ of the orbit, which can be expressed by
the obvious relations

$\varphi_{\\!{}_{\scriptscriptstyle
osc}}=2J\\!\int\limits_{\,r_{\\!\scriptscriptstyle
min}}^{\>r_{\\!\scriptscriptstyle
max}}\\!\frac{\,\mathrm{e}^{F}}{\,r^{2}\sqrt{E^{2}-V_{\\!\scriptscriptstyle
eff}\,}}\,dr\,,\quad\Delta\varphi\,=\,\varphi_{\\!{}_{\scriptscriptstyle
osc}}-2\pi\,,$ (18)

where $r_{\\!\scriptscriptstyle min}$ and $r_{\\!\scriptscriptstyle max}$ are
solutions of the equation $E^{2}-V\vphantom{\underline{A}}_{\\!\\!eff}=0$
which are uniquely determined by the requirement of being located most nearest
to (and on the opposite sides of) the minimum of the effective potential.
Thus, a bound orbit of the general type oscillates near a stable circular
orbit and $\varphi_{\\!{}_{\scriptscriptstyle osc}}$ is the angle between two
successive oscillations.

There is a key difference between bound orbits around a compact black holes
and the Schwarzschild black hole of the same mass: in the former case, the
angle between closest pericentre points is either negative or, at least,
(much) less than that in the latter case. Note that a bound orbit around a
vacuum black hole always has the advance of its pericentre. Fig. 2 represents
typical trajectories of a test particle. An orbit with the negative angle of
precession $\Delta\varphi=-\pi$ is shown in the top-left panel. Two prolate
bound orbits (with one and the same apocentre radius) around an extremely
compact black hole are represented in the top-right and bottom-left panels. A
bound orbit (with approximately the same apocentre radius) around the
Schwarzschild black hole of the same mass is shown in the bottom-right panel.
It should be stressed that these shapes of bound orbits and their precession
angle deficits are typical for compact black holes.

##  5 Conclusions

We consider compact spherically symmetric black holes with minimally coupled
scalar fields and concentrate our attention on asymptotically flat self-
gravitating configurations having the event horizons located at radii much
smaller than $2m$. Such configurations can be thought of as rigorous
mathematical models of the strongly gravitating objects, surrounded by dark
matter, in the centres of normal galaxies. We show that the radius of the
event horizon of a scalar field black hole always less (possibly much less)
than the Schwarzschild radius of vacuum black hole of the same mass and can be
arbitrary close to zero. We also consider bound orbits around compact black
holes and find that the angle between closest pericentre points of such orbits
is either negative or, at least, less than that for the Schwarzschild black
hole of the same mass. The numerical experiments show that the corresponding
radius of the innermost stable circular orbit is also less than in the vacuum
case. These features play key roles in distinguishing between scalar field
black holes, vacuum black holes, wormholes, naked singularities, and boson
stars. Note that boson stars have no innermost stable circular orbits, and
scalar field wormholes and naked singularities (of the same positive mass)
have marginal (or degenerated) stable circular orbits with the zero angular
momentum of test particles. These result can possibly be applied to
interpretations of future astronomical observations. In particular, one can
hope that direct observation of the central region of $SgrA^{*}$, will soon be
obtained by the Event Horizon Telescope Collaboration.

## References

## References

  * [1] Akiyama K et al 2015 Astrophys. J. 820, 90 (arXiv: 1602.05527)
  * [2] Fish V L et al 2016 Astrophys. J. 820, 150 (arXiv: 1505.03545)
  * [3] Goddi G et al 2017 Int. J. Mod. Phys. D20, 1730001-239 (arXiv: 1606.088879)
  * [4] Li Z and Bambi C 2014 Phys. Rev. D 90, 024071 (arXiv: 1405.1883)
  * [5] Vieira R S S, Schee J, Kluźniak W, Stuchlík Z, and Abramowicz M 2014 Phys. Rev. D 90, 024035
  * [6] De Laurentis M, Younsi Z, Porth O, Mizuno Y, and Rezzolla L 2018 Phys. Rev. D 97, 104024
  * [7] Stashko O S and Zhdanov V I 2018 Gen. Relat. Gravit. 50, Issue 5, 105 (arXiv: 1702.02800)
  * [8] Mishra A and Chakraborty S 2018 Eur. Phys. J. C 78, Issue 5, 374 (arXiv: 1710.06791)
  * [9] Willenborg F, Grunau S, Kleihaus B, and Kunz J 2018 Phys. Rev. D 97 124002 (arXiv: 1801.09769)
  * [10] Kratovitch P V, Potashov I M, Tchemarina Ju V, and Tsirulev A N 2017 Journal of Physics: Conference Series 934, Issue 1, 012047 (arXiv: 1805.04447)
  * [11] Nikonov V V, Tchemarina Ju V and Tsirulev A N 2008 Class. Quantum Grav. 25 138001
  * [12] Bechmann O and Lechtenfeld O 1995 Class. Quantum Grav. 12 1473 (arXiv: gr-qc/9502011)
  * [13] Bronnikov K A and Shikin G N 2002 Grav. Cosmol. 8 107 (arXiv: gr-qc/0109027)
  * [14] Tchemarina Ju V and Tsirulev A N 2009 Grav. Cosmol. 15 94
  * [15] Azreg-Aïnou M 2010 Gen. Rel. Grav. 42 1427 (arXiv: gr-qc/0912.1722)
  * [16] Solovyev D A and Tsirulev A N 2012 Class. Quantum Grav. 29 055013