# Text Classification in Memristor-based Spiking Neural Networks

Jinqi Huang∗1, Alex Serb2, Spyros Stathopoulos1, Themis Prodromakis2 1
Electronics and Computer Science, University of Southampton  2 School of
Engineering, University of Edinburgh<EMAIL_ADDRESS>
###### Abstract

Memristors, emerging non-volatile memory devices, have shown promising
potential in neuromorphic hardware designs, especially in spiking neural
network (SNN) hardware implementation. Memristor-based SNNs have been
successfully applied in a wide range of applications, including image
classification and pattern recognition. However, implementing memristor-based
SNNs in text classification is still under exploration. One of the main
reasons is that training memristor-based SNNs for text classification is
costly due to the lack of efficient learning rules and memristor non-
idealities. To address these issues and accelerate the research of exploring
memristor-based spiking neural networks in text classification applications,
we develop a simulation framework with a virtual memristor array using an
empirical memristor model. We use this framework to demonstrate a sentiment
analysis task in the IMDB movie reviews dataset. We take two approaches to
obtain trained spiking neural networks with memristor models: 1) by converting
a pre-trained artificial neural network (ANN) to a memristor-based SNN, or 2)
by training a memristor-based SNN directly. These two approaches can be
applied in two scenarios: offline classification and online training. We
achieve the classification accuracy of 85.88% by converting a pre-trained ANN
to a memristor-based SNN and 84.86% by training the memristor-based SNN
directly, given that the baseline training accuracy of the equivalent ANN is
86.02%. We conclude that it is possible to achieve similar classification
accuracy in simulation from ANNs to SNNs and from non-memristive synapses to
data-driven memristive synapses. We also investigate how global parameters
such as spike train length, the read noise, and the weight updating stop
conditions affect the neural networks in both approaches. This investigation
further indicates that the simulation using statistic memristor models in the
two approaches presented by this paper can assist the exploration of
memristor-based SNNs in natural language processing tasks.

Keywords: Text classification, sentiment analysis, natural language
processing, spiking neural networks, memristors, online learning, offline
classification

##  1 Introduction

SNNs [1], as the third generation of ANNs [2], have become a widely-researched
topic recently. Their spike-based nature and temporal information coding
scheme bring advantages in biological plausibility [3, 4], sparsity [5], and
stochasticity [6], further making them power-efficient tools compared with
traditional ANNs. Therefore, they have been widely used in applications
spanning across different areas, such as hand gesture detection [7, 8], gait
detection [9], financial time series prediction [10], music composer
classification [11], real-time signal processing [12], disease detection [13],
and robotics [14, 15, 16].

Despite the promising potential, hardware implementations of SNNs in
conventional Von-Neumann processors can hardly compare with the biological
systems regarding energy efficiency. This performance gap mainly stems from
the frequent data movements between processing units and the memory to fetch
SNN-specific internal variables such as membrane voltages [17, 18]. In-memory
computing allows an alternative approach to implementing SNNs in hardware. The
computation is executed locally without frequent data movements by combining
the processing units and the memory. This working scheme removes the
restriction set by the limited communication bandwidth that commonly exists in
Von-Neumann processors. An SRAM-based in-memory computing system [19] has
reported having the energy efficiency 10x better than google’s TPU [20]. Based
on this idea, many dedicated SNN hardware designs have been developed, from
fully digital [21, 22] to mixed-signal [23], from large-scaled computation
platforms [24] to systems-on-wafer [25]. Some designs have been reported to
achieve biological speed [22] or multiple orders of magnitude superior to
their conventional Von-Neumann predecessors [21].

Meanwhile, the development of emerging non-volatile technology has also seen
growth in the past few years. Memristors are known as one of the most mature
emerging resistive memory devices among all. The extreme downscaled size [26],
the large-scale integrability [27], multi-bit storage [28], crossbar
arrangement [29], ns-scale switching speed [30] and high power efficiency [31]
make memristors good candidates to replace SRAM in in-memory computing
applications [32, 33, 34]. Moreover, memristors behave similarly to biological
synapses in neural networks resulting from their inherent physical properties
inspiring their use in spiking neural networks. As examples, References [35,
36] incorporated memristors as synapses in spiking neural networks to perform
unsupervised learning using spike-timing-dependent plasticity (STDP) [37].
References [38, 39, 40, 41] demonstrated pattern recognition and image
classification in memristor-based spiking neural networks.

However, only a few works reported memristor-based spiking neural networks for
text processing tasks: Design [42] employed phase-change memory (PCM), a type
of memristor device, to accelerate spiking neural network for language
modelling, and Design [43] performed a pattern generation task in recurrent
spiking neural networks with PCM models using e-prop learning rule [44]. One
of the main challenges is the costly training of memristor-based spiking
neural networks for text processing tasks. Firstly, gradient-based learning
rules, regardless of local (e.g. e-prop) or non-local (e.g. surrogate gradient
descent [45, 46]), require averaging the accumulated errors over the spike
train window that is used to represent a single continuous value. This method
considers the effects of every single spike when updating weights. It is
inefficient regarding the computational speed and the space efficiency,
especially when the length of spike trains for representing a single numerical
value is very long and when memristors are involved in the design. Besides,
compared to one-hot vectors, word embeddings, such as GloVe [47] and Word2vec
[48], are more commonly used to represent words in text classification tasks,
to improve space efficiency and to map the relations between words. In
conventional ANNs, a word embedding layer serves as a lookup table to obtain
the dense representation of the input texts. Pre-trained word embeddings are
usually fine-tuned later with other network parameters in training to improve
classification accuracy. During the training process, they are treated as
linear layers using backpropagation (BP) [49]. However, in SNNs, the dense
word representation needs to be converted to spike trains as the input to the
network [50, 51]. So far, there is no theoretical support regarding how to
train word embeddings in spiking neural networks. Lastly, reading and writing
memristor arrays can also add variations due to the read noise and the write
variation. One solution to these issues is to run the training in a software
simulator with realistic memristor models before starting the hardware design.
Compared with training in hardware, this solution offers a fast data process
and intuition for the system’s performance. Reference [52] has introduced an
algorithm-level simulator for memristor-based spiking neural networks in
memristor models, and Design [38] also used this solution to simulate
recurrent spiking neural networks in PCM models.

In this paper, we develop a simulation framework in Pytorch [53] with an
empirical memristor model proposed by [54] to perform text classification
tasks. To use this framework, we firstly present two approaches to obtaining
trained memristor-based spiking neural networks for performing text
classification tasks. One is converting a pre-trained ANN to its equivalent
SNN, and another is training an SNN directly using an ANN-based learning rule.
Both approaches train the network using gradient descent-based learning rules
without considering the effect of every single spike to improve the
computational speed and the space efficiency. In this way, we can also treat a
word embedding layer as a linear layer in a conventional ANN and train it
using BP. Then, we showcase a sentiment analysis task in the IMDB movie
reviews dataset [55] to validate both approaches. To the best of our
knowledge, this is the first demonstration of a text classification task
performed in a spiking neural network with a realistic memristor model.
Lastly, we summarise the classification accuracy from both approaches and
investigate how global parameters affect the system performance in both
approaches.

The contributions of this work are summarised below:

  1. 1.

Developing a simulation framework with an empirical memristor model to
demonstrate the first text classification task in spiking neural networks with
a realistic memristor model.

  2. 2.

Presenting and benchmarking two approaches to obtain trained memristor-based
spiking neural networks for text classification tasks.

  3. 3.

Investigating system sensitivity to global parameters.

To achieve these goals, we first introduce the overview of the two approaches
in section 2.1 and explain the details in section 2.2 and 2.3. Following are
the method to add the memristor model in the simulation framework in section
2.4 and the weight updating scheme in section 2.5. Section 3 displays all
experimental results for the text classification task in the simulation
framework with memristor models, and section 4 summarises the paper.

Figure 1: Two approaches to obtaining trained SNNs: approach 1 trains ANNs
first (coloured orange), converts trained ANNs to equivalent SNNs (coloured
green), and then maps SNN weights to memristors (coloured red); approach 2
converts ANNs to SNNs (coloured green), adds a memristor model (coloured red),
and then trains memristor-based SNNs directly using ANN-based learning rules
(coloured orange).

##  2 Methods

###  2.1 Methodology Overview

The main difference between ANNs and SNNs is the conveyance method: ANNs
employ continuous values, whereas SNNs use 0-or-1 spikes. As a result, the
basic idea to bridge the gap between two neural network architectures is to
find the relations to map continuous values to spikes. Reference [56]
suggested that the spike rates of an SNN are proportional to the ReLU [57]
activation outputs of the ANN equivalents with an error term that is ignorable
in shallow networks. Given that inputs and weights are restricted within [0,
1] in a memristor-based SNN, the ReLU activation outputs are always non-
negative. Therefore, the spike rates of an SNN are proportional to the input
currents. Reference [58] further proved that for neurons reset by subtraction,
the larger the membrane voltage threshold, the smaller the error term.
Additionally, weight normalisation [59] and adaptive threshold [60] are also
commonly used in ANN-to-SNN conversion to match the accuracy. These past works
lay the theoretical foundation of the two approaches to obtaining trained
memristor-based SNNs proposed in this paper.

Figure 1 gives the diagram of the two approaches: Starting from an untrained
ANN, approach 1 trains the ANN using typical ANN-based learning rules, then
converts the trained ANN to a trained SNN; approach 2 firstly converts the ANN
to its equivalent SNN before training with ANN-based learning rules. After
obtaining the trained SNN, memristor models are added to store the weights. In
approach 1, this step is as simple as mapping the weights in memristor arrays
as conductance; in approach 2, memristor resistance states (RSs) are trained
alongside the whole SNN. Therefore, the evolution of RSs can be observed
during the training process when using approach 2.

Figure 2: Workflow of the simulation framework for performing the sentiment
analysis task in the IMDB review dataset. Inputs, ANN structure, SNN-specified
structure, and memristor-related modules are coloured blue, orange, green and
red, respectively. $b$, $s$, $e$, and $o$ represent batch size, sentence
length, word embedding dimensions, and output dimensions, respectively. We use
a movie review, ’A wonderful little production’, as an example to explain the
input pre-process procedure: firstly, input sentences are tokenised to lists
of words (’[’A’, ’wonderful’, ’little’, ’production’]’ in the example).
Secondly, the word lists are converted to lists of word IDs (’[6, 385, 120,
370]’ in the example). Sentences in the same batch are padded to have the same
length. The lists of words are further transformed to a one-hot vector, and
padded one-hot vectors of the sentences in the same batch are packed into an
input representation tensor with the size of ($b\times s\times v$). Notably,
in Pytorch, the word ID lists do not need to be converted to one-hot
representation.

Now we walk through the workflow of the designed simulation framework
performing the sentiment analysis task to explain the two approaches. Figure 2
shows the network we employed in the sentiment analysis task in the IMDB
reviews dataset. After being tokenised using the basic English tokeniser and
transformed to word IDs, sentences in the same batch are padded to have the
same sentence length and packed into an input one-hot representation tensor
with the size of ($b\times s\times v$) (where $b$, $s$, $v$ are batch size,
sentence length, and vocabulary size, respectively). The one-hot
representation tensor is then fed to a word embedding layer to generate a
dense word representation whose size is ($b\times s\times e$) (where $e$ is
the word embedding dimension), followed by an average pooling layer to obtain
an averaged sentence representation across the whole sentence by squeezing the
sentence dimension.

In approach 1, the averaged sentence representation with the size of ($b\times
e$) is sent straight to a linear layer whose weights’ size is ($e\times o$)
(where $o$ is the output dimension). Specifically, there is only one output
neuron in the sentiment analysis task because we only have two categories:
’positive’ or ’negative’. The binary cross-entropy loss is calculated once the
inference results are obtained, and the errors are back-propagated to update
weights in the linear and word embedding layers. After the network is trained,
the conversion to an SNN starts. Firstly, the averaged sentence representation
is transformed to Poisson spike trains which are delivered to the single-layer
spiking neural network afterwards. We use the leaky integrate-and-fire neuron
model [61] in this task. Next, the weights in the SNN are mapped in memristor
arrays. When executing an inference, weights are loaded from the memristor
arrays. The output neuron membrane voltages are calculated according to the
mathematical expression of the neuron model, and the firing states are
determined by comparing the voltage with the threshold. The spike rates within
a time window give the final inference results: if the spike rate is higher
than 50%, it yields a ’positive’ result; otherwise, a ’negative’ result is
given.

In approach 2, the averaged sentence representation is likewise transformed to
Poisson spike trains, and memristor arrays are randomly initialised. After an
inference completes, the errors are back-propagated according to the same ANN-
based learning rule used in approach 1. The weight changes are converted to
the memristive RSs to be used to trigger RS updates. The errors are further
back-propagated to the word embedding layer to generate new spike trains in
the following iterations.

In the next two subsections, more details that need to be taken into
consideration in implementing two approaches in the simulation framework will
be introduced.

###  2.2 Approach 1: Converting a Pre-trained ANN to an SNN

Weights in the linear and word embedding layers need to be constrained in [0,
1] when training the ANN using approach 1. With this restriction, when
converting the ANN to the SNN, weights can be represented according to the
following equation:

$w=\frac{G-G_{min}}{G_{max}-G_{min}}$ (1)

Where $G$, $G_{min}$, $G_{max}$ are the current, the minimum, and the maximum
memristor conductance, respectively. The averaged sentence representation is
also restricted within [0, 1]. Therefore, when converting the ANN to the SNN,
spike trains can be generated using the continuous values in the averaged
sentence representation as firing rates. Reference [50] proposed a method to
determine whether to fire a spike at a certain time step to generate a Poisson
spike train using the equation shown below:

$P_{t}=\cases{1&if $x_{c}>x_{random,t}$\\\ 0&otherwise\\\ },\quad 0\leq t\leq
T$ (2)

Where $P_{t}$ gives the spike generation results of the averaged sentence
representation in the given time step $t$ in the spike train that uses $T$
time steps to represent a single continuous value, $x_{c}$ is the continuous
value in the averaged sentence representation, and $x_{random,t}$ is the
random value generated in the time step $t$ within [0, 1]. A matrix of random
values with the size of ($e\times T$) is initialised before the binary word
embedding conversion. The averaged sentence representation is then compared to
the random-value matrix to generate Poisson spike trains stored in a matrix
with the size of ($e\times T$) as inputs of the single-layer SNN.

Because the weights in the linear layer $W_{s}$ and the averaged sentence
representation $x_{c}$ are both non-negative when training the ANN, the
outputs of the linear layer calculated by the equation $y=\sum W_{s}x_{c}$ are
never negative. A negative constant offset $C$ need to be added to the ANN
outputs when using the sigmoid function to map the continuous-valued outputs
to probability. If the negative offset $C$ is not used, the non-negative
sigmoid function input will result in an output probability never smaller than
50%. Thus, to make the probability mapping functional, the term passed to the
sigmoid function becomes $(y+C)$ instead of $y$. When the term $(y+C)$ equals
0, the inference result is right on the boundary of dividing two categories.
Therefore, $C=-\sum\bar{W_{s}}\bar{x_{c}}=-0.5\times 0.5\times e$ where
$\bar{W_{s}}$ and $\bar{x_{c}}$ are both the theoretical averaged values (0.5
and 0.5, respectively). However, when converted to the SNN, the criterion for
dividing categories is the output spike rates rather than the sigmoid output.
Therefore, a sigmoid function is not used in SNN inference, and no offset is
needed in the converted SNN.

The membrane voltage threshold is an extra parameter not used in an ANN. The
threshold value needs to be carefully chosen to achieve comparable
classification accuracy in the converted SNN. Membrane voltages that are
either much higher or much lower than the threshold can cause information loss
[60]. A threshold value ’in the middle’ can potentially cause the least
information loss. In converted SNNs, the membrane voltages range is known
since both weights and word embeddings are constrained within [0, 1]. When all
inputs and weights are equal to ’1’ along the whole spike train, this gives
the highest theoretical firing rate. As explained in [58], the output firing
rates of an SNN are proportional to the inputs. Therefore, when we keep the
weights equal to ’1’ and change the inputs to half ’1’ and half ’0’, the
firing rate of the output neuron is supposed to be 50%—setting the membrane
voltage of the output neuron in this scenario as the threshold with the value
of $0.5\times e$ is a good starting point. It guarantees maximum firing when
the inputs are maximised and that under any initial conditions, the output
neuron will losslessly accumulate membrane potential in a manner commensurate
to the strength of the inputs without saturation or clipping. The threshold
value needs further fine-tuning based on the approximated value to achieve
better performance.

###  2.3 Approach 2: Training an SNN directly

Training SNNs with ANN-based gradient descent learning rules [62] can be
problematic due to the non-differential property of membrane voltages. A
common solution is to employ surrogate derivatives [45, 46] by using straight-
through estimators [63, 64, 65]. This work proposed an alternative solution,
which stems from ANN-based gradient-based algorithms, to avoid directly
dealing with non-differentiable functions. In this work, we use the leaky
integrate-and-fire (LIF) neuron model [61] with reset by subtraction
mechanism:

$\displaystyle V_{t}=V_{t-1}+\sum W_{s}x_{t}-Z_{t-1}V_{th}$ (3) $\displaystyle
Z_{t}=h(V_{t}-V_{th})$ (4)

Where $V_{t}$, $x_{t}$, and $Z_{t}$ denote membrane voltages, input spikes and
output spikes at time step $t$, respectively. $V_{th}$ is the membrane voltage
threshold, $W_{s}$ represents the weights in the single-layer SNN, and $h(x)$
is the Heaviside step function. For this two-category classification task, we
use the binary cross-entropy cost function with a sigmoid layer:

$E=-(\hat{y}\log y+(1-\hat{y})(1-\log y))$ (5)

Where $E$ represents the binary cross-entropy loss, $\hat{y}$ and $y$ are
labels and the probabilities given by the network, respectively.

As the learning rule, we choose Adaptive Gradient Algorithm (Adagrad) [66], a
gradient-based learning rule with an adaptive learning rate incorporating the
past gradients. The learning rule is explained by the equations below:

$\displaystyle s\leftarrow s+g^{2}$ (6)
$\displaystyle\theta\leftarrow\theta-\eta\frac{g}{\sqrt{s}+\epsilon}$ (7)

Where $g$ is the gradient of the objective function over some variable of
interest $\theta$, $s$ is a state variable that accumulates the square of the
gradients, and $\epsilon$ is a small value to avoid dividing by 0.

Now we walk through the steps to acquire the final expression of the weight
updates. For the same reason explained in subsection 2.2, the output neuron
firing rate is always non-negative. Therefore, the output firing rate needs to
add a negative offset before being passed to a sigmoid function. The offset
chosen here is -0.5. The output probabilities of two categories used in
calculating the loss are given by passing the output firing rate with a
negative offset to a sigmoid function:

$\displaystyle a=r+C$ (8) $\displaystyle y=\frac{1}{1+e^{-a}}$ (9)

Where $r$ is the output neuron firing rate, $C$ is the negative offset, and
$a$ is a temporary variable to save the sum of the firing rate and the offset.
The output firing rates of an SNN are proportional to the equivalent ANN
activation outputs in the equivalent ANN, as proven in [58]. Given that the
outputs are always non-negative, we further proved that in a memristor-based
SNN, the output firing rates are proportional to the equivalent ANN outputs as
follows (please see the derivation in supplementary material section 1):

$r=\frac{V_{c}}{V_{th}}-\frac{V_{t}-V_{0}}{V_{th}\cdot T}$ (10)

Suppose the equivalent ANN has the mathematical expression shown below:

$V_{c}=\sum W_{s}x_{c}$ (11)

Where $V_{c}$ and $x_{c}$ are continuous-valued ANN outputs and inputs, and
$V_{0}$ is the initial membrane voltage of a spiking neuron when processing a
new sample. $V_{0}$ is 0 in this work because we reset the membrane voltage
accumulator before a new inference. The error term $V_{t}/(V_{th}\cdot T)$ in
Equation 10 depends on the final membrane voltage $V_{t}$ at the end of the
spike window $T$, and $V_{c}$ is the continuous value that the spike train in
the spike window $T$ represents. Intuitively speaking, $V_{t}$ implicitly
depends on $V_{c}$. Therefore, the error term $V_{t}/(V_{th}\cdot T)$ is
weakly correlated to $V_{c}$, and this weak correlation can potentially
introduce an error in the conversion of an ANN into an SNN. This error can be
accumulated in a deep neural network, but it is ignorable in a shallow network
[58]. If we ignore the weak dependency of the error term to $V_{c}$ and use
two constants $\alpha$ and $\beta$ to replace the error term and the constant-
coefficient $1/V_{th}$, we can describe the output firing rates differently:

$r=\alpha V_{c}+\beta$ (12)

Therefore, the derivative of $r_{t}$ with respect to $V_{c}$ equals to
$dr/dV_{c}=\alpha$. As a result, the gradient of the cost function over the
weights in the single-layer SNN is expressed as follows using the chain rule:

$\eqalign{\frac{dE}{dW_{s}}&=\frac{dE}{da}\cdot\frac{da}{dV_{c}}\cdot\frac{dV_{c}}{dW_{s}}\cr&=(y-\hat{y})\cdot\alpha\cdot
x_{c}}$ (13)

Similarly, the word embedding matrix can be updated by further back-
propagating the errors:

$\eqalign{\frac{dE}{dW_{e}}&=\frac{dE}{da}\cdot\frac{da}{dV_{c}}\cdot\frac{dV_{c}}{dx_{c}}\cdot\frac{dx_{c}}{dW_{e}}\cr&=(y-\hat{y})\cdot\alpha\cdot
W_{s}x_{e}}$ (14)

Where $W_{s}$ denotes the parameters in the word embedding matrix, and $x_{e}$
is the one-hot input word representation.

In Adagrad, the change amounts are decided by the term
$g/(\sqrt{s}+\epsilon)$. Both $g$ and $\sqrt{s}$ are proportional to $\alpha$.
Therefore, $\alpha$ can be cancelled out in the expression. The final
expression to be used in Adagrad without $\alpha$ is simplified as follows:

$\displaystyle g_{s}=(y-\hat{y})\cdot x_{c}$ (15) $\displaystyle
g_{e}=(y-\hat{y})\cdot W_{s}x_{e}$ (16)

Compared to accumulating errors through the whole spike train using surrogate
gradient descent, which can be expressed as the equation below, our method
improves the computational speed and the space efficiency by skipping the step
of accumulating errors across the spike window.

$\eqalign{\frac{dE}{dW_{s}}&=\frac{dE}{da}\cdot\frac{da}{dZ_{t}}\cdot\frac{dZ_{t}}{dy_{t}}\cdot\frac{dV_{t}}{dW_{s}}\cr&=(y-\hat{y})\cdot\frac{1}{T}\sum^{T}h^{\prime}(V_{t}-V_{th})x_{t}}$
(17)

If using a binary activation [65] shown in Equation 18 to replace the
Heaviside step function and choosing a threshold equal or larger than the half
of the maximum possible membrane voltage to avoid information loss, the
gradient calculated by surrogate gradient descent (Equation 17) can be
simplified to $(y-\hat{y})\cdot x_{c}/(2V_{th})$. Cancelling out the constant
coefficient $1/(2V_{th})$ gives the same simplified gradient expressions shown
in Equation 15 and 16.

$h^{\prime}(V_{t}-V_{th})=\cases{\frac{1}{2V_{th}}&if $0<V_{t}<2V_{th}$\\\
0&otherwise\\\ }$ (18)

###  2.4 Adding Memristor Models

This work employs the statistical memristor model proposed by [54] to predict
the behaviour of memristors during the training or inference process in
memristor-based SNNs. The switching dynamics of memristor devices are shown
below for convenience:

$\frac{dR}{dt}=\cases{A_{p}(-1+\exp(|v|/{t_{p}}))h(r_{p}(v)-R)(r_{p}(v)-R)^{2}&if
$v>0$\\\ A_{n}(-1+\exp(|v|/{t_{n}}))h(R-r_{n}(v))(R-r_{n}(v))^{2}&otherwise\\\
}$ (19)

Where $A_{p,n}$ denotes the scaling factor, the term
$(-1+\exp(|v|/{t_{p,n}}))$ reflects the exponential increment between the
switching rate and the bias voltage $v$, and the rest of the equation portrays
the dependency on the current memristive RSs and the dynamic upper/lower
boundaries $r_{p,n}$ represented as the equation below with fitting parameters
$a_{0p,n}$ and $a_{1p,n}$:

$r_{p,n}=a_{0p,n}+a_{1p,n}v$ (20)

In practice, all parameters can be extracted using the method proposed by [67]
with the aid of characterisation tools [68, 69]. This work takes these
parameters to predict the memristive RSs given the current memristive RSs and
the bias voltages. Using Pytorch, all memristive RSs are updated and accessed
simultaneously. Given the pulse voltage v and the pulse-width pw, the duration
is converted to iteration times given by
$\mathbf{N}=\mathbf{int}(\mathbf{pw}/dt)$ and memristive RSs are updated in
loops with equation $\Delta R=\sum\limits^{N}(dR/dt)$.

###  2.5 Weight Updating Scheme

After obtaining the expected weights to map in memristor arrays, the essential
step is to convert the expected weights to expected memristive RSs. Rewriting
Equation 2 gives the expression of the expected memristive RSs given the
expected weights:

$R=\frac{1}{w(\frac{1}{R_{min}}-\frac{1}{R_{max}})+\frac{1}{R_{max}}}$ (21)

Where $R$ denotes the expected memristive RS, and $R_{min}$ and $R_{max}$ are
the lower and the upper RS boundaries that are the minimal/maximal of a set of
$r_{p,n}$ values with different bias voltage $v$. As shown in Equation 19, the
switching dynamics of memristor devices are non-linear and governed both by
bias voltages and the current memristive RSs. Therefore, it is challenging to
trigger memristors to reach certain RSs. Reference [52] proposed a ’predict-
write-verify’ loop to find the parameter sets (including bias voltages and the
pulse-width) that can lead to memristive RSs closest to expected. Given
multiple sets of bias voltages and the pulse-width, firstly, predict all
resulting memristor RSs, and find the one that contributes to the shortest
distance to the expected RSs. Next, apply the pulse with this chosen parameter
set, and verify the actual new RSs. Repeat the same procedure until the
differences between the actual RSs and the expected RSs are within an
acceptable range defined by parameter $\mathbf{R\
tolerance}=|R_{expected}-R_{actual}|/R_{expected}$. Considering a possible
scenario of infinite loops caused by the read noise, another parameter maxN is
also supplied to restrict the maximum iteration times.

When implementing this weight updating scheme, a critical design choice is to
choose pulsing parameter sets. Suppose the provided pulsing parameter can only
lead to small changes in RSs. In that case, it will take multiple steps to
reach the expected RSs and potentially lead to a significant error when
reaching the maximum iteration times set by the maxN. If parameter sets only
allow large-step RS changes, the targeted RSs may be ’missed’. Therefore, a
practical solution is to match the required RS changes and the estimated RS
changes that the selected parameter sets can result in. For approach 1, the
required RS changes usually are relatively large. Therefore, parameter sets
that can lead to large-step and small-step RS changes are both needed: the
former is used to make RSs converge quickly, and the latter is for fine-
tuning. Approach 2, however, only needs parameter sets that lead to relatively
smaller-step RS changes because the required RS changes in each training step
are relatively small. To make sure the weight updating can meet the stop
condition set by the R tolerance, for both approaches, parameter sets that can
lead to RS update smaller than the RS margin defined by the R tolerance need
to be chosen.

##  3 Results

###  3.1 Experiment Overview

Figure 3: (a) Building a vocabulary using the training samples. Step 1: Create
a look-up table for storing word frequency. Step 2: Omit the words that appear
less than 10x to reduce the vocabulary size. Step 3: Iterate through all
samples to transform them into word IDs using the vocabulary. Words are tagged
as ’[unk]’ if they cannot be found in the vocabulary. When packed into
batches, samples in the same batch are padded with the ’[pad]’ token to have
the same sentence length. (b) Training and validation accuracy evolution
curves in the ANN (approach 1), the SNN (approach 2), and the SNN with
memristor models (approach 2).

We demonstrate a sentiment analysis task with the IMDB movie reviews dataset
to validate the two approaches to obtaining a trained memristor-based SNN. The
dataset includes 25k highly polar movie reviews training samples and 25k test
samples. Inspired by the GitHub project [70], we firstly use the training
samples to build a vocabulary to pre-process the input data (Figure 3 (a)). We
count the word frequency in the training samples and put those words that
appear over 10x times into the vocabulary. This step creates a vocabulary of
20473 words. Next, we split the original training set into a sub-training set
with 17.5k samples and a validation set with 7.5k samples. After that, an
input review is tokenised using the basic English tokeniser and further
converted into a word ID vector. The one-hot representation of the word ID
vector is then fed to the neural network to start an inference. The word
embedding layer contains pre-trained GloVe word embeddings [47] with the
embedding dimension of 100. For the training and the inference processes,
please refer to Section ’Methods’. In this experiment, the network is trained
for 5 epochs. For a single training epoch, after training using the entire
training set with 17.5k samples, the training is switched off, and 7.5k
samples from the validation set are fed into the network to validate the
training effect for this epoch. The training and the validation set samples
are shuffled before starting another epoch. The trainable network parameters
that lead to the smallest validation loss are reloaded to the network for
testing after training. 25k samples from an independent test set are then sent
to the network. Experiment parameters for the two approaches are summarised in
Table 2, accuracy evolutions during the training process in Figure 3 (b), and
the final test accuracy in Table 1. The baseline test accuracy is 86.02%,
given by the ANN in approach 1. According to [70], to achieve higher
classification accuracy, more complex architectures are needed such as
transformers [71]. We can notice from Table 1 that, in approach 1, the test
accuracy degradations from the ANN to the SNN and the SNN to the SNN with
memristor models are very small (0.13% and 0.01% respectively). In approach 2,
the degradation from the ANN to the SNN is 0.1%, whilst the test accuracy
drops 1.06% by adding memristor models.

Table 1: Test accuracy for approach 1 and 2.

  * | ANN | SNN | SNN with memristors | SNN | SNN with memristors
---|---|---|---|---|---
approach | 1 | 1 | 1 | 2 | 2
accuracy (%) | 86.02 | 85.89 | 85.88 | 85.92 | 84.86

Table 2: Parameters used in the sentiment analysis tasks baseline
configuration for Approach 1 and 2. Different parameter values for different
approaches are highlighted in light grey.

  * Parameters | Approach 1 | Approach 2
---|---|---
Training set size | 17.5k | 17.5k
Validation set size | 7.5k | 7.5k
Test set size | 25k | 25k
Vocabulary size | 20473 | 20473
Embedding dimension | 100 | 100
Output dimension | 1 | 1
Batch size | 1 | 1
Offset | -25 | -0.5
T | 1k | 1k
$V_{th}$ | 50 | 56.75
$\eta$ | 0.05 | 0.05
$\epsilon$ | $1\times 10^{-8}$ | $1\times 10^{-8}$
Array size | 10 $\times$ 10 | 10 $\times$ 10
$A_{p}$ | 0.21389 | 0.21389
$A_{n}$ | -0.81302 | -0.81302
$a_{0p}$ | 37087 | 37087
$a_{0n}$ | 43430 | 43430
$a_{1p}$ | -20193 | -20193
$a_{1n}$ | 34333 | 34333
$t_{p}$ | 1.6591 | 1.6591
$t_{n}$ | 1.5148 | 1.5148
Positive pulse magnitude (V) | 0.9 | 0.9
Positive pulse duration (us) | 1, 2, 10, 20, 50, 100 | 1, 2, 10, 20, 50
Negative pulse magnitude (V) | -1.2 | -1.2
Negative pulse duration (us) | 1, 2, 10, 20, 100, 1000, 2000, 5000 | 1, 2, 10, 20, 100
dt (ms) | 1 | 1
R tolerance | 0.05% | 0.05%
MaxN | 5 | 5

###  3.2 Results Analysis of Approach 1

Figure 4: Result analysis for converting an ANN to a memristor-based SNN. (a)
Test accuracy vs spike train length, where a train represents a single
continuous value. (b) Test accuracy vs the R tolerance values when mapping
weights to memristor RSs. (c) Weight standard deviations vs R tolerance when
mapping weights to memristor RSs. (d) The measured weights along with the
actual weights with the read noise of 10%. (e) The test accuracy dependency on
the read noise. (f) Weight standard deviations vs read noise.

Now we look into the results from approach 1. Figure 4 (a) shows how the spike
train length for representing a single numerical value affects the
classification accuracy. From the curve, we can see that the accuracy shows an
’increased-and-saturated’ tendency along with the increase of the spike train
length because the spike trains are stochastic. The longer the spike trains,
the more accurately they can represent the numerical values. However,
increasing the length will not change the spike representation accuracy when
the spike train is long enough. Therefore, system performance plateaus. Longer
spike trains also require longer runtime to iterate through the whole spike
train for processing a single sample. Therefore, we choose T = 1k for
experiments as the standard configuration in this section, though longer spike
trains can result in more accurate representations.

System sensitivity to the R tolerance when using approach 1 is given in Figure
4 (b). In theory, the larger the R tolerance, the larger the maximum errors
between the expected memristor RSs and the actual values. Therefore, a large R
tolerance value could lead to performance degradation. However, as shown in
the bar chart, the system is robust to large R tolerance values. To
investigate the reasons, we plot the Euclidean distance between the expected
and the actual weights along with the increase of the read noise (Figure 4
(c)). Even when the R tolerance is as large as 20%, the Euclidean distance is
not higher than 0.06. This indicates that the system can tolerate variations
in exporting weights into memristor arrays when converting a trained ANN to an
SNN with memristors.

Similarly, approach 1 also exhibits robustness to the read noise of the read
operation for memristor arrays (Figure 4 (e)). In practice, the instruments
used to read memristive RSs introduce read noise, potentially affecting system
performance. The read noise in the framework is simulated by adding a random
number that does not exceed a certain percentage of the memristor RSs. In
other words, when we say the read noise is 10%, it means that the errors
caused by the read noise are no more than $\pm$10% of the memristor RSs
(Figure 4 (d)). The read noise is both included in the prediction stage of the
weight updating scheme and the inference. In this simulation, even with 20%
read noise, which means the error between the memristor RSs shown by the
instruments and the actual values is up to 20% in the worst case, the system
can still achieve a classification accuracy of 81.79%. Figure 4 (f) shows the
Euclidean distance between the actual and the measured weights with different
read noise values. The Euclidean distance caused by 20% read noise is around
0.105 (Figure 4 (f)).

###  3.3 Results Analysis of Approach 2

Figure 5: Result analysis for training the memristor-based SNN directly. (a)
Weight evolution curve for synapse 0 in the single-layer SNN during training
epoch 1. The curve is plotted every 175 samples. The calculated expected
weights are plotted in blue, and the actual weights mapped into memristor
arrays are in orange. Notably, memristive devices stop updating once the
requested updates drop below the R tolerance. (b)-(c) Training/validation and
test accuracy with different R tolerance values. (d)-(e) Measured weight
evolution for synapse 0 in the training epoch 1 with different read noise
values when the R tolerance is 0.5% and 3.0%. The curve in blue is the
baseline without the read noise. Weights are plotted every 175 samples. (f)
Shmoo plot for the test accuracy (%) with different read noise values (x axis)
and different R tolerance (y axis). Positive and negative results are coloured
blue and yellow, respectively.

Now we investigate how the R tolerance affects weight updating when using
approach 2. Figure 5 (a) plots the expected weight values evolution in the
synapse 0 in epoch 1 when training the memristor-based SNN directly (in the
blue curve) and the actual weights written in memristors (as orange dots). The
curves are plotted every 175 samples. The expected weight (in the blue curve)
progressively converges to the optima. However, the actual weights written in
memristors stop updating after a certain point because the weight updating is
cut off by the R tolerance when the weight change is smaller than the value
set by the R tolerance. We envisage that the larger the R tolerance value, the
earlier the weight updating stops. Figure 2.3 (b) and (c) give the
training/validation accuracy curves and the final test accuracy with different
R tolerance values. The accuracy shows no noticeable change when the R
tolerance is smaller than 4.0%. However, when the R tolerance is increased to
5%, the network cannot give reasonable classification predictions due to the
early cut-off.

Contrary to the read noise-robustness of the system built using the ANN
conversion, directly training the memristor-based SNN is sensitive to the read
noise. Figure 5 (d) depicts the measured weight evolution in synapse 0 during
the training epoch 1 with different read noise values when the R tolerance is
0.5%. Weights are plotted every 175 samples. The ideal situation gives the
baseline without the read noise in blue. The weight in the baseline converges
to an optimal state with the value around 0.328, and it converges to around
0.330 with 0.1% read noise. However, the weight evolutions during epoch 1 with
the read noise higher than the R tolerance (specifically, with 0.5% and 1.0%
read noise) fail to converge. We speculate this is because weights cannot be
correctly updated when the variations introduced by the read noise exceed the
acceptable range defined by the R tolerance. The effects can be multi-folded:
the relatively large measured errors confuse the weight updating scheme to
determine the weight change magnitudes or even directions, introduce weight
wobblings, and make the R tolerance unable to be used as an effective stop
condition. With 3.0% of R tolerance, which is higher than the largest read
noise we have tested, weights in all conditions converge (shown in Figure 5
(e)). When the R tolerance is higher than the read noise, the major errors
during the weight updating are introduced by the R tolerance rather than the
read noise so that the impact of the read noise can be ignored. Two features
can be noticed in Figure 5 (e)): firstly, more significant wobblings exist
along with the increase of the read noise, with the largest wobbling of
$\sim\pm 0.005$; secondly, the larger the read noise, the farther the measured
convergent state apart from the baseline, with the largest distance of 0.03.
However, both the fluctuations and the distances seem to be within acceptable
ranges. Figure 5 (f) gives the test accuracy with different read noise and R
tolerance. We conclude that when the R tolerance is higher than the read
noise, weights converge, and the test accuracy shows no noticeable difference
from the baseline; otherwise, weights fail to converge, and the final test
accuracy shows negative results.

##  4 Conclusion

In this paper, we developed a simulation framework incorporating an empirical
memristor model to demonstrate the first text classification task in
memristor-based spiking neural networks with the IMDB movie reviews dataset.
We took two paths to obtain trained memristor-based spiking neural networks by
1) converting a pre-trained ANN to a memristor-based SNN or 2) training a
memristor-based SNN directly using the equivalent ANN-based learning rules. We
elaborated on the specific details of estimating the critical parameters to
guide users to achieve comparable performance when taking the presented
methods. We achieved the test accuracy of 85.88% and 84.86% from the two
approaches with only 0.14% and 1.16% degradations, respectively, given the
baseline accuracy of 86.02% from the equivalent ANN. Finally, we investigated
how the SNN stochasticity and memristor non-idealities affect the system
performance. We concluded that it is possible to achieve very similar
performance in simulation from ANNs to memristor-based SNNs and from non-
memristive synapses to data-driven virtual memristive synapses. Proven by
demonstrating a sentiment analysis task, using the simulation framework
presented in this paper can accelerate research in employing memristor-based
SNNs in text classification tasks by overcoming the lack of learning rules and
memristor non-idealities.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Data Availability Statement

The original code of this work can be found in https://github.com/hjq310/text-
classification-in-memristorsnn.

## References

## References

  * [1] Abhronil Sengupta, Yuting Ye, Robert Wang, Chiao Liu, and Kaushik Roy.  Going deeper in spiking neural networks: Vgg and residual architectures.  Frontiers in Neuroscience, 13:95, 2019.
  * [2] Wolfgang Maass.  Networks of spiking neurons: The third generation of neural network models.  Neural Networks, 10(9):1659–1671, 1997.
  * [3] Lyle Long and Guoliang Fang.  A review of biologically plausible neuron models for spiking neural networks.  AIAA Infotech@ Aerospace 2010, page 3540, 2010.
  * [4] Hadar Cohen Duwek and Elishai Ezra Tsur.  Biologically plausible spiking neural networks for perceptual filling-in.  In Proceedings of the Annual Meeting of the Cognitive Science Society, volume 43, 2021.
  * [5] Ping Tak Peter Tang, Tsung-Han Lin, and Mike Davies.  Sparse coding by spiking neural networks: Convergence theory and computational results, 2017.  arXiv:1705.05475.
  * [6] K. Ahmed, A. Shrestha, Q. Qiu, and Q. Wu.  Probabilistic inference using stochastic spiking neural networks on a neurosynaptic processor.  In 2016 International Joint Conference on Neural Networks (IJCNN), pages 4286–4293, 2016.
  * [7] Long Cheng, Yang Liu, Zeng-Guang Hou, Min Tan, Dajun Du, and Minrui Fei.  A rapid spiking neural network approach with an application on hand gesture recognition.  IEEE Transactions on Cognitive and Developmental Systems, 13(1):151–161, 2021.
  * [8] Seijoon Kim, Seongsik Park, Byunggook Na, and Sungroh Yoon.  Spiking-yolo: Spiking neural network for energy-efficient object detection.  In AAAI, 2020.
  * [9] Wang Wei Lee, Sunil L. Kukreja, and Nitish V. Thakor.  Cone: Convex-optimized-synaptic efficacies for temporally precise spike mapping.  IEEE Transactions on Neural Networks and Learning Systems, 28(4):849–861, 2017.
  * [10] David Reid, Abir Jaafar Hussain, and Hissam Tawfik.  Financial time series prediction using spiking neural networks.  PloS one, 9(8):e103656–e103656, 2014.
  * [11] Chaitanya Prasad N, Krishnakant Saboo, and Bipin Rajendran.  Composer classification based on temporal coding in adaptive spiking neural networks.  In 2015 International Joint Conference on Neural Networks (IJCNN), pages 1–8, 2015.
  * [12] Martin J. Pearson, A. G. Pipe, B. Mitchinson, K. Gurney, C. Melhuish, I. Gilhespy, and M. Nibouche.  Implementing spiking neural networks for real-time signal-processing and control applications: A model-validated fpga approach.  IEEE Transactions on Neural Networks, 18(5):1472–1487, 2007.
  * [13] Samanwoy Ghosh-Dastidar and Hojjat Adeli.  A new supervised learning algorithm for multiple spiking neural networks with application in epilepsy and seizure detection.  Neural Networks, 22(10):1419–1431, 2009.
  * [14] Elmar Rueckert, David Kappel, Daniel Tanneberg, Dejan Pecevski, and Jan Peters.  Recurrent spiking networks solve planning tasks.  Scientific Reports, 6(1):21142, 2016.
  * [15] Xiuqing Wang, Zeng-Guang Hou, Feng Lv, Min Tan, and Yongji Wang.  Mobile robots’ modular navigation controller using spiking neural networks.  Neurocomputing, 134:230–238, 2014.  Special issue on the 2011 Sino-foreign-interchange Workshop on Intelligence Science and Intelligent Data Engineering (IScIDE 2011) Learning Algorithms and Applications.
  * [16] Zhenshan Bing, Claus Meschede, Florian Röhrbein, Kai Huang, and Alois C. Knoll.  A survey of robotics control based on learning-inspired spiking neural networks.  Frontiers in Neurorobotics, 12, 2018.
  * [17] Ulrich Roth, Axel Jahnke, and Heinrich Klar.  Hardware requirements for spike-processing neural networks.  In From Natural to Artificial Neural Computation, pages 720–727, Berlin, Heidelberg, 1995. Springer Berlin Heidelberg.
  * [18] John Backus.  Can programming be liberated from the von neumann style? a functional style and its algebra of programs.  Commun. ACM, 21(8):613–641, August 1978.
  * [19] S. K. Gonugondla, M. Kang, and N. Shanbhag.  A 42pj/decision 3.12tops/w robust in-memory machine learning classifier with on-chip training.  In 2018 IEEE International Solid - State Circuits Conference - (ISSCC), pages 490–492, 2018.
  * [20] Norman P. Jouppi, Cliff Young, Nishant Patil, David Patterson, Gaurav Agrawal, Raminder Bajwa, Sarah Bates, Suresh Bhatia, et al.  In-datacenter performance analysis of a tensor processing unit, 2017.
  * [21] Mike Davies, Narayan Srinivasa, Tsung-Han Lin, Gautham Chinya, Yongqiang Cao, Sri Harsha Choday, Georgios Dimou, Prasad Joshi, Nabil Imam, Shweta Jain, Yuyun Liao, et al.  Loihi: A neuromorphic manycore processor with on-chip learning.  IEEE Micro, 38(1):82–99, 2018.
  * [22] Filipp Akopyan, Jun Sawada, Andrew Cassidy, Rodrigo Alvarez-Icaza, John Arthur, Paul Merolla, Nabil Imam, Yutaka Nakamura, Pallab Datta, Gi-Joon Nam, et al.  Truenorth: Design and tool flow of a 65 mw 1 million neuron programmable neurosynaptic chip.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 34(10):1537–1557, 2015.
  * [23] Ben Varkey Benjamin, Peiran Gao, Emmett McQuinn, Swadesh Choudhary, Anand R. Chandrasekaran, Jean-Marie Bussat, Rodrigo Alvarez-Icaza, John V. Arthur, Paul A. Merolla, and Kwabena Boahen.  Neurogrid: A mixed-analog-digital multichip system for large-scale neural simulations.  Proceedings of the IEEE, 102(5):699–716, 2014.
  * [24] Eustace Painkras, Luis A. Plana, Jim Garside, Steve Temple, Francesco Galluppi, Cameron Patterson, David R. Lester, Andrew D. Brown, and Steve B. Furber.  Spinnaker: A 1-w 18-core system-on-chip for massively-parallel neural network simulation.  IEEE Journal of Solid-State Circuits, 48(8):1943–1953, 2013.
  * [25] Stefan Scholze, Holger Eisenreich, Sebastian Höppner, Georg Ellguth, Stephan Henker, Mario Ander, Stefan Hänzsche, Johannes Partzsch, Christian Mayr, and René Schüffny.  A 32 gbit/s communication soc for a waferscale neuromorphic system.  Integration the VLSI Journal, 45:61–75, 01 2012.
  * [26] Ali Khiat, Peter Ayliffe, and Themis Prodromakis.  High density crossbar arrays with sub- 15 nm single cells via liftoff process only.  Scientific Reports, 6(32614):1–8, September 2016.
  * [27] Can Li, Yunning Li, Hao Jiang, Wenhao Song, Peng Lin, Zhongrui Wang, J. Joshua Yang, Qiangfei Xia, Miao Hu, Eric Montgomery, et al.  Large memristor crossbars for analog computing.  In 2018 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–4, 2018.
  * [28] Spyros Stathopoulos, Ali Khiat, Maria Trapatseli, Simone Cortese, Alexander Serb, Ilia Valov, and Themis Prodromakis.  Multibit memory operation of metal-oxide bi-layer memristors.  Scientific Reports, 7, 12 2017.
  * [29] Maheswari Sivan, Yida Li, Hasita Veluri, Yunshan Zhao, Baoshan Tang, Xinghua Wang, Evgeny Zamburg, Jin Leong, Jessie Niu, Umesh Chand, et al.  All wse2 1t1r resistive ram cell for future monolithic 3d embedded memory integration.  Nature Communications, 10, 11 2019.
  * [30] Byung Joon Choi, Antonio Torrezan, John William Strachan, P. Kotula, Andrew Lohn, Matthew Marinella, Zhiyong Li, Stan Williams, and Jianhua Joshua Yang.  High-speed and low-energy nitride memristors.  Advanced Functional Materials, 26, 05 2016.
  * [31] L. Goux, A. Fantini, G. Kar, Y.-Y. Chen, N. Jossart, R. Degraeve, S. Clima, B. Govoreanu, G. Lorenzo, G. Pourtois, et al.  Ultralow sub-500na operating current high-performance tin $\backslash$ al2o3 $\backslash$ hfo2 $\backslash$ hf $\backslash$ tin bipolar rram achieved through understanding-based stack-engineering.  In 2012 Symposium on VLSI Technology (VLSIT), pages 159–160, 2012\.
  * [32] Alexantrou Serb, Edoardo Manino, Ioannis Messaris, Long Tran-Thanh, and Themis Prodromakis.  Hardware-level bayesian inference.  In Neural Information Processing Systems, December 2017.
  * [33] Peng Yao, Huaqiang Wu, Bin Gao, Jianshi Tang, Qingtian Zhang, Wenqiang Zhang, J Joshua Yang, and He Qian.  Fully hardware-implemented memristor convolutional neural network.  Nature, 577(7792):641–646, 2020.
  * [34] Stefano Ambrogio, Pritish Narayanan, Hsinyu Tsai, Robert M Shelby, Irem Boybat, Carmelo di Nolfo, Severin Sidler, et al.  Equivalent-accuracy accelerated neural-network training using analogue memory.  Nature, 558(7708):60–67, jun 2018.
  * [35] Alexander Serb, Johannes Bill, Ali Khiat, Radu Berdan, Robert Legenstein, and Themis Prodromakis.  Unsupervised learning in probabilistic neural networks with multi-state metal-oxide memristive synapses.  Nature Communications, 7, 2016.
  * [36] Erika Covi, Stefano Brivio, Alexander Serb, Themis Prodromakis, Marco Fanciulli, and Sabina Spiga.  Analog memristive synapse in spiking networks implementing unsupervised learning.  Frontiers in Neuroscience, 10:482, 2016.
  * [37] Henry Markram, Joachim Lübke, Michael Frotscher, and Bert Sakmann.  Regulation of synaptic efficacy by coincidence of postsynaptic aps and epsps.  Science (New York, N.Y.), 275:213–5, 02 1997.
  * [38] Melika Payvand, Mohammed E Fouda, Fadi Kurdahi, Ahmed M Eltawil, Senior Member, and Emre O Neftci.  On-Chip Error-Triggered Learning of Multi-Layer Memristive Spiking Neural Networks.  10(4):522–535, 2020.
  * [39] Irem Boybat, Manuel Le Gallo, S R Nandakumar, Timoleon Moraitis, Thomas Parnell, Tomas Tuma, Bipin Rajendran, Yusuf Leblebici, Abu Sebastian, and Evangelos Eleftheriou.  Neuromorphic computing with multi-memristive synapses.  Nature Communications, 9(1):1–15, 2018.
  * [40] Beiye Liu, Yiran Chen, Bryant Wysocki, and Tingwen Huang.  Reconfigurable neuromorphic computing system with memristor-based synapse design.  Neural Processing Letters, 41, 04 2013.
  * [41] Chenchen Liu, Bonan Yan, Chaofei Yang, Linghao Song, Zheng Li, Beiye Liu, Yiran Chen, Hai Li, Qing Wu, and Hao Jiang.  A spiking neuromorphic design with resistive crossbar.  In 2015 52nd ACM/EDAC/IEEE Design Automation Conference (DAC), pages 1–6, 2015.
  * [42] Stanisław Woźniak, Angeliki Pantazi, Thomas Bohnstingl, and Evangelos Eleftheriou.  Deep learning incorporating biologically inspired neural dynamics and in-memory computing.  Nature Machine Intelligence, 2(6):325–336, 2020.
  * [43] Yigit Demirag, Charlotte Frenkel, Melika Payvand, and Giacomo Indiveri.  Online Training of Spiking Recurrent Neural Networks with Phase-Change Memory Synapses.  2021\.  arXiv:1410.5093.
  * [44] Guillaume Bellec, Franz Scherr, Anand Subramoney, Elias Hajek, Darjan Salaj, Robert Legenstein, and Wolfgang Maass.  A solution to the learning dilemma for recurrent networks of spiking neurons.  Nature Communications, 11(1):3625, 2020.
  * [45] Emre O. Neftci, Hesham Mostafa, and Friedemann Zenke.  Surrogate gradient learning in spiking neural networks, 2019.  arXiv: 1901.09948.
  * [46] Peter O’Connor, Daniel Neil, Shih-Chii Liu, Tobi Delbruck, and Michael Pfeiffer.  Real-time classification and sensor fusion with a spiking deep belief network.  Frontiers in Neuroscience, 7, 2013.
  * [47] Jeffrey Pennington, Richard Socher, and Christopher D. Manning.  Glove: Global vectors for word representation.  In Empirical Methods in Natural Language Processing (EMNLP), pages 1532–1543, 2014.
  * [48] Tomas Mikolov, Kai Chen, Greg Corrado, and Jeffrey Dean.  Efficient estimation of word representations in vector space, 2013.  arXiv:1301.3781.
  * [49] D. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams.  Learning representations by back-propagating errors.  Nature, 323:533–536, 1986.
  * [50] Yuwei Wang, Yi Zeng, Jianbo Tang, and Bo Xu.  Biological neuron coding inspired binary word embeddings.  Cognitive Computation, 11(5):676–684, 2019.
  * [51] Marcin Białas, Marcin Michał Mirończuk, and Jacek Mańdziuk.  Biologically plausible learning of text representation with spiking neural networks.  In Parallel Problem Solving from Nature – PPSN XVI, pages 433–447, Cham, 2020. Springer International Publishing.
  * [52] Jinqi Huang, Spyros Stathopoulos, Alexantrou Serb, and Themis Prodromakis.  Neuropack: An algorithm-level python-based simulator for memristor-empowered neuro-inspired computing.  Frontiers in Nanotechnology, 4, 2022.
  * [53] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.  Pytorch: An imperative style, high-performance deep learning library.  In Advances in Neural Information Processing Systems 32, pages 8024–8035. Curran Associates, Inc., 2019.
  * [54] Yiannis Messaris, Alexander Serb, Ali Khiat, Spyridon Nikolaidis, and Themis Prodromakis.  A compact verilog-a reram switching model.  03 2017.
  * [55] Andrew L. Maas, Raymond E. Daly, Peter T. Pham, Dan Huang, Andrew Y. Ng, and Christopher Potts.  Learning word vectors for sentiment analysis.  In Proceedings of the 49th Annual Meeting of the Association for Computational Linguistics: Human Language Technologies, pages 142–150, Portland, Oregon, USA, June 2011. Association for Computational Linguistics.
  * [56] Yongqiang Cao, Yang Chen, and Deepak Khosla.  Spiking deep convolutional neural networks for energy-efficient object recognition.  International Journal of Computer Vision, 113(1):54–66, 2015.
  * [57] Abien Fred Agarap.  Deep learning using rectified linear units (relu).  2018\.  arXiv:1803.08375.
  * [58] Bodo Rueckauer, Iulia-Alexandra Lungu, Yuhuang Hu, Michael Pfeiffer, and Shih-Chii Liu.  Conversion of continuous-valued deep networks to efficient event-driven networks for image classification.  Frontiers in Neuroscience, 11, 2017.
  * [59] Peter U. Diehl, Daniel Neil, Jonathan Binas, Matthew Cook, Shih-Chii Liu, and Michael Pfeiffer.  Fast-classifying, high-accuracy spiking deep networks through weight and threshold balancing.  In 2015 International Joint Conference on Neural Networks (IJCNN), pages 1–8, 2015.
  * [60] Davide Zambrano and Sander M Bohte.  Fast and efficient asynchronous neural computation with adapting spiking neural networks.  2016\.  arXiv:1609.02053.
  * [61] L.F Abbott.  Lapicque’s introduction of the integrate-and-fire model neuron (1907).  Brain Research Bulletin, 50(5):303–304, 1999.
  * [62] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  2016\.  arXiv:1609.04747.
  * [63] Geoffrey Hinton.  Coursera - neural networks for machine learning - geoffrey hinton.  2012\.
  * [64] Yoshua Bengio, Nicholas Léonard, and Aaron C. Courville.  Estimating or propagating gradients through stochastic neurons for conditional computation.  2013\.  arXiv:1308.3432.
  * [65] Shihui Yin, Shreyas K. Venkataramanaiah, Gregory K. Chen, Ram Krishnamurthy, Yu Cao, Chaitali Chakrabarti, and Jae-sun Seo.  Algorithm and hardware design of discrete-time spiking neural networks based on back propagation with binary activations.  In 2017 IEEE Biomedical Circuits and Systems Conference (BioCAS), pages 1–5, 2017.
  * [66] Agnes Lydia and Sagayaraj Francis.  Adagrad - an optimizer for stochastic gradient descent.  Volume 6:566–568, 05 2019.
  * [67] Yiannis Messaris, Spyridon Nikolaidis, Alexander Serb, Spyros Stathopoulos, Isha Gupta, Ali Khiat, and Themis Prodromakis.  A tio2 reram parameter extraction method.  pages 1–4, 05 2017.
  * [68] Radu Berdan, Alexander Serb, Ali Khiat, Anna Regoutz, Christos Papavassiliou, and Themis Prodromakis.  A $\mu$-controller-based system for interfacing selectorless rram crossbar arrays.  IEEE Transactions on Electron Devices, 62(7):2190–2196, 2015.
  * [69] Patrick Foster, Jinqi Huang, Alex Serb, Spyros Stathopoulos, Christos Papavassiliou, and Themis Prodromakis.  An fpga-based system for generalised electron devices testing.  2022\.  arXiv:2202.00499.
  * [70] Ben Trevett.  pytorch-sentiment-analysis, 2021.  https://github.com/bentrevett/pytorch-sentiment-analysis.
  * [71] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.


# Electric Penrose process: high-energy acceleration of ionized particles by
non-rotating weakly charged black hole

Arman Tursunov<EMAIL_ADDRESS>Bakhtinur Juraev  Zdeněk
Stuchlík  Martin Kološ  Research Centre for Theoretical Physics and
Astrophysics, Institute of Physics, Silesian University in Opava, Bezručovo
nám.13, CZ-74601 Opava, Czech Republic

###### Abstract

In many astrophysical scenarios the charge of the black hole is often
neglected due to unrealistically large values of the charge required for the
Reissner-Nordström spacetime metric. However, black holes may possess small
electric charge due to various selective accretion mechanisms. In this paper
we investigate the effect of a small hypothetical electric charge of a
Schwarzschild black hole on the ionization of a freely falling neutral
particle and subsequent escape of the ionized particle from the black hole. We
show that the energy of ionized particle can grow ultra-high and discuss
distinguishing signatures of particle acceleration by weakly charged black
holes. We also discuss a possible application of the proposed mechanism as an
alternative cosmic ray acceleration scenario. In particular we show that the
Galactic centre supermassive black hole is capable to act as a PeVatron of
protons. The presented mechanism can serve as a simple toy model of a non-
rotating compact object acting as a particle accelerator with a potential
astrophysical implementations related to the cosmic ray physics and beyond.

††preprint: APS/123-QED

##  I Weakly charged black hole

Recently, it has been pointed out that the ionization or decay of neutral
particles in the vicinity of rotating Kerr black hole immersed into external
magnetic field can lead to the acceleration of ionized particles to ultra-high
energies, with the Lorentz $\gamma$-factors of particles that may exceeding
$10^{12}$ near supermassive black holes in realistically plausible conditions
[1, 2]. The formalism of the mentioned acceleration mechanism is based on the
magnetic Penrose process [3, 4] in its novel, ultra-efficient regime, in which
the energy of ionized particle drives away the rotational energy of the black
hole through electromagnetic interaction. It has been claimed that the
mechanism might be responsible for the production of the highest-energy cosmic
ray particles [1] with energy exceeding $10^{20}$eV when applied to realistic
supermassive black hole candidates. The process, however, requires the
rotation of the black hole and the presence of external magnetic field.

In this paper we investigate whether the acceleration of ionized particles can
be achieved in a more simplified settings, namely, in the vicinity of a non-
rotating Schwarzschild black hole with a radial test electric field. By a test
electric field we denote the field, whose energy-momentum tensor can be
neglected in the description of the gravitational field of the black hole.
This implies that the electric field influences on the dynamics of charged
particles only, being negligible for the geodesics of neutral particles. Such
a simplified setup is motivated by the following reasons.

First of all, the no-hair theorem of black hole physics states that the
spacetime around black holes can be fully described by at most three metric
parameters – black hole mass, spin and electric charge [5]. The later is
usually neglected in astrophysical scenarios justified, on the one hand, by
unrealistically large values of the charge required for its visible effect on
the spacetime metric and on the other hand by quick discharge of any charge
excess by an accretion of a plasma surrounding black hole. Indeed, one can
compare the gravitational radius of a black hole with the characteristic
length of the charge $Q_{\rm G}$ of the Reissner-Nordström black hole, which
gives the maximum charge of the black hole per solar mass111Hereafter in this
section we use the esu-cgs system of units, in which the electrostatic unit of
charge is given by 1 Fr $\equiv$ 1 esu = 1 cm3/2 g1/2 s-1. In SI system of
units, 1 C $=3\times 10^{9}$ Fr.

$\displaystyle\sqrt{\frac{Q_{\rm G}^{2}G}{c^{4}}}=\frac{2GM}{c^{2}},$ (1)
$\displaystyle\Rightarrow$ $\displaystyle Q_{\rm G}=2G^{1/2}M\approx
10^{30}\frac{M}{M_{\odot}}~{}\mbox{Fr}.$ (2)

This value of the charge is unattainable in any known astrophysically relevant
scenario. Moreover, a neutralization of such a hypothetical charge $Q_{\rm G}$
would require an accretion of a net charge with the total mass of accreted
charged particles of $M_{Q}=m_{p,e}Q_{\rm G}/e$, where the indices denote
proton and electron, respectively, and $e$ is an elementary charge. Luminosity
of black hole surrounded by plasma or accretion disk can be derived from
infalling matter as $L=\epsilon\dot{M}c^{2}$, where $\dot{M}$ accretion rate
and $\epsilon$ is the fraction of the rest mass energy radiated away. On the
other hand, from the balance of gravitational force and radiation pressure in
the vicinity of a black hole one can derive the Eddington luminosity for fully
ionized hydrogen plasma surrounding a black hole in the form

$L_{\rm Edd}=\frac{4\pi GMm_{p}c}{\sigma_{\rm T}}\approx 1.26\times
10^{38}\left(\frac{M}{M_{\odot}}\right){\rm erg/s}.$ (3)

Defining charged matter accretion rate as the fraction of total accretion
rate, $\dot{M}_{Q}=\delta\cdot\dot{M}$, one can derive the neutralization
timescale of the maximally charged black hole in the following form

$t_{Q,{\rm
acc.}}=\frac{4}{3}\frac{e^{3}\,\epsilon\,m_{p,e}}{G^{1/2}\,c^{3}\,\delta\,m_{p}m_{e}^{2}}\approx
2.5\times
10^{-2}\left(\frac{m_{p,e}}{m_{p}}\right)\left(\frac{\epsilon}{\delta}\right){\rm
s},$ (4)

which is estimated for a positive black hole charge. In case of a negative
charge of the black hole, the timescale is $\approx 1835$ times faster. Both
${\epsilon}$ and $\delta$ have values in the range $(0\,,\,1)$ and in many
cases are of similar orders of magnitude. This implies that in all
astrophysically relevant settings any net charge of a black hole would be
neutralized relatively quickly unless there is a mechanism preventing the
black hole from neutralization. Some of such mechanisms are briefly mentioned
below. Thus, the Reissner-Nordström spacetime metric is interesting, but
astrophysically not viable. An exception can be represented by the Reissner-
Nordström spacetime having a zero electric and non-zero magnetic charge due to
special character of interaction with electrically charged matter [6].

There exist several astrophysical scenarios based on a selective accretion, in
which an astrophysical black hole may possess a small electric charge [7, 8,
9, 10, 11, 12, 13, 14]. Since protons are about 1836 times more massive than
electrons, the balance between the gravitational and Coulombic forces for the
particles close to the surface of the compact object is obtained when the
black hole acquires a positive net electric charge of the order of $Q\sim
3\times 10^{11}$Fr per solar mass [11, 8]. Moreover, matter surrounding black
hole can be ionized and charged by the irradiating photons taking away some
electrons [9]. In that case, following the above described argument, the
charge of the black hole is likely positive, being of the order of $Q\sim
10^{11}$Fr per solar mass. Another famous mechanism of charging of black holes
is the Wald’s mechanism [10], in which the charge is naturally induced by the
twisting of magnetic field lines due to the frame-dragging effect of the
rotation of black hole. As a result, both the black hole and surrounding
magnetosphere should acquire equal and opposite charge of the order of $Q\sim
10^{18}$Fr per solar mass [see, e.g. 7, 10, 15]. In all cases, the charge of
the black hole is much weaker than its maximal theoretical limit (2) by many
orders of magnitude. Therefore, astrophysical black holes can be considered as
weakly charged, i.e. gravitational effect of the charge on the spacetime
metric can be rightly neglected. Thus, depending on whether the black hole is
spinning or not, the realistic value of the black hole’s charge may vary
between the values

$10^{11}\frac{M}{M_{\odot}}\,{\rm Fr}\lesssim Q_{\rm BH}\lesssim
10^{18}\frac{M}{M_{\odot}}\,{\rm Fr}.$ (5)

For more details and estimates, the reader may refer to the works [12, 11] and
references therein, where various black hole charging scenarios are compared
and the results are applied to the Galactic centre black hole.

The presence of induced electric field and corresponding electric charge of a
black hole plays an important role in the mechanisms of the rotational energy
extraction from black holes, such as Blandford–Znajek mechanism [16] and
magnetic Penrose process [2]. Discharge of the induced charge by accreting
charged matter drives away the black hole’s rotational energy. The character
of the electric field around magnetized Kerr black hole is not spherically
symmetric, it is rather of quadrupole character. Therefore, in most of the
realistic particle acceleration scenarios, such as the production of
relativistic jets, it is assumed that the black hole is spinning.

In this paper, we raise a question whether the charged particles can be
accelerated to large $\gamma$-factors $>>1$ by non-spinning Schwarzschild
black holes or in the case when the rotation of the black hole is very slow
for the operation of the rotational energy extraction mechanisms? For that we
assume that the black hole is non-rotating with a small electric charge of
above described limits, whose gravitational contribution is negligible. We are
aimed to calculate the energy of the charged particle after the ionization of
a neutral particle in the black hole’s vicinity.

Hereafter we use the signature $(-,+,+,+)$, and the system of geometric units,
in which $G=1=c$, unless the constants are written explicitly and given in
esu-cgs system of units.

##  II Dynamics of charged particle

###  II.1 Background setup & equations of motion

We start from the Schwarzschild spacetime metric

$ds^{2}=-f(r)dt^{2}+f^{-1}(r)dr^{2}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\phi^{2}),$
(6)

where $f(r)$ is the lapse function parametrized by the black hole mass $M$ as
follows

$f(r)=1-\frac{2M}{r}.$ (7)

Let us assume the presence of the radial electric field with corresponding
small electric charge $Q$ at the center of the coordinate. In this case, the
only non-zero covariant component of the electromagnetic four-potential
$A_{\mu}=(A_{t},0,0,0)$ has the following simple form

$A_{t}=-\frac{Q}{r}.$ (8)

The anti-symmetric tensor of the electromagnetic field
$F_{\alpha\beta}=A_{\beta,\alpha}-A_{\alpha,\beta}$ has the only one
independent nonzero component

$F_{tr}=-F_{rt}=-\frac{Q}{r^{2}}.$ (9)

Figure 1: Radial dependence of the effective potential $V_{eff}$ for a charged
particle around weakly charged non-rotating black hole in the equatorial plane
$\theta=\pi/2$ for different values of the parameters ${\cal L}$ and ${\cal
Q}$.

Let us now consider the motion of a charged particle of mass $m$ and charge
$q$ in the combined background gravitational and electric fields. The motion
of a charged particle is governed by the Lorentz equation in curved spacetime

$\frac{\mathrm{d}u^{\mu}}{\mathrm{d}\tau}+\Gamma^{\mu}_{\alpha\beta}u^{\alpha}u^{\beta}-\frac{q}{m}F^{\mu}_{\,\,\,\nu}u^{\nu}=0,$
(10)

where $u^{\mu}$ is the four-velocity of the particle, $\tau$ is the proper
time of the particle and $\Gamma^{\mu}_{\alpha\beta}$ – Christoffel symbols.

Due to symmetries of the background Schwarzschild metric one can introduce two
integrals of motion, corresponding to temporal and spatial components of the
canonical four-momentum of the charged particle
$P_{\alpha}=mu_{\alpha}+qA_{\alpha}$:

$\displaystyle\frac{P_{t}}{m}$ $\displaystyle=$ $\displaystyle-{\cal
E}\equiv-\frac{E}{m}=u_{t}-\frac{qQ}{mr},$ (11)
$\displaystyle\frac{P_{\phi}}{m}$ $\displaystyle=$ $\displaystyle{\cal
L}\equiv\frac{L}{m}=u_{\phi},$ (12)

where ${\cal E}$ and ${\cal L}$ denote specific energy and specific angular
momentum of the charged particle. Since both gravitational and electric fields
are spherically symmetric and there is no preferred plane of the motion, one
can fix the motion of the charged particle to the equatorial plane
($\theta=\pi/2$), without loss of generality. Thus, three non-vanishing
components of the equation of motion (10) can be found in the form

$\displaystyle\frac{\mathrm{d}u^{t}}{\mathrm{d}\tau}$ $\displaystyle=$
$\displaystyle\frac{u^{r}\left[Qr-2M(er+Q)\right]}{r(r-2M)^{2}}$ (13)
$\displaystyle\frac{\mathrm{d}u^{r}}{\mathrm{d}\tau}$ $\displaystyle=$
$\displaystyle\frac{eQ}{r^{2}}+\frac{\mathcal{L}^{2}(r-2M)}{r^{4}}-\frac{M\left[e^{2}-(u^{r})^{2}\right]}{r(r-2M)},$
(14) $\displaystyle\frac{\mathrm{d}u^{\phi}}{\mathrm{d}\tau}$ $\displaystyle=$
$\displaystyle-\frac{2~{}\mathcal{L}~{}u^{r}}{r^{3}},$ (15) $\displaystyle{\rm
where}$ $\displaystyle e={\cal E}-\frac{qQ}{mr}.$ (16)

Equations (13) - (15) are ordinary differential equations, which can be easily
solved numerically.

###  II.2 Effective potential

Using the normalization condition for a massive particle $u^{\mu}u_{\mu}=-1$,
one can derive the effective potential for the charged particle moving around
a weakly charged Schwarzschild black hole in the form

$V_{\rm
eff}(r)=\frac{\mathcal{Q}}{r}+\sqrt{f(r)\left(1+\frac{\mathcal{L}^{2}}{r^{2}}\right)},$
(17)

where $\mathcal{Q}=Qq/m$ is a parameter characterizing the electric
interaction between the charges of the particle and black hole. Without loss
of generality we set the mass of the black hole to be equal to unity, i.e.
$M=1$.

Since the right hand side of the effective potential (17) is always positive
one can distiguish two qualitatively different situations depending on the
sign of the parameter $\mathcal{Q}$. When $\mathcal{Q}>0$, the charges of the
particle and black hole have the same sign, so the electric interaction is
repulsive. In opposite case, when $\mathcal{Q}<0$, the charges of the particle
and black hole have different signs, so the electric interaction is
attractive. The term ${\cal L}^{2}$ under the root of Eq.(17) means that the
clockwise and counter-clockwise directions of the motion are equivalent.

The radial profile of the effective potential is shown in Figure 1. One can
see that the effect of the charge parameter $\cal Q$ is similar to those of
the angular momentum $\cal L$, i.e. increasing (or decreasing) both parameters
$\cal Q$ and $\cal L$ one can increase (or decrease) the value of the
effective potential. It is interesting to note that taking into account the
parameter $\cal Q$ can mimic the effect of angular momentum (compare, e.g. red
curve in the middle plot with a very similar blue curve on the right plot of
the Figure 1).

Figure 2: Left: position of the ISCO of charged particle in the dependence on
the charge parameter $\cal Q$. Middle: Angular momentum of charged particle at
ISCO against ISCO position. Right: velocity of charged particle at ISCO. In
all plots the red lines correspond to the positive charge parameter ${\cal
Q}>0$, while black curves correspond to the negative charge parameter ${\cal
Q}<0$.

The stationary points of the effective potential $V_{\rm eff}(r)$ is given by
the equation

$\partial_{r}V_{\rm eff}(r)=0.$ (18)

Note that in the case of the weakly charged Schwarzschild black hole all the
local extrema of the effective potential $V_{\rm eff}$ are located in the
equatorial plane $\theta=\pi/2$. Equation (18) leads to a polynomial equation
of the fourth order in the radial coordinate

$r^{2}(J-1)+\mathcal{L}^{2}(r-3)=0,$ (19)

$\mbox{where}\quad
J=\frac{\mathcal{Q}}{r}\sqrt{\frac{(r-2)(\mathcal{L}^{2}+r^{2})}{r}}.$ (20)

The solution of equation (19) has four roots of $\mathcal{L}$ and two of them
being independent

$\begin{split}\mathcal{L}_{\pm}^{2}=\frac{r}{(r-3)^{2}}\Big{[}-\mathcal{Q}^{2}-3r+\frac{\mathcal{Q}^{2}r}{2}+r^{2}\\\
\pm\mathcal{Q}\sqrt{\mathcal{Q}^{2}-12r+4r^{2}}\left(1-\frac{r}{2}\right)\Big{]},\end{split}$
(21)

###  II.3 Angular velocity measured at infinity

Noticing that in the equatorial plane the four velocity takes the form
$u^{\alpha}=u^{t}(1,v,0,\Omega)$, where $v=dr/dt$, $\Omega=d\phi/dt$ and using
the normalization condition $u^{\alpha}u_{\alpha}=-k$, where $k=1$ for massive
particle and $k=0$ for massless particle, we can obtain the following equation

$(u^{t})^{2}(f^{-1}(r)v^{2}-f(r)+\Omega^{2}r^{2})=-k.$ (22)

Simplifying equation above, we can easily derive equation for angular velocity
measured by a static observer at infinity $\Omega=d\phi/dt$

$\Omega=\pm\frac{1}{u_{t}r}\sqrt{(u_{t})^{2}(f(r)-f^{-1}(r)v^{2})-kf^{2}(r)}.$
(23)

The possible values of $\Omega$ are limited to

$\Omega_{-}\leq\Omega\leq\Omega_{+},\quad\Omega_{\pm}=\pm\frac{\sqrt{f(r)}}{r}.$
(24)

corresponding to the photon motion.

###  II.4 Innermost stable circular orbit

The innermost stable circular orbit (ISCO) in the Schwarzschild spacetime is
located at $r_{\rm ISCO}=6M$. In the case when the electric charge is
included, it will be shifted from $r_{\rm ISCO}=6M$. Local extremum of the
function $\mathcal{L}_{\pm}$ determine the innermost stable circular orbit
(ISCO), its radius, angular momentum and energy. The ISCO can also be found
from the condition of $\partial^{2}_{r}V_{\rm
eff}(r,\mathcal{L},\mathcal{Q})=0$, which gives

$\begin{split}\mathcal{L}^{2}r^{2}(J(r-2)+2)+r^{4}(J(r-2)-r+3)\\\
+\mathcal{L}^{4}((r-3)r+3)=0.\end{split}$ (25)

Solving this equation with respect to $r$ gives us four solutions for the ISCO
with only two of them being real and independent.

One can also calculate the velocity $v$ of the charged particle at the ISCO,
which is given by the formula

$v=\sqrt{\frac{1}{1+{r^{2}_{\rm ISCO}}/{\mathcal{L}^{2}_{\rm ISCO}}}}.$ (26)

Dependence of the ISCO position $r_{\rm isco}$ on the charge parameter ${\cal
Q}$ and the change of the values of ${\cal L}_{\rm ISCO}$ and $v$ on the ISCO
position are shown in Figure 2. ISCO is increasing for both positive and
negative ${\cal Q}$. Similar results have been also obtained recently by [17],
where the ISCO in a similar setting is properly discussed.

Figure 3: Ratio of energies of ionized and neutral particles plotted against
the black hole charge $Q$ (left) and the ionization rate $Z/A$ (right). The
curves denote different positions of the ionization point: $r_{\rm
ion}=2GM/c^{2}$ (solid), $r_{\rm ion}=4GM/c^{2}$ (dashed) and $r_{\rm
ion}=6GM/c^{2}$ (dot-dashed).

##  III Energy of ionized particle

###  III.1 Conservation laws

Let us now consider the decay of a particle 1 into two fragments 2 and 3 close
to the event horizon of a weakly charged Schwarzschild black hole at the
equatorial plane. One can write the following conservation laws before and
after decay

$\displaystyle E_{1}=E_{2}+E_{3},\quad L_{1}=L_{2}+L_{3},\quad
q_{1}=q_{2}+q_{3},$ (27) $\displaystyle
m_{1}\dot{r_{1}}=m_{2}\dot{r_{2}}+m_{3}\dot{r_{3}},\quad m_{1}\geqslant
m_{2}+m_{3},$ (28)

where dot indicates derivatives with respect to the particle’s proper time
$\tau$. Using above conservation laws, one can find the equation

$m_{1}u_{1}^{\phi}=m_{2}u^{\phi}_{2}+m_{3}u^{\phi}_{3}.$ (29)

Noticing that $u^{\phi}=\Omega u^{t}=\Omega e/f(r)$, where
$e_{i}=(E_{i}+q_{i}A_{t})/m_{i}$, with $i=1,2,3$ indicating the particle’s
number, the equation (29) will take the following from

$\Omega_{1}m_{1}e_{1}=\Omega_{2}m_{2}e_{2}+\Omega_{3}m_{3}e_{3}.$ (30)

Solving the above equation with respect to the energy of one of the fragments,
e.g. $E_{3}$ we find

$E_{3}=\frac{\Omega_{1}-\Omega_{2}}{\Omega_{3}-\Omega_{2}}\left(E_{1}+q_{1}A_{t}\right)-q_{3}A_{t},$
(31)

where $\Omega_{i}=d\phi_{i}/dt$ is an angular velocity of $i$th particle,
given by (23), with restricted values (24).

###  III.2 Maximum energy of ionized particle

To maximize the energy of ionized particle we choose the particle 1 to be
neutral, i.e. $q_{1}=0$. We are also free to choose the energy of the particle
1, which we set to its rest mass energy, i.e. $E_{1}=m_{1}$ or ${\cal E}=1$.
In this case, the angular velocity (23) for the particle 1 will take the
following simple form

$\Omega_{1}=\frac{1}{r^{2}}\sqrt{2(r-2)}.$ (32)

Without loss of generality we choose the ionized particle to be the particle
3. The energy of the ionized particle is maximal, when the term
$(\Omega_{1}-\Omega_{2})/(\Omega_{3}-\Omega_{2})$ is maximized. This occurs
when we set the angular momentum of fragments to their limiting values. Then
we find

$\frac{\Omega_{1}-\Omega_{2}}{\Omega_{3}-\Omega_{2}}\Big{|}_{\rm
max}=\frac{1}{\sqrt{2~{}r_{\rm ion}}}+\frac{1}{2},$ (33)

where $r_{\rm ion}$ is the ionization radius. We see that the ratio (33)
decreases with increasing $r_{\rm ion}$ and maximal, when $r_{\rm ion}$
coincides with the event horizon. Thus, at $r_{\rm ion}=2$, the ratio (33) is
equal to unity. Finally, we write the expression for the energy of ionized
particle in the form

$E_{3}=\left(\frac{1}{\sqrt{2~{}r_{\rm
ion}}}+\frac{1}{2}\right)E_{1}+\frac{q_{3}Q}{r_{\rm ion}}.$ (34)

One can see that the energy of ionized particle is maximal when $q_{3}$ and
$Q$ have the same sign, which is also expected result – the charged particle
is accelerated due to the Coulombic repulsion force acting between the black
hole and particle. It is useful to define the ratio between the energies of
ionized and neutral particles, which would represent the efficiency of the
acceleration process. Writing the black hole mass and the speed of light
explicitly and substituting $q_{3}=Ze$ and $m_{1}\approx A\,m_{n}$, where $Z$
and $A$ are the atomic and mass numbers, $e$ is an elementary charge and
$m_{n}$ is the nucleon mass, we find

$\frac{E_{3}}{E_{1}}=\frac{1}{2}+\sqrt{\frac{GM}{2\,c^{2}\,r_{\rm
ion}}}+\frac{Z\,e\,Q}{A\,m_{n}\,c^{2}\,r_{\rm ion}}.$ (35)

The ionized particle is accelerated only when the right hand side of the
equation (35) is greater than unity. If the ionization point appears near the
event horizon, $r_{\rm ion}\approx 2GM/c^{2}$, then the condition
$E_{3}>E_{1}$ is satisfied for arbitrary positive values of the black hole
charge, $Q>0$. If the ionization point occurs at the ISCO radius, i.e. $r_{\rm
ion}=6GM/c^{2}$, for the energy of ionized particle to be greater than the
energy of infalling neutral particle, $E_{3}>E_{1}$, the charge of the black
hole has to satisfy the following condition

$Q\gtrsim 5.8\times 10^{11}\,\,\frac{A}{Z}\,\,\frac{M}{M_{\odot}}\,\,\,{\rm
Fr}\,,$ (36)

which is slightly greater than the lower limit of the estimated realistic
limits of the black hole charge, given by (5).

In Figure 3 (left) we plot the efficiency of the acceleration mechanism (ratio
of the energies of ionized and infalling particles with $Z/A=1$) against the
value of the black hole’s charge per solar mass at various ionization points.
The efficiency grows considerably with increasing the value of the black
hole’s charge and slightly decreases with increasing the distance between the
black hole and ionization point. In Figure 3 (right) we show the dependence of
the energy ratio on the ionization rate $Z/A$ for a fixed value of the black
hole charge. Let us estimate the maximal energy of ionized particle, which can
be accelerated by the non-rotating weakly charged black hole. From (35) and
using the uppermost realistic limit of the charge (5) we get

$E_{\rm ion}^{\rm max}\approx 1.01\times
10^{6}\,\,\,Z\,\,\frac{Q}{10^{18}\,{\rm Fr}}\,\,\frac{M_{\odot}}{M}\,\,{\rm
GeV},$ (37)

or, equivalently $\approx 1620$ erg. The ratio of energies of ionized and
neutral particles, in this case, is equal to $E_{\rm ion}^{\rm
max}/E_{n}\approx 10^{6}$. In sharp contrast to the magnetic Penrose process
[1, 2], where the energy of ionized particle is increasing with increasing the
black hole mass, in the case of a non-rotating weakly charged black hole the
energy of a charged particle is inversely proportional to the mass of the
black hole. Therefore, the maximal energy is determined by the limiting value
of the charge-to-mass ratio of the black hole $Q/M$ (see, the limits (5) and
the charge of the ionized particle $Ze$. This implies that the maximal energy
of ionized particle accelerated by the weakly charged non-rotating black hole
is similar for both stellar mass and supermassive black holes222Note that in
the weakly magnetized rotating black hole case [2], the energy of ionized
particle grows proportionally to the black hole’s mass and magnetic field
strength $E_{\rm ion}\sim BM$..

##  IV Possible application: high-energy cosmic rays

Figure 4: Constraints on the black hole mass and its charge to accelerate
protons of various energies ($1-10^{12}$GeV), marked by black lines. The red
line corresponds to the maximum theoretical limit of the Reissner-Nordström
charge. The blue vertical line shows the acceleration capability of an example
source, corresponding to a supermassive black hole Sgr A* located at the
Galactic centre, whose charge constrains are taken from [12].

In [1] it has been proposed that the ultra-efficient regime of the magnetic
Penrose process, which requires the presence of an external magnetic field in
the vicinity of rotating black hole can be relevant for the explanation of the
origin of the ultra-high-energy cosmic rays (UHECRs). UHECRs are the phenomena
composed from individual charged particles with detected energies exceeding
$10^{20}$eV, whose production mechanism remains unknown. In this section we
show that similar acceleration can be achieved in a more simplified setup of a
weakly charged non-rotating black hole. In Figure 4 we demonstrate the
constraints on the black hole mass and charge to serve as an accelerator of
charged particles (protons) of certain energy. As it is expected, increasing
the black hole charge for a given black hole mass one can reach the UHECR
orders of energies. A central charge of the black hole, in this case, is still
smaller than the maximal theoretical charge limit by many orders of magnitude.

As an example, in Figure 4 we depict the acceleration capability of the
Galactic centre supermassive black hole Sgr A*, whose charge has been
constrained in [12]. In particular, the constraint shows that Sgr A* is
capable to act as a PeVatron of charged particle with the energy of
accelerated protons being of the order of $10^{15}$eV. It is interesting to
note that this energy coincides with the knee of the cosmic ray spectrum,
above which the particle’s flux is considerably suppressed, suggesting the
change in the cosmic ray source. The result is also similar to that obtained
in the magnetic field and rotating black hole case [1], thus the presented
mechanism can serve as another cosmic ray acceleration mechanism and
alternative explanation of the cosmic ray knee when applied to the Galactic
supermassive centre black hole Sgr A*.

##  V Conclusions

In this paper we proposed a simple mechanism of acceleration of particles to
high-energies by the ionization of neutral particles in the vicinity of a non-
rotating Schwarzschild black hole carrying small electric charge, whose
gravitational effect on the spacetime metric is negligible. We started from
the description of the motion of charged particle and showed that the
effective potential in the case of a weakly charged black hole increases (or
decreases) with increasing (or decreasing) the electric interaction parameter
${\cal Q}$. We found that the innermost stable circular orbit (ISCO) of the
charged particle increases for both positive and negative values of the
parameter ${\cal Q}$. The results are in accord with previous similar studies
by [18, 11, 17].

We have found that the energy of ionized particle can be much greater than
initial energy of the neutral particle if both charges of the ionized particle
and the black hole have the same sign. Thus, similar acceleration process
occurring in the magnetized Kerr black hole spacetime and studied by [1, 19]
works also in the weakly charged non-rotating black hole case. For the
realistic upper limit of the black hole charge given by (5), which is at least
12 orders of magnitude smaller than the maximal theoretical Reissner-Nordström
limit (2), we have found that the charged particle with the charge $Ze$ can be
accelerated to the Lorentz $\gamma$-factor exceeding $Z\times 10^{6}$.

It is necessary to note that the energy of accelerated charged particles comes
in expense of the electrostatic energy of the black hole given by its electric
charge in contrast to many black hole energy extraction mechanisms, which use
the rotational energy of the black hole (see, e.g. [20, 16, 3, 21]). The
electric Penrose process studied in the current paper is similar to the
magnetic Penrose process, proposed in mid 1980’s in [3, 4, 22], where it is
the rotational energy of the black hole that is extracted. Connection of the
magnetic Penrose process with another famous black hole energy extraction
process – the Blandford-Znajek mechanism [16] was discussed in [23]. It has
been shown in [2] that the energy of accelerated ionized particle in the
magnetic Penrose process is proportional to the product of the magnetic field
strength and the black hole mass, i.e. $E_{\rm ion}\sim B\,M_{\rm BH}$.
However, in the electric Penrose process studied in the current paper, this
energy is proportional to the black hole’s charge and inversely proportional
to the black hole’s mass, i.e. $E_{\rm ion}\sim Q/M_{\rm BH}$. Therefore, the
maximal energy of ionized particle is restricted by the upper limit of the
charge-to-mass ratio of a black hole.

Further we also discussed possible astrophysical application of the presented
results for the production and acceleration of the ultra-high-energy cosmic
rays. In Figure 4 we presented the constraints on the maximum energy of the
cosmic ray particle accelerated by a non-rotating black hole with given mass
and charge. It has been shown that the ultra-high-energy acceleration does not
require the presence of a strong black hole charge, in a sense that the
Schwarzschild black hole spacetime is sufficient. We have applied the model to
the Galactic centre supermassive black hole Sgr A*, which is the best known
black hole candidate, and found that the energy of accelerated protons can
slightly exceed $10^{15}$eV, which coincides with the knee of the cosmic ray
spectrum. The presented results, however, are quite general and can
potentially operate also in neutron stars. It is especially interesting to
look for a similar acceleration scenarios in the electrospheres of pulsars.
Nevertheless, we leave this discussion for further studies.

Despite the similarities between the acceleration processes described in the
current paper and the magnetic Penrose process of particle acceleration by
rotating black hole in a magnetic field [1, 2], one of the main differences
between the two is that the motion of a charged particle in the former case is
always regular, while in the later case the chaotic behaviour of escaping
charged particle is usually observed [24, 19]. In the presence of external
magnetic field, it is expected that significantly larger number of charged
particles escape along the the directions given by magnetic field lines. The
character of the induced electric field around magnetized rotating Kerr black
hole has no spherical symmetry, it is rather of quadrupole character.
Meanwhile in the field of magnetized Schwarzschild black holes only
redirections in chaotic motion are observed [25], but no acceleration is
possible as no electric part of the field is induced. In a weakly charged non-
rotating black hole case, the combined gravitational and electric field is
spherically symmetric, therefore one would expect isotropic statistics of
escaping charged particles with no preferred direction of the motion. In
general, this can affect the statistics and interpretation of observed events,
being of the distinguishing observational signature of the electric Penrose
process.

###### Acknowledgements.

The authors would like to acknowledge the Research Centre for Theoretical
Physics and Astrophysics and Institute of Physics of Silesian University in
Opava for institutional support. Z. S. acknowledges the support of the Grant
No. 19-03950S of Czech Science Foundation (GACR).

## References

  * Tursunov  _et al._ [2020a] A. Tursunov, Z. Stuchlík, M. Kološ, N. Dadhich, and B. Ahmedov, Supermassive Black Holes as Possible Sources of Ultrahigh-energy Cosmic Rays, The Astrophysical Journal 895, 14 (2020a), arXiv:2004.07907 [astro-ph.HE] .
  * Tursunov and Dadhich [2019] A. Tursunov and N. Dadhich, Fifty Years of Energy Extraction from Rotating Black Hole: Revisiting Magnetic Penrose Process, Universe 5, 125 (2019), arXiv:1905.05321 [astro-ph.HE] .
  * Wagh  _et al._ [1985] S. M. Wagh, S. V. Dhurandhar, and N. Dadhich, Revival of the Penrose process for astrophysical applications, The Astrophysical Journal 290, 12 (1985).
  * Parthasarathy  _et al._ [1986] S. Parthasarathy, S. M. Wagh, S. V. Dhurandhar, and N. Dadhich, High efficiency of the Penrose process of energy extraction from rotating black holes immersed in electromagnetic fields, The Astrophysical Journal 307, 38 (1986).
  * Misner  _et al._ [1973] C. W. Misner, K. S. Thorne, and J. A. Wheeler,  _Gravitation_ (1973).
  * Stuchlik [1983] Z. Stuchlik, The Motion of Test Particles in Black-Hole Backgrounds with Non-Zero Cosmological Constant, Bulletin of the Astronomical Institutes of Czechoslovakia 34, 129 (1983).
  * Ruffini and Wilson [1975] R. Ruffini and J. R. Wilson, Relativistic magnetohydrodynamical effects of plasma accreting into a black hole, Phys. Rev. D 12, 2959 (1975).
  * Bally and Harrison [1978] J. Bally and E. R. Harrison, The electrically polarized universe., The Astrophysical Journal 220, 743 (1978).
  * Weingartner  _et al._ [2006] J. C. Weingartner, B. T. Draine, and D. K. Barr, Photoelectric Emission from Dust Grains Exposed to Extreme Ultraviolet and X-Ray Radiation, The Astrophysical Journal 645, 1188 (2006), arXiv:astro-ph/0601296 [astro-ph] .
  * Wald [1974] R. M. Wald, Black hole in a uniform magnetic field, Phys. Rev. D 10, 1680 (1974).
  * Zajacek and Tursunov [2019] M. Zajacek and A. Tursunov, The Electric Charge of Black Holes: Is It Really Always Negligible, The Observatory 139, 231 (2019), arXiv:1904.04654 [astro-ph.GA] .
  * Zajaček  _et al._ [2018] M. Zajaček, A. Tursunov, A. Eckart, and S. Britzen, On the charge of the Galactic centre black hole, Mon. Not. R. Astron Soc. 480, 4408 (2018), arXiv:1808.07327 [astro-ph.GA] .
  * Levin  _et al._ [2018] J. Levin, D. J. D’Orazio, and S. Garcia-Saenz, Black hole pulsar, Phys. Rev. D 98, 123002 (2018), arXiv:1808.07887 [astro-ph.HE] .
  * Tursunov  _et al._ [2020b] A. Tursunov, M. Zajaček, A. Eckart, M. Kološ, S. Britzen, Z. Stuchlík, B. Czerny, and V. Karas, Effect of Electromagnetic Interaction on Galactic Center Flare Components, The Astrophysical Journal 897, 99 (2020b), arXiv:1912.08174 [astro-ph.GA] .
  * Tursunov  _et al._ [2016] A. Tursunov, Z. Stuchlík, and M. Kološ, Circular orbits and related quasiharmonic oscillatory motion of charged particles around weakly magnetized rotating black holes, Phys. Rev. D 93, 084012 (2016), arXiv:1603.07264 [gr-qc] .
  * Blandford and Znajek [1977] R. D. Blandford and R. L. Znajek, Electromagnetic extraction of energy from Kerr black holes, Mon. Not. R. Astron Soc. 179, 433 (1977).
  * Hackstein and Hackmann [2020] J. P. Hackstein and E. Hackmann, Influence of weak electromagnetic fields on charged particle ISCOs, General Relativity and Gravitation 52, 22 (2020), arXiv:1911.07645 [gr-qc] .
  * Pugliese  _et al._ [2011] D. Pugliese, H. Quevedo, and R. Ruffini, Motion of charged test particles in Reissner-Nordström spacetime, Phys. Rev. D 83, 104052 (2011), arXiv:1103.1807 [gr-qc] .
  * Stuchlík  _et al._ [2020] Z. Stuchlík, M. Kološ, J. Kovář, P. Slaný, and A. Tursunov, Influence of Cosmic Repulsion and Magnetic Fields on Accretion Disks Rotating around Kerr Black Holes, Universe 6, 26 (2020).
  * Penrose [1969] R. Penrose, Gravitational Collapse: the Role of General Relativity, Nuovo Cimento Rivista Serie 1, 252 (1969).
  * Kološ  _et al._ [2021] M. Kološ, A. Tursunov, and Z. Stuchlík, Radiative Penrose process: Energy gain by a single radiating charged particle in the ergosphere of rotating black hole, Phys. Rev. D 103, 024021 (2021), arXiv:2010.09481 [gr-qc] .
  * Bhat  _et al._ [1985] M. Bhat, S. Dhurandhar, and N. Dadhich, Energetics of the Kerr-Newman black hole by the Penrose process, Journal of Astrophysics and Astronomy 6, 85 (1985).
  * Dadhich  _et al._ [2018] N. Dadhich, A. Tursunov, B. Ahmedov, and Z. Stuchlík, The distinguishing signature of magnetic Penrose process, Mon. Not. R. Astron Soc. 478, L89 (2018), arXiv:1804.09679 [astro-ph.HE] .
  * Stuchlík and Kološ [2016] Z. Stuchlík and M. Kološ, Acceleration of the charged particles due to chaotic scattering in the combined black hole gravitational field and asymptotically uniform magnetic field, European Physical Journal C 76, 32 (2016), arXiv:1511.02936 [gr-qc] .
  * Kološ  _et al._ [2015] M. Kološ, Z. Stuchlík, and A. Tursunov, Quasi-harmonic oscillatory motion of charged particles around a Schwarzschild black hole immersed in a uniform magnetic field, Classical and Quantum Gravity 32, 165009 (2015), arXiv:1506.06799 [gr-qc] .