# Ergodicity, lack thereof, and the performance of
reservoir computing with memristive networks

Valentina Baccetti  School of Science, RMIT University, Melbourne, Victoria
3000, Australia  Ruomin Zhu  Zdenka Kuncic  School of Physics, University of
Sydney, Sydney, NSW 2006, Australia  Francesco Caravelli  Theoretical Division
(T4), Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA

###### Abstract

Networks composed of nanoscale memristive components, such as nanowire and
nanoparticle networks, have recently received considerable attention because
of their potential use as neuromorphic devices. In this study, we explore the
connection between ergodicity in memristive and nanowire networks, showing
that the performance of reservoir devices improves when these networks are
tuned to operate at the edge between two global stability points. The lack of
ergodicity is associated with the emergence of memory in the system. We
measure the level of ergodicity using the Thirumalai-Mountain metric, and we
show that in the absence of ergodicity, two memristive systems show improved
performance when utilized as reservoir computers (RC). In particular, we
highlight that it is also important to let the system synchronize to the input
signal in order for the performance of the RC to exhibit improvements over the
baseline.

††preprint: APS/123-QED

##  I  Introduction

Memristive networks are electronic circuits that use memristive devices, a
type of resistive switching memory element. These networks store and recall
information by modifying the device’s resistance. They find applications in
computer memory, neuromorphic computing, and signal processing [1]. Memristive
networks offer benefits such as high density, low power consumption, and
potential for high-speed operation. Simultaneously, there is growing interest
in alternative approaches to computation and optimization in response to the
rapidly increasing demands on computing [2]. Various proposals have emerged to
address this challenge, some of which involve the use of oscillators or
frequency domain encoding [3, 4, 5, 6, 7], leveraging near- or in-memory
computation [8, 9, 10, 11, 12, 13], and exploring memcomputing [12, 13]. These
innovative techniques aim to provide more efficient solutions for complex
optimization problems [14, 5, 12, 15, 16, 3, 4, 17, 18]. However,
understanding large assemblies of memristive devices in non-equilibrium
statistical mechanics remains a challenge. Recent research has focused on
exploring the geometric and statistical properties of nanowire networks, where
electrical junctions exhibit memristive behavior due to the interplay between
tunneling and filament formation phenomena[19, 20].

Conductive nano-filament formation at nanowire-nanowire junctions creates a
memristive device, while quantum tunneling contributes to additional
nonlinearity in switching behavior, making the dynamics more complex [21]. As
most of the voltage drop occurs at the junctions, a basic model for the
conductance evolution of these networks involves an assembly of memristive
devices with voltage or current generators. The dynamic behavior of these
systems is currently being studied, especially regarding bias-induced
conductance transitions [22, 23, 24, 21].

Various conductance transitions have been observed in memristive devices,
often characterized by transient unstable dynamics of the memristive
components and their internal memory parameter. One well-known transition,
which has been defined as “rumbling”, has been analytically identified in the
simplest model of memristive networks [23]. A similar transition has also been
observed in nanowire networks [21]. This transition involves the system
effectively moving between different minima of an effective potential and is
marked by bursts of transient positive Lyapunov exponents. It arises from the
coexistence of multiple low-dimensional equilibrium points in the dynamics.
The mean-field theory for such systems can be ensured by mapping them to a
PEDS (projective embedding of dynamical systems) [25]. This is relevant to our
study as we aim to identify non-ergodic behavior associated with these
transitions.

In certain memristive systems, as the applied voltage increases, an effective
mean potential develops with multiple minima. As one minimum becomes dominant,
the system undergoes a rapid chaotic transition towards this emerging stable
fixed point. This transition is distinct from the conventional Landau picture
of symmetry breaking with bifurcation, as it arises from the competition
between two minima.

While the understanding of these transitions is clearer in simplified
memristive device models, it becomes less evident in more realistic systems.
In this paper, we aim to describe the dynamic behavior of these transitions in
two systems using ergodicity measures. Ergodicity, a concept in statistical
mechanics and thermodynamics, pertains to the long-term behavior of a system.
It states that the time average of a system over an extended period is
equivalent to the ensemble average of the system. In essence, it implies that
the system reaches a steady state, and its long-term average behavior aligns
with the average behavior across many instances of the same initial
conditions.

Ergodicity plays a crucial role in various natural sciences, such as physics
and chemistry, enabling an understanding of system evolution, equilibrium
attainment, and the prediction of long-term behavior based on short-term
observations. In the context of memory-dependent systems, the lack of ergodic
behavior signifies (hard) ergodicity breaking, which occurs when symmetry
breaking transpires in thermodynamic systems [26]. Typically, physical systems
operate within a regime where ergodicity holds for a subset of possible phase
space states. Quantifying ergodicity and gaining insight into the dynamical
state of a physical system is vital for comprehending the system’s operational
regime. One of the goals of this paper will also be to test the hypothesis
that computation near a transition point improves.

To test this hypothesis, we will use reservoir computing (RC) as our
computational model [27], which has been recently shown to be universal [28].
It has been reported for instance that the “edge of chaos” may be important
for the performance of RC [29], but it is important to stress that critical
states have been reported both for biological neuronal networks in the brain
[30] and in other artificial neural networks [31].

This article is organized as follows: in section II we introduce the general
definition of ergodicity and how it is related to memory; in section III we
introduce the two memristive systems we are considering in this study; in
section IV we give the definition of the TM metric for the two memristive
systems we are considering, and test the RC task results for both memristive
systems in terms of ergodicity breaking. Conclusions follow.

##  II Memory, ergodic Convergence and the Thirumalai-Mountain metric

Ergodicity and the emergence of memory are closely intertwined concepts within
the realms of statistical physics and complexity science. Ergodicity refers to
a system’s property of uniformly exploring all possible states in its phase
space over an extended period. In essence, it characterizes a system’s
behavior as it traverses its accessible states in a time-averaged manner. On
the other hand, memory denotes the influence of past states on the current
state of a system.

In many systems, the emergence of memory is closely connected to the violation
of ergodicity. When a system is ergodic, its behavior remains independent of
its past history, exhibiting no memory effect. However, when ergodicity is
violated, the system may display persistent or long-term correlations,
resulting in the emergence of memory.

A notable example of the relationship between ergodicity and memory can be
observed in spin-glasses, which are disordered magnetic systems that exemplify
complex dynamics with persistent or long-term correlations that foster memory
emergence. The behavior of spin glasses is often described through two-time
correlation functions, which measure correlations between spin configurations
at different points in time.

Studying two-time correlation functions in spin glasses sheds light on system
dynamics, including the emergence of memory and ageing phenomena. These
functions measure correlations between spin configurations at different times,
while the auto-correlation function reveals details about the system’s
relaxation dynamics. Slow decays observed in these functions signify the
violation of ergodicity and the emergence of memory, while the dependence on
both $t$ and the waiting time $t_{w}$ highlights the ageing phenomenon in spin
glasses [32] and other frustrated systems [33, 34]. Ergodicity refers to a
system’s property where the long-term behavior can be deduced from a single,
extended observation of the system, rather than relying on multiple
independent realizations.

According to Boltzmann’s hypothesis, the trajectories of any dynamical system
in its phase space eventually evolve into regions where macroscopic properties
reach thermodynamic equilibrium [35]. The ergodic hypothesis states that
ensemble averages and time averages coincide as time progresses. This means
that the ensemble-averaged value of an observable, denoted as $\langle
g\rangle$, can be obtained by averaging its values over time during the
observable’s evolution. Mathematically for the observable $g(t)$ we have:

$\displaystyle\langle
g\rangle=\lim_{T\to\infty}\frac{1}{T}\int_{{0}}^{T}g(t)dt.$ (1)

It is important to note that various definitions of ergodicity exist in the
Physics and Mathematics literature [36, 37, 38], with significant
distinctions. For example, in Markov chains, ergodic behavior requires a
strongly connected transition graph, which is rarely the case. In practice, a
system can exhibit ergodicity within a specific subset of its phase space. The
concept of “effective ergodicity” was introduced to describe scenarios where
the system rapidly samples coarse-grained regions [37].

A measure of effective ergodic convergence is based on the observation that
certain components of a system exhibit identical average characteristics at
thermal equilibrium. These characteristics are determined by an observable
defined on the system’s phase space, denoted as $\Gamma$. To assess the
effective ergodic behavior of the observable, it is necessary to estimate its
average value using an ensemble approach, such as a thermal ensemble. This
estimation is typically performed using the Thirumalai-Mountain (TM)
$g$-fluctuating metric, denoted as $\Omega_{g}(t)$.

The TM $g$-fluctuating metric, introduced by Thirumalai and Mountain [37, 39,
40], quantifies the difference between the ensemble-averaged value of the
observable, $g(\Gamma)$, and the sum of the instantaneous values of $g(t)$ for
each component of the system. At a given time $t$, the TM $g$-fluctuating
metric is expressed as:

$\displaystyle\Omega_{g}(t)=\frac{1}{N}\sum_{j=1}^{N}\big{[}\bar{g}_{j}(t)-\langle
g(t)\rangle\big{]}^{2},$ (2)

Here, $\bar{g}_{j}(t)$ represents the time-averaged value per component, and
$\langle g(t)\rangle$ denotes the instantaneous ensemble average, defined as:

$\displaystyle\bar{g}_{j}(t)=\frac{1}{t}\sum_{i=0}^{t}g_{j}(t_{i}),\ \langle
g(t)\rangle=\frac{1}{N}\sum_{j=1}^{N}g_{j}(t).$ (3)

This definition assumes that $g(\Gamma)$ serves as a suitable physical order
parameter, effectively characterizing the system’s behavior. The definition of
instantaneous ensemble average will depend on the considered system, as it
will be shown in Sec. IV.

The rate of ergodic convergence can be quantified using the derivative of the
effective ergodic convergence, denoted as $\Omega_{g}^{{}^{\prime}}$. This
quantity is given by:

$\displaystyle\Omega_{g}^{{}^{\prime}}=\frac{\Omega_{G}(t)}{\Omega_{G}(0)}\to\frac{1}{tD_{g}},$
(4)

in the diffusive regime. When the system is instead in a sub-diffusive or
super-diffusive regime, the power of the relaxation in time changes from 1 to
a different exponent. Coarse graining of the phase space leads to the
clustering of the system’s accessible states, making the concept of effective
ergodicity more applicable. Effective ergodicity is achieved when the system
uniformly explores the coarse-grained regions within a finite time [39]. Here,
$D_{G}$ represents the diffusion coefficient associated with the property
being studied, and $\Omega_{g}$ refers to the effective ergodic convergence.
This definition aligns with other notions of ergodicity in cases where
diffusion follows a power law [41]. The rate of ergodic convergence,
determined using the TM metric, provides an estimate of the system’s
ergodicity. The behavior of the rate, as described by Eq. (4), indicates the
system’s attainment of effective ergodicity. For example, if the inverse of
the rate scales linearly with time, the system reaches ergodicity in a
diffusive manner: $1/{\Omega}_{g}\rightarrow D_{G}t$, where
$D_{G}=\Omega_{G}(0)$ represents the diffusion coefficient of the property
$G$. It is generally expected that the rate scales with time as
$\Omega^{\prime}_{G}(t)\sim t^{-p}$.

When the inverse of $\Omega_{g}^{{}^{\prime}}$ exhibits a linear relationship
with time, it indicates that all points in the phase space are equally likely,
resembling the behavior of Brownian motion. This approach has been applied in
various contexts, such as simple liquids [37], earthquake fault networks [42,
43], and the Ising model [44].

##  III Models of memristive devices

We would like to introduce the notion of memristive device that we will use in
the following. We will consider two models of memristive device, a resistive
current-controlled memristive device, which can be described by the equations

$\displaystyle V(t)$ $\displaystyle=$ $\displaystyle R\big{(}x\big{)}I(t),$
(5) $\displaystyle\frac{dx}{dt}$ $\displaystyle=$ $\displaystyle f(I,x).$

The second model we consider is a conductance based device, of the form

$\displaystyle I(t)$ $\displaystyle=$ $\displaystyle
G\big{(}\lambda\big{)}V(t),$ (6) $\displaystyle\frac{d\lambda}{dt}$
$\displaystyle=$ $\displaystyle f(V,\lambda).$ (7)

Both models satisfy the pinched hysteresis property for the $I-V$ curve,
characteristic of memristive devices [45].

###  III.1 Memristive network toy model

In [46], the dynamical equation for a circuit of memristors was derived under
the assumption of a resistance current-controlled device. The resistance
function used, $R(x)=R_{on}(1-x)+xR_{off}$, approximates $TiO_{2}$ memristors,
where $R_{on}$ and $R_{off}$ represent the limiting resistances, and
$x\in[0,1]$ is internal memory parameter, the state variable describing the
size of the oxygen-deficient conducting layer. At the lowest order, the
evolution of the internal memory parameter can be described by a simple
equation with hard boundaries:

$\frac{dx}{dt}=\frac{R_{off}}{\beta}I-\alpha
x=\frac{R_{off}}{\beta}\frac{V}{R(x)}-\alpha x$ (8)

Here, $\alpha$ and $\beta$ are the decay constant and the effective activation
voltage per unit of time, respectively, which determine the timescales of the
system.

While the model presented above provides a simple description of a polar
resistive device, various extensions have been explored in the literature. For
instance, to account for diffusive effects near the boundaries, some models
remove the hard boundaries and introduce a window function [47, 48]. Although
these models better capture the detailed IV curves of physical devices, they
still exhibit the fundamental pinched hysteresis behavior observed in the
linear model.

In adimensional units ($\tau=\alpha t$), the equation for $x(t)$ in a single
memristor device under an applied voltage $S$ can be derived using Ohm’s law
($S=RI$), [49, 23]. The resulting equation is:

$\frac{d}{d\tau}x=\frac{S}{\alpha\beta}\frac{1}{1-\chi
x}-x=-\partial_{x}V(x,s),$ (9)

Here, $\chi=\frac{R_{off}-R_{on}}{R_{off}}$ and $s=\frac{S}{\alpha\beta}$,
with $0\leq\chi\leq 1$ in relevant physical cases. The dynamics of the system,
represented by a single memristor device, are fully characterized by the
gradient descent in the effective potential given by [23]:

$V(x,s)=\frac{1}{2}x^{2}+\frac{s}{\chi}\log(1-\chi x),$ (10)

Figure 1: Top: Evolution of the mean-field potential for the toy model as a
function of voltage. Bottom: Phase diagram of the toy model, and the
appearance and disappearance of stable points. The intermediate region in
voltage between $Vc$ and $Vc^{*}$ is characterized by the coexistence of two
stable points. The buffer region corresponds to a metastable state which acts
as the boundary between the attracting stable points for a single memristive
device, while the shading represents the tunneling region for the network of
devices.

The potential exhibits two minima separated by a barrier, as depicted in Fig.
1 (top), for $s=0.15$ and $s=0.24$ with $\chi=0.9$. The range of $s$ for the
existence of a barrier is limited, and when $\chi$ approaches 1, the local
minimum can move inside the domain $[0,1]$, leading to the emergence of an
unstable fixed point (i.e., the peak of the barrier). Consequently, two basins
of attraction and locally stable minima are formed. A pictorial phase diagram
illustrating this behavior is shown in Fig. 1 (bottom), highlighting the
critical voltage points at which such behavior occurs. Below a critical
voltage point $V_{c}$, only one fixed point is present. At the $V_{c}$ value,
a new fixed point emerges, but it becomes metastable in the presence of noise
(dashed red line). At an intermediate point $Vc<V_{m}<Vc^{*}$, the two
metastable points have equal energy, representing the switching point where
the original stable fixed point becomes metastable. At higher values of
$V=Vc^{*}$, the first fixed point disappears by merging with the barrier and
becoming flat. Ultimately, only one fixed point remains at higher values.

For a network of memristors of this type, the mean-field potential resembles
exactly the single memristor defined above. For details, see Supp. Mat. A. The
key difference is that the mean-field potential is only an approximation when
the system is high dimensional, e.g. now there is an effective metastable
region, as shown in Fig. 1, instead of two stable minima, when one of the
minima is lower than the other. Such effective “tunneling” can be explained
using the theory developed in [25], i.e. the fact that local mean field maxima
become saddle points in high dimensions. Here we will use this fact to study
effective ergodicity breaking in memristive networks.

In fact, these two regimes can be characterized via the Thirumalai-Mountain
(TM) metric introduced before. In Fig. 2 we show the TM metric as a function
of time for $\bar{s}=0.18$ (top), $\bar{s}=0.24$ (center), and $\bar{s}=0.25$
(numerically calculated for a circuit of N = 50 memristors with parameters
$\alpha=\beta=1$ and $\chi=0.9$). As we can see, in the former case the TM
metric relaxes as a power law with an approximate intermediate exponent of
$p\in[-2,-1.5]$, typical of a diffusive regime. For $\bar{s}=0.24$, we see
instead a typical behavior of weak ergodicity breaking, e.g. a transient non-
monotonicity of the TM metric, associated with the trajectories effectively
“tunneling” through the mean field barrier. This behavior is indeed transient
as for $\bar{s}=0.25$ the TM metric relaxes again as a power law.

Using the toy model, we can then pinpoint such non-monotonicity of the TM
metric to a transient transition between two stable asymptotic states and the
effective symmetry breaking of the potential. Thus, we can immediately
identify the symmetry breaking of the potential as the culprit of such
transient non-ergodic behavior.

Figure 2: Thirumalai-Mountain metric as a function of time and amplitude,
calculated for the stochastic memristive model with noise, with the sample
average calculated using the mean field value introduced in (13). We observe
that, in the regime in which the potential has a single minimum ($V<Vc$), the
TM metric decays rapidly. Instead, for higher voltages in which two minima are
present ($Vc<V<Vc^{*}$), the TM metric does not decay. This is a symptom of
non-ergodic behavior.

###  III.2 Nanowire networks

We now consider a more realistic model of memristive networks, which can be
associated with self-assemblies of silver nanowires. Via established bottom-up
self-assembly techniques, one can readily synthesize nanowire networks (NWNs)
[22, 50]. These NWNs typically have a 2D spatial distribution of randomly
oriented nanowires that are interconnected by cross-point MIM junctions. The
NWN we consider have densities of 10 junctions/$\mu$m2 and 0.5
nanowires/$\mu$m2. Device electrodes can be deposited onto the substrate using
a mask, from which the conductance measurements can be readily performed. This
bio-inspired structure is difficult to design and fabricate using top-down
techniques. Self-assembled NWNs exhibit topological properties, such as small-
world propensity and modularity, that are similar to biological neural
networks and distinct from random and grid-like networks [51]. Unlike fully
connected bipartite networks in artificial neural networks, small-world
networks have local connectivity and short path lengths, making them
relatively sparse. Although small-worldness is necessary for important
functional properties, such as synchronizability and information flow, it
alone cannot explain the diverse range of dynamics across networks that
exhibit this structural property.

A model for the simulation of realistic NWNs has been introduced previously in
the literature [52, 24, 21]. Fig. 3 (a) shows a visualization of a simulated
nanowire network containing 1000 nanowires and 6877 junctions. Self-assembly
is modeled by distributing nanowires on a $3\times 3\,\mu{\rm m}^{2}$ 2D
plane, with their centers uniformly sampled from $[0,3]$ and orientation
uniformly sampled from $[0,\pi]$. The lengths of the nanowires are sampled
from a gamma distribution (mean = $100\,$nm, standard deviation $10\,$nm),
based on experimental measurements [22]. In theoretical studies, as
illustrated in Fig. 3 (b), the NWN is transformed to the corresponding
graphical representation, where the nodes represent nanowires and the edges
are the cross-point junctions.

Figure 3: Example of a nanowire network generated with the random wire model.
(a) Simulated NWN with 1024 nanowires and 6877 junctions. (b) Graphical
representation of the NWN in (a).

In this work, all simulation results for nanowire networks are generated using
a network comprised of 1000 nanowires and 6877 junctions. All variables,
except the adjacency matrix $A$, are time-dependent. A model for the
conductance of a single junction, associated with the filament length is
provided in the supplementary material, see Sec. A. Each junction evolves as
voltage bias is continuously applied to the network. The modified nodal
analysis approach is applied to the graphical representation to solve
Kirchhoff’s voltage and current conservation laws at each time step [53]. This
is equivalent to the method used for the derivation of the exact network
equation for the toy model, eqn. (17). Although the NWN model is based on
polymer-coated Ag nanowires, with memristive junction internal dynamics that
differ from that of the toy model (based on metal-oxide memristors), the
network dynamics are similar and one should think of the two models as
equivalent from a physical perspective.

Figure 4: (a) Collective conductance of the network between the source and
drain nodes with a DC input of varying pulse width. Dashed red line indicates
the edge of formation of a high conductance current pathway between the source
and drain nodes. (b) Visualization of network activation levels for different
pre-initilization pulse widths $T_{0}$.

For the purpose of using a NWN as a reservoir, a Mackey-Glass time-series
signal with delay parameter $\tau=17$ is delivered to a source electrode as
the input voltage signal. Before implementing the time-series prediction task
using RC, a DC input of varying duration is applied to the NWN to initialize
the internal state of the network and prepare it for RC. We refer to this pre-
initialization protocol as “priming the system” [54]. Fig. 4(a) shows the
reservoir’s conductance (blue curve) as a function of the DC input length
$T_{0}$. The shaded region represents the general conductance transition
regime, identified from previous studies [21, 24], and the dashed line at
$T_{0}=2.17\,$s represents when the first conductance pathways form between
the source and drain nodes. The internal state of the network for different
$T_{0}$ is visualized in Fig. 4(b). In cases where the reservoir is under-
activated ($T_{0}<2\,$s), the majority of memristive components remain
inactive, resulting in insufficient dynamics from the network. When the
reservoir is over-activated ($T_{0}>8\,$s), the internal dynamics of the
system become saturated, limiting the system’s capacity to process additional
information. The conductance transition regime from Fig. 4(a) corresponds to
an intermediate dynamical state of the reservoir, where conductance paths
first span the network and the internal state of the system produces dynamical
features that are more diverse than at other activation times.

##  IV Ergodicity, reservoir computing and bistability

###  IV.1 Toy model

Reservoir computing using [55] was studied for the first time in [56], among
other passive circuits. A brief recap of the procedures behind reservoir
computing is provided in Supp. Mat. B.

Here, we have used a similar scheme to the reservoir, using a memristive
circuit comprising $N=50$ idealized memristors.

We have considered parameter values $\alpha=\beta=1$, and $\chi=0.9$, for
which the system experiences a symmetry break, and where we expect its
dynamics to become non-ergodic, see III.1. This is the regime where the
reservoir is at the edge of stability, as defined in [29], where in general,
reservoirs can have optimal performance, although not always guaranteed. Weak
ergodicity breaking in a dynamical system is associated with a strongly
chaotic regime [41]. We will come back to this point later. We have also added
a small noise value $\sigma=0.01$. The equations (18) were numerically
integrated with a time step $dt=0.1$.

For the input signal $V(t)$ we have used a Mackey-Glass time series with delay
parameter $r=0.2$, $\gamma=0.1$, $\tau=17$, and $n=10$. Our input signal was
given by

$\displaystyle S(t)=\mathrm{b}+\mathrm{a}V(t),$ (11)

with $a\in[0,10]$ a multiplicative factor of the input $V(t)$ and
$b\in[-0.4,1]$ a parameter bias. An example of the input signal is shown in
Fig. 6 (top). The parameter $b$ represents the bias input to the network,
while $a$ is the amplitude.

For each value of $a$ and $b$, we performed a Mackey-Glass reconstruction task
using the internal memory states as the dynamical variables for the RC. The
rMSE, which we use as a measure of the performance of the RC, is shown in Fig.
6. We can see from the figure that varying $a$ leads to a reduced value of the
rMSE, right where the rumbling transition is present, approximately for
$b=0.2$.

The difference in the quality of the task can be observed, for the Mackey-
Glass time series, in Fig. 5, for the toy model described in App. A.

Figure 5: Reservoir computing prediction task with the toy model, for $b=0.2$,
and $a=0.01$ (a) and $a=10$ (b)respectively. As we can see, the regime of
optimal prediction corresponds both the regime in which the system’s
Thirumalai-Mountain metric does not converge to zero.

At a value of $b=0.2$, a plot of the rMSE is shown in Fig. 7. As we can see,
the value of the rMSE decreases as a function of $a$, meaning that the larger
the magnitude of the input signal the better the performance, but this happens
only when $b$ is carefully chosen near the transition point. The reason why
this occurs can be inferred by analyzing the response of the internal memory
values, which are reported in the Appendix, as a function of the parameter $b$
and for $a=0.2$, below, near, and above the transition point. When $a$ is such
that the potential has a single minimum, the value of the memory oscillates in
the vicinity of that minimum. When however $a$ is such that the potential has
two minima, e.g. in the symmetry-breaking regime, at sufficiently large values
of $b$ the memory values start to oscillate between the two minima. This
implies that the time evolution of the system is much more dynamic in the
symmetry-broken phase, and this effectively results in an improvement in the
performance of the RC.

Figure 6:  MSE of the Mackey-Glass reconstruction as a function of the
parameters $a$ (Amp) and $b$ (bias). Figure 7: (a): Reservoir computing
prediction error as a function of the amplitude of the input signal, while the
bias is fixed on the transition point. We can see that rMSE decreases as a
function of the amplitude. (a): Thirumalai-Mountain metric as a function of
time for various amplitudes. We see that for increasing amplitudes, the metric
ceases to converge to zero, indicating an effective non-ergodic behavior. In
the intermediate regime, we see oscillations due to the fact that the system
is jumping from one minimum to the other.

As we can see, increasing the value of the amplitude leads to the internal
memory values fluctuating more prominently between the two stable states. This
is also shown in Fig. 7 (b), where we plot the TM metric evaluated on the
response of the system in the small $Amp=1$ and large $Amp=10$ amplitude
regimes. In one first case, the system is still ergodic, while in the case
$Amp=10$ the TM metric is not converging to a value zero.

Furthermore, we have calculated the $TM_{x}$ metric near the transition bias
$a=0.2$ for varying values of the amplitude $b$ as

$\displaystyle\Omega_{x}(t)=\frac{1}{N}\sum^{N}_{j}\left[\bar{x}_{j}(t)-\langle
x(t)\rangle\right]^{2}.$ (12)

As the ensemble average $\langle x(t)\rangle$ we have consider the mean field
value $x_{cg}$ defined as [23]

$\displaystyle x_{cg}=\frac{1}{N}\sum^{N}_{ij}{\mathcal{P}}_{ij}x_{j}.$ (13)

We have chosen the mean field value $x_{cg}$ as a natural order parameter that
resembles, in definition, the general meaning of ensemble average.

This is actually one of the advantages of studying this model first, as in
this case, we know the details of the order parameter for the whole memristive
network.

###  IV.2 Nanowire networks

The nanowire connectome dynamics has been extensively explored across various
studies. Prior research has highlighted that nanowire Networks (NWNs) showcase
brain-like dynamics, demonstrating their optimal information storage and
processing capabilities at conductance transition points [21, 24, 23]. More
recently, a dynamical mean-field theoretical technique for polymer-coated Ag
nanowires has uncovered emergent dynamical features [57] such as transitions.
In the context of the two-terminal setup, these transitions are commonly
observed within a regime termed as the ’edge of formation’ [21]. This ’edge of
formation’ is delineated by the activation of memristive components, but
precedes an exponential surge in the formation of parallel paths. As depicted
in Fig. 4, this regime establishes a select few high-conductance current paths
between the two electrodes (a phenomenon known as the ’winner-takes-all’)
[58]. Within this state, the internal dynamics of the NWN intricately map the
input signal to a diverse feature space.

In the case of the toy model, previous analytical studies provide a
comprehensive understanding of both the system dynamics and the order
parameters to be employed [49, 23], which allows us to apply a simplification
for the TM metric. Nevertheless, the same technique cannot be utilized for the
nanowire model since it is more realistic and cannot be characterized in the
same way. For that reason, the collective conductance of the nanowire network
is used as the order parameter, which is determined based on the individual
conductances of individual memristive junctions and the underlying circuitry
shaped by the connectivity. In the meantime, the findings derived from the toy
model will play an important role in interpreting the results of nanowires.

####  IV.2.1 Thirumalai-Mountain metric

We now wish to link computation performance and ergodicity, also in the case
of NWNs, we used the TM metric to understand the effective ergodicity of these
dynamical systems; in particular, eqn. 12 is also used to calculate the TM
metric for NWNs. However, the meaning associated with each quantity is closer
to the definition of the metric used in the original work on supercooled
liquid, [37]. The overall observable of interest is the two-probe conductance,
defined earlier for nanowires, which is a scalar. Then, the time average is
calculated using the global conductance, while the collective conductance is
employed as an ensemble average across various realizations of the initial
conditions:

$\displaystyle\bar{g}(t)$ $\displaystyle=\frac{1}{T}\int_{0}^{T}G^{i}(t)dt,$
(14) $\displaystyle\langle g(t)\rangle$
$\displaystyle=\frac{1}{N}\sum_{i}G^{i}(t),$ (15)

where $G^{i}(t)$ is the two-probe conductance of a particular realization at
time $t$ and $G(t)$ is the collective and effective conductance of the network
between the two points. For the time average quantity, note that we select a
single element of the ensemble to evaluate the time average.

The ensemble of different realizations is generated by randomly perturbing the
filament levels of junctions in the network:

$\displaystyle\Lambda^{i}=(1+\delta)\Lambda_{0},$ (16)

where $\delta$ is randomly sampled from a flat distribution, while
$\Lambda^{i}$ is the parameter controlling the length of the filaments for all
junctions, while $\Lambda_{0}$ is the initial condition of the filament. The
random variable $\delta$ is sampled from a uniform distribution over
$(-0.1,0.1)$. Thus, effectively we are sampling over the initial conditions of
the system.

The TM metric, as a function of the bias, is shown in Fig. 8 (b). As we can
see, and consistently with the case of the toy model, for small and large
values of the bias, the metric decays. At intermediate values of the bias, at
the edge of the transition between low and high conductance states, where the
conductance synchronizes with the input, the TM metric fails to decay. This is
exactly the regime in which the conductance transition occurs and is a
signature of a special state for the nanowire network, highly synchronized
with the input.

Figure 8:  (a) Nanowire system simulated conductance, using the Mackey-Glass
time series as input, for various values of the bias voltage $b$. As we can
see, for small values and larger values of the bias, the conductance decays
and grows to the asymptotic value. At intermediate values of the bias, the
conductance oscillates. (b) Thirumalai-Mountain metric for the effective
conductance of the nanowire network.

####  IV.2.2 Reservoir Computing

For the realistic nanowire network model above, the implementation of RC
involved the two probe conductance, with designated input and readout nodes.
Such construction has already been considered in the literature [54, 59].

The fitting task of the Mackey-Glass time series was performed using NWNs
under the RC framework, with two nanowire nodes in the network selected as the
source and the drain. The MG signal was linearly transformed as described in
Eq. 11 and delivered to the network as input, while the same MG signal with 5
steps ahead was employed as the target. The task can be breakdown into three
phases:

  1. 1.

Priming: A 2 $V$ DC input of varying length $T_{0}$ was applied to drive the
internal state of the network. The first 1000 data points of the MG time-
series were delivered subsequently to wash out the influence from the initial
state.

  2. 2.

Training: The effective conductance of the network corresponding to
$t=1000-4000$ ms was measured and multiplexed using the virtual node technique
[60] to provide training features (see details in Appendix). The readout layer
was trained using a linear regression with a ridge parameter $r=0.01$.

  3. 3.

Testing: The effective conductance during $t=4000-5000$ ms was measured and
multiplexed in the same fashion as the training phase, while the trained
readout layer was applied to make predictions. The performance of the
reservoir can thus be evaluated accordingly.

We consider two regimes: one in which the system is initialized away from the
voltage transition point between the high and low conductance state, and one
in which the system sits at the boundary between the two. As we can see from
Fig. 9 (top), the RMSE of the RC model has a behavior very similar to the one
observed for the toy model of Fig. 6. At values of the amplitude close to the
transition point, in which the Thirumalai-Mountain fails to converge to zero,
the physical RC performance peaks. Our intuition is that the system tuned at
the point in which it is at the edge of a transition, is more prone to
synchronization [29] with the input signal. This, combined with existing
knowledge [54] on the average “priming” time that it takes for the system to
reach synchronization with the signal, shows that the combination of input
time and choice of input voltage leads to optimal performance. We can see the
difference between the tuned and non-tuned network in Fig. 10.

Figure 9: RMSE of RC prediction task with the simulated nanowire system as a
function of driving. Top: RMSE as a function of the bias and the amplitude of
the input signal. As we can see, the minimum is located exactly at the
numerically observed transition point in bias. Bottom: $\log$RMSE as a
function of the priming time $T0$ and the amplitude of the input. We see that
the optimal results are obtained after a certain initial time $T0\approx 2$ s.

.

Figure 10: Fitting result for different biases, for amp = 4.9 V, $T_{0}$ = 2
s. As we can see, near the transition the performance of the RC prediction
task increases dramatically. The red dashed line represents the divide between
training and test sets.

.

##  V Conclusions

Memory effects are a critical aspect of many physical and biological systems,
and they have been shown to play a vital role in the behavior of complex
systems. Meanwhile, ergodicity is a property of systems that describes the
degree to which they explore their phase space. In recent years, in
particular, physical systems with memory such as nanowires or nanoparticle
connectomes, memristive, and other nanoscale devices, have become increasingly
essential candidates as substrates for synthetic intelligent devices, e.g.
brain-like physical material. In particular, there are strong indications,
both in theoretical models and experiments, that these devices exhibit
conductance transitions both of the first and second order. These transitions
have been linked to the edge of chaos, a concept that refers to the boundary
between ordered and chaotic behavior in complex systems.

The present study explores the interplay between ergodicity breaking and
memory in two models of memristive devices. The first model we studied was a
toy model introduced in the literature to understand analytically the
properties of purely memristive networks and has been instrumental in
understanding their non-equilibrium properties such as Lyapunov functions [61]
or conductance transitions [23]. The second is a more realistic model for
memristive networks composed of nanowires [21, 62, 24, 54], used to understand
the properties of polymer-coated Ag nanowires [22, 63, 64]. In both cases, we
studied the Thirumalai-Mountain [37] metric to understand how the systems
relax when driven by different inputs, and in particular the ergodic
properties of these systems. Both in the case of the toy model, for which the
conductance transition has been studied analytically using a variety of
techniques [65, 66, 23, 67, 25], and in the case of the nanowire model [21],
we found that the Thirumalai-Mountain metric signals lack of ergodicity near
the voltage bias value where the conductance transitions are expected to be.

This is, in particular, interesting in view of the fact that it has been
recently observed that edge of instability can be linked to a computational
advantage [68, 29]. A similar result is observed in this paper. We did observe
that, in particular, this is not necessarily true unless the dynamical system
under study is synchronized to the input signal, as suggested in [29]. In
fact, this improved performance occurs only after a transient period
(“priming”) which allows the nanowire network to synchronize to the input
signal, as shown in Fig. 9. The results described above are consistent across
two different models: the toy model in which conductance transitions can be
understood quantitatively and analytically, and the more realistic nanowire
model able to capture the experimentally observed conductance in Ag nanowires.
Thus, these results suggest that there might be an underlying common theory to
explain these transitions.

In conclusion, by connecting memory effects, ergodicity, and the edge of
chaos, we have identified a set of principles that can be used to create more
effective computational models. Our research has shown for the first time that
non-ergodic behavior can be linked to the effectiveness of reservoir
computing, leading to new approaches for developing more advanced and
efficient computational tools. We believe that our findings will inspire
further investigations into the connections between memory effects,
ergodicity, and the edge of chaos, and will lead to new and exciting
developments in the field of machine learning and beyond. In particular, in
future work we will discuss the application of the ideas developed in this
paper to meta-plasticity with memristive systems, [69, 70], in particular in
view of the recent results on the meta-plasticity of nanowire networks shown
in [64].

###### Acknowledgements.

The work of FC was carried out under the auspices of the NNSA of the U.S. DoE
at LANL under Contract No. DE-AC52-06NA25396, and in particular support from
LDRD via 20230338ER and 20230627ER. RZ is supported by a Postgraduate Research
Excellence Award scholarship from the University of Sydney. VB acknowledges
funding through the RMIT Vice-Chancelor’s Research Fellowship.

## References

  * Mehonic  _et al._ [2020] A. Mehonic, A. Sebastian, B. Rajendran, O. Simeone, E. Vasilaki, and A. J. Kenyon, Memristors—from in-memory computing, deep learning acceleration, and spiking neural networks to the future of neuromorphic and bio-inspired computing, Advanced Intelligent Systems 2, 2000085 (2020).
  * Oliver [2019] W. D. Oliver, Quantum computing takes flight, Nature 574, 487 (2019).
  * Böhm  _et al._ [2019] F. Böhm, G. Verschaffelt, and G. V. der Sande, A poor man’s coherent ising machine based on opto-electronic feedback systems for solving optimization problems, Nature Comm. 10, 10.1038/s41467-019-11484-3 (2019).
  * Pierangeli  _et al._ [2019] D. Pierangeli, G. Marcucci, and C. Conti, Large-scale photonic ising machine by spatial light modulation, Phys. Rev. Lett. 122, 213902 (2019).
  * Vadlamani  _et al._ [2020] S. K. Vadlamani, T. P. Xiao, and E. Yablonovitch, Physics successfully implements lagrange multiplier optimization, Proc. of the Nat. Aca. of Sci. 117, 26639 (2020).
  * Csaba and Porod [2020] G. Csaba and W. Porod, Coupled oscillators for computing: A review and perspective, App. Phys. Rev. 7, 011302 (2020).
  * Goto and et al. [2021] H. Goto and et al., High-performance combinatorial optimization based on classical mechanics, Sci. Adv. (2021).
  * Singh and et. al. [2019] G. Singh and et. al., Near-memory computing: Past, present, and future, Micro. and Micro. 71, 102868 (2019).
  * Ielmini and Wong [2018] D. Ielmini and H.-S. P. Wong, In-memory computing with resistive switching devices, Nature Ele. 1, 333 (2018).
  * Traversa and et. al. [2014] F. L. Traversa and et. al., Dynamic computing random access memory, Nanotechnology 25, 285201 (2014).
  * Sebastian and et. al. [2020] A. Sebastian and et. al., Memory devices and applications for in-memory computing, Nature Nano. 15, 529 (2020).
  * Traversa and Di Ventra [2015] F. L. Traversa and M. Di Ventra, Universal memcomputing machines, IEEE Trans. on Neur. Net. and Learn. Sys. 26, 2702 (2015).
  * Di Ventra and Traversa [2018] M. Di Ventra and F. L. Traversa, Perspective: Memcomputing: Leveraging memory and physics to compute efficiently, J. of App. Phys. 123, 180901 (2018).
  * Hennessy and Patterson [2019] J. L. Hennessy and D. A. Patterson, A new golden age for computer architecture, Comm. of the ACM 62, 48 (2019).
  * Sutton and et. al. [2017] B. Sutton and et. al., Intrinsic optimization using stochastic nanomagnets, Sci. Rep. 7, 10.1038/srep44370 (2017).
  * Torrejon and et. al. [2017] J. Torrejon and et. al., Neuromorphic computing with nanoscale spintronic oscillators, Nature 547, 428 (2017).
  * Kirkpatrick  _et al._ [1983] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671 (1983).
  * Dorigo and Stützle [2004] M. Dorigo and T. Stützle,  _Ant Colony Optimization_ (MIT, 2004).
  * Milano  _et al._ [2019] G. Milano, S. Porro, I. Valov, and C. Ricciardi, Recent developments and perspectives for memristive devices based on metal oxide nanowires, Advanced Electronic Materials 5, 1800909 (2019).
  * Kuncic and Nakayama [2021] Z. Kuncic and T. Nakayama, Neuromorphic nanowire networks: principles, progress and future prospects for neuro-inspired information processing, Advances in Physics: X 6, 1894234 (2021).
  * Hochstetter and et al. [2021] J. Hochstetter and et al., Avalanches and edge-of-chaos learning in neuromorphic nanowire networks, Nature Comm. 12, 10.1038/s41467-021-24260-z (2021).
  * Diaz-Alvarez  _et al._ [2019] A. Diaz-Alvarez, R. Higuchi, P. Sanz-Leon, I. Marcus, Y. Shingaya, A. Z. Stieg, J. K. Gimzewski, Z. Kuncic, and T. Nakayama, Emergent dynamics of neuromorphic nanowire networks, Scientific Reports 9, 14920 (2019).
  * Caravelli  _et al._ [2021] F. Caravelli, F. Sheldon, and F. L. Traversa, Global minimization via classical tunneling assisted by collective force field formation, Science Advances 7, 1542 (2021), https://www.science.org/doi/pdf/10.1126/sciadv.abh1542 .
  * Zhu  _et al._ [2021] R. Zhu, J. Hochstetter, A. Loeffler, A. Diaz-Alvarez, T. Nakayama, J. T. Lizier, and Z. Kuncic, Information dynamics in neuromorphic nanowire networks, Scientific Reports 11, 13047 (2021).
  * Caravelli  _et al._ [2023a] F. Caravelli, F. Traversa, M. Bonnin, and F. Bonani, Projective embedding of dynamical systems: Uniform mean field equations, Physica D: Nonlinear Phenomena 450, 133747 (2023a).
  * Palmer [1982] R. Palmer, Broken ergodicity, Advances in Physics 31, 669 (1982), https://doi.org/10.1080/00018738200101438 .
  * Jaeger and Haas [2004] H. Jaeger and H. Haas, Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication, Science 304, 78 (2004).
  * Grigoryeva and Ortega [2018] L. Grigoryeva and J.-P. Ortega, Echo state networks are universal, Neural Networks 108, 495 (2018).
  * Carroll [2020] T. L. Carroll, Do reservoir computers work best at the edge of chaos?, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 10.1063/5.0038163 (2020), 121109, https://pubs.aip.org/aip/cha/article-pdf/doi/10.1063/5.0038163/14106721/121109_1_online.pdf .
  * Tagliazucchi  _et al._ [2012] E. Tagliazucchi, P. Balenzuela, D. Fraiman, and D. R. Chialvo, Criticality in large-scale brain fMRI dynamics unveiled by a novel point process analysis, Frontiers in Physiology 3, 10.3389/fphys.2012.00015 (2012).
  * Morales and Muñoz [2021] G. B. Morales and M. A. Muñoz, Optimal input representation in neural systems at the edge of chaos, Biology 10, 702 (2021).
  * Vincent [2007] E. Vincent, Ageing, rejuvenation and memory: The example of spin-glasses, in  _Ageing and the Glass Transition_, edited by M. Henkel, M. Pleimling, and R. Sanctuary (Springer Berlin Heidelberg, Berlin, Heidelberg, 2007) pp. 7–60.
  * Saccone  _et al._ [2022] M. Saccone, F. Caravelli, K. Hofhuis, S. Parchenko, Y. A. Birkhölzer, S. Dhuey, A. Kleibert, S. van Dijken, C. Nisoli, and A. Farhan, Direct observation of a dynamical glass transition in a nanomagnetic artificial hopfield network, Nature Physics 18, 517 (2022).
  * Saccone  _et al._ [2023] M. Saccone, F. Caravelli, and e. al, Real-space observation of ergodicity transitions in artificial spin ice, Nature Communications 14 (2023).
  * Dorfman [1999] J. R. Dorfman,  _An Introduction to Chaos in Nonequilibrium Statistical Mechanics_ (Cambridge university Press, 1999).
  * Ma [1985] S.-K. Ma,  _Statistical Mechanics_ (World Scientific, 1985).
  * Thirumalai  _et al._ [1989] D. Thirumalai, R. D. Mountain, and T. R. Kirkpatrick, Ergodic behavior in supercooled liquids and in glasses, Phys. Rev. A 39, 3563 (1989).
  * Petersen [1983] K. Petersen,  _Ergodic theory_ (World Scientific, 1983).
  * Mountain and Thirumalai [1989] R. D. Mountain and D. Thirumalai, Measures of effective ergodic convergence in liquids, The Journal of Physical Chemistry 93, 6975 (1989), publisher: American Chemical Society.
  * Thirumalai and Mountain [1993] D. Thirumalai and R. D. Mountain, Activated dynamics, loss of ergodicity, and transport in supercooled liquids, Phys. Rev. E 47, 479 (1993).
  * Bel and Barkai [2006] G. Bel and E. Barkai, Ergodicity breaking in a deterministic dynamical system, Europhys. Lett. 74, 15 (2006).
  * Tiampo  _et al._ [2003] K. F. Tiampo, J. B. Rundle, W. Klein, J. S. S. Martins, and C. D. Ferguson, Ergodic dynamics in a natural threshold system, Phys. Rev. Lett. 91, 238501 (2003).
  * Tiampo  _et al._ [2007] K. F. Tiampo, J. B. Rundle, W. Klein, J. Holliday, J. S. Sá Martins, and C. D. Ferguson, Ergodicity in natural earthquake fault networks, Phys. Rev. E 75, 066107 (2007).
  * Süzen [2014] M. Süzen, Effective ergodicity in single-spin-flip dynamics, Phys. Rev. E 90, 032141 (2014).
  * Chua [1971] L. Chua, Memristor-the missing circuit element, IEEE Transactions on Circuit Theory 18, 507 (1971).
  * Caravelli [2017a] F. Caravelli, Locality of interactions for planar memristive circuits, Phys. Rev. E 96, 052206 (2017a).
  * Joglekar and Wolf [2009] Y. N. Joglekar and S. J. Wolf, The elusive memristor: properties of basic electrical circuits, Eur. J. of Phys. 30, 661 (2009).
  * Prodromakis and et. al. [2011] T. Prodromakis and et. al., A versatile memristor model with nonlinear dopant kinetics, IEEE Trans. on El. Dev. 58, 3099 (2011).
  * Caravelli  _et al._ [2017] F. Caravelli, F. L. Traversa, and M. Di Ventra, Complex dynamics of memristive circuits: Analytical results and universal slow relaxation, Phys. Rev. E 95, 022140 (2017).
  * Li  _et al._ [2020] Q. Li, A. Diaz-Alvarez, R. Iguchi, J. Hochstetter, A. Loeffler, R. Zhu, Y. Shingaya, Z. Kuncic, K.-i. Uchida, and T. Nakayama, Dynamic electrical pathway tuning in neuromorphic nanowire networks, Advanced Functional Materials 30, 2003679 (2020).
  * Loeffler  _et al._ [2020] A. Loeffler, R. Zhu, J. Hochstetter, M. Li, K. Fu, A. Diaz-Alvarez, T. Nakayama, J. M. Shine, and Z. Kuncic, Topological properties of neuromorphic nanowire networks, Frontiers in Neuroscience 14, 10.3389/fnins.2020.00184 (2020).
  * Kuncic  _et al._ [2020] Z. Kuncic, O. Kavehei, R. Zhu, A. Loeffler, K. Fu, J. Hochstetter, M. Li, J. M. Shine, A. Diaz-Alvarez, A. Stieg, J. Gimzewski, and T. Nakayama, Neuromorphic Information Processing with Nanowire Networks, in  _2020 IEEE International Symposium on Circuits and Systems (ISCAS)_ (2020) pp. 1–5.
  * Nilsson and Riedel [2011] J. W. Nilsson and S. Riedel,  _Electric Circuits (9th ed),_ (Pearson Education, Saddle River NJ, 2011).
  * Fu  _et al._ [2020] K. Fu, R. Zhu, A. Loeffler, and e. al, Reservoir computing with neuromemristive nanowire networks, Proceedings of the International Joint Conference on Neural Networks (IJCNN) 20006228 (2020).
  * Caravelli and Sheldon [2020] F. Caravelli and F. C. Sheldon, Phases of memristive circuits via an interacting disorder approach, https://arxiv.org/abs/2009.00114 (2020).
  * Sheldon  _et al._ [2022] F. C. Sheldon, A. Kolchinsky, and F. Caravelli, The computational capacity of lrc, memristive and hybrid reservoirs, Phys. Rev. E 106 (2022).
  * Caravelli  _et al._ [2023b] F. Caravelli, G. Milano, C. Ricciardi, and Z. Kuncic, Mean field theory of self-organizing memristive connectomes, Annalen der Physik 535, 10.1002/andp.202300090 (2023b).
  * Manning  _et al._ [2018] H. Manning, F. Niosi, C. da Rocha, and e. al., Emergence of winner-takes-all connectivity paths in random nanowire networks, Nat Commun 9 (2018).
  * Daniels  _et al._ [2022] R. Daniels, J. Mallinson, Z. Heywood, P. Bones, M. Arnold, and S. Brown, Reservoir computing with 3D nanowire networks, Neural Networks 154, 122 (2022).
  * Appeltant  _et al._ [2011] L. Appeltant, M. Soriano, G. V. der Sande, J. Danckaert, S. Massar, J. Dambre, B. Schrauwen, C. Mirasso, and I. Fischer, Information processing using a single dynamical node as complex system, Nature Communications 2, 10.1038/ncomms1476 (2011).
  * Caravelli [2019] F. Caravelli, Asymptotic behavior of memristive circuits, Entropy 21, 789 (2019).
  * Zhu  _et al._ [2020] R. Zhu, J. Hochstetter, A. Loeffler, A. Diaz-Alvarez, A. Stieg, J. Gimzewski, T. Nakayama, and Z. Kuncic, Harnessing adaptive dynamics in neuro-memristive nanowire networks for transfer learning, in  _2020 International Conference on Rebooting Computing (ICRC)_ (2020) pp. 102–106.
  * Milano  _et al._ [2022] G. Milano, G. Pedretti, K. Montano, and e. al., In materia reservoir computing with a fully memristive architecture based on self-organizing nanowire networks, Nat. Mater. , 195–202 (2022).
  * Loeffler  _et al._ [2023] A. Loeffler, A. Diaz-Alvarez, R. Zhu, N. Ganesh, J. M. Shine, T. Nakayama, and Z. Kuncic, Neuromorphic learning, working memory, and metaplasticity in nanowire networks, Science Advances 9, 10.1126/sciadv.adg3289 (2023).
  * Bartolucci  _et al._ [2020] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, Inversion-free leontief inverse: statistical regularities in input-output analysis from partial information, https://arxiv.org/abs/2009.06350 (2020).
  * Bartolucci  _et al._ [2023] S. Bartolucci, F. Caccioli, F. Caravelli, and P. Vivo, Ranking influential nodes in networks from aggregate local information, Physical Review Research 5, 10.1103/physrevresearch.5.033123 (2023).
  * Caravelli [2023] F. Caravelli, Cycle equivalence classes, orthogonal weingarten calculus, and the mean field theory of memristive systems, arXiv:2304.14890 (2023).
  * Langton [1990] C. G. Langton, Computation at the edge of chaos: phase transitions and emergent computation, Physica D 42, 12 (1990).
  * Chialvo and Bak [1999] D. Chialvo and P. Bak, Learning from mistakes, Neuroscience 90, 1137 (1999).
  * Carbajal  _et al._ [2022] J. P. Carbajal, D. A. Martin, and D. R. Chialvo, Learning by mistakes in memristor networks, Physical Review E 105, 10.1103/physreve.105.054306 (2022).
  * Caravelli [2017b] F. Caravelli, The mise en scéne of memristive networks: effective memory, dynamics and learning, Int. J. of Par., Em. and Dist. Sys. 33, 350 (2017b).
  * Bollobás [1998] B. Bollobás,  _Modern Graph Theory_ (Springer New York, 1998).
  * Sheldon  _et al._ [2020] F. C. Sheldon, F. Caravelli, and C. Coffrin, Fully analog memristive circuits for optimization tasks: a comparison, To appear in Handbook of Unconventional Computing, Ed. A. Adamatzky (2020).
  * Pershin  _et al._ [2013] Y. V. Pershin, V. A. Slipko, and M. Di Ventra, Complex dynamics and scale invariance of one-dimensional memristive networks, Physical Review E 87, 022116 (2013).
  * Sillin  _et al._ [2013] H. O. Sillin, R. Aguilera, H.-H. Shieh, A. V. Avizienis, M. Aono, A. Z. Stieg, and J. K. Gimzewski, A theoretical and experimental study of neuromorphic atomic switch networks for reservoir computing, Nanotechnology 24, 384004 (2013).
  * Simmons [1963] J. G. Simmons, Generalized formula for the electric tunnel effect between similar electrodes separated by a thin insulating film, Journal of Applied Physics 34, 1793 (1963).
  * Bellew  _et al._ [2015] A. T. Bellew, H. G. Manning, C. Gomes da Rocha, M. S. Ferreira, and J. J. Boland, Resistance of Single Ag Nanowire Junctions and Their Role in the Conductivity of Nanowire Networks, ACS Nano 9, 11422 (2015).

##  Appendix A Memristive models

###  $TiO_{2}$ toy model

For a circuit with multiple memristors having varying resistances $R(x_{i})$
and constant voltage generators $S_{i}=\tilde{S}$ connected in series, the
equation (9) for a single memristor extends to a system of coupled nonlinear
ordinary differential equations. This system describes the network dynamics of
the memory elements $x_{i}(t)$ and can be expressed as [49, 61]:

$\displaystyle\frac{d}{dt}\vec{x}=\frac{1}{\beta}(I-\chi{\mathcal{P}}X)^{-1}{\mathcal{P}}\vec{S}-\alpha\vec{x},$
(17)

Here, $\chi=\frac{R_{off}-R_{on}}{R_{off}}<1$ and
$X_{ij}(t)=x_{i}(t)\delta_{ij}$. The matrix ${\mathcal{P}}$ represents the
projection operator on the vector space of cycles in the circuit’s graph
$\mathcal{G}$ [49], and it can be obtained from the directed incidence matrix
$B$ of $\mathcal{G}$ as ${\mathcal{P}}=I-B(B^{t}B)^{-1}B^{t}$ [46, 71, 53,
72]. The property of ${\mathcal{P}}$ being a projection operator reflects
Kirchhoff’s circuit laws. In this paper, we consider $B$ as a random matrix to
abstract the dynamical system from any specific circuit topology.

It is interesting to see that there is a range of voltage in which the system
experiences bistability. This is not a pitchfork bifurcation, but it is
nonetheless the emergence of two (meta)stable points in the range
$[Vc,Vc^{*}]$ of Fig. 1 (bottom). We can see the relevance of this feature in
the study of ergodicity as a function of voltage, as we explain below.

In order to study the ergodic vs. non-ergodic behavior of such toy model, we
extend the system to a stochastic one, introducing noise. The dynamical
equations of a memristive device can incorporate noise through stochastic
terms, introducing randomness into the system’s behavior. This stochasticity
is manifested in the memory component of the device. This is a diffusive
memristive model of the form

$\displaystyle d\vec{x}$
$\displaystyle=\frac{1}{\beta}\left[I-\chi{\mathcal{P}}X\right]^{-1}{\mathcal{P}}\vec{S}dt-\alpha\vec{x}dt$
$\displaystyle\ +\sigma(\vec{x},t)d\vec{\eta},$ (18)

the vector $d\vec{\eta}$ contains $n$ Wiener stochastic processes with
$\sigma$ the diffusion coefficient.

It is known that for the memristive toy model of eqn. (17) and in the limit of
dense networks, for values of $\chi=0.9$ and
$\bar{s}=\frac{1}{\alpha\beta}\frac{1}{N}\sum_{i=1}^{N}S_{i}\approx 0.24$ in
which the system develops two mean-field minima, induces bulk dynamical
transition characterized by a short transient of positive Lyapunov exponents
[23]. For smaller or larger values of $\bar{s}$ instead, the memory parameters
are in a laminar regime, characterized by the presence of only a global
minimum in the potential $V(x,s)$.

We discuss the properties of the effective potential. For $\alpha=\beta=1$ and
$\chi=0.9$, the range is $1/10<s<5/18$. The emergence of the two minima is
non-perturbative requiring $\chi$ to be close to 1. In the case of $\chi=0$,
the network consists of regular resistors, and the behavior of the single-
element dynamical system is characterized by simple basins of attraction
without any exotic features. To ensure the validity of the equations of motion
for the single variables, cutoff functions, such as
$\frac{d}{d\tau}x=-W(x)f(x,y)$, are applied, where $W(x)=1$ for
${0<x<1}\cup{x=1,f(x,y)>0}\cup{x=0,f(x,y)<0}$, and zero otherwise. An implicit
analytical solution for this differential equation has been derived in [73].

### Junction Model

Nanowire-nanowire cross-points are modeled as electrically insulating,
ionically conducting junctions [74, 75, 22, 52]. The junctions exhibit
voltage-threshold memristive switching due to electro-chemical metallization
and the influence of electron tunnelling transport [21]. However, the model
neglects junction-level stochastic fluctuations, which introduces further
nonlinearities and lead to more complex dynamics [21]. Junction conductance,
$G=G(\lambda)$, is determined by a state variable $\lambda(t)$ which
characterizes the length of the conducting filament responsible for memristive
switching:

$\displaystyle
G(\lambda)=\frac{1}{R_{\mathrm{t}}(\lambda)+R_{\mathrm{on}}+\xi}+\frac{1}{R_{\mathrm{off}}},$
(19)

where

$\displaystyle\xi=\frac{R_{\mathrm{on}}^{2}}{R_{\mathrm{off}}-R_{\mathrm{on}}},$
(20)

and the tunnelling conductance $G_{t}(\lambda)$ is calculated using the low
voltage Simmon’s formula (Eq. 21) for a MIM junction [76]:

$\displaystyle G_{t}(\lambda)$ $\displaystyle=$
$\displaystyle[R_{\mathrm{t}}(\lambda)]^{-1}$ (21) $\displaystyle=$
$\displaystyle\frac{3(2m_{*})^{1/2}e^{5/2}(\phi/e)^{1/2}}{A2h^{2}s^{2}}$
$\displaystyle\hskip
2.84544pt\cdot\exp\left[-\frac{4\pi(2m_{*}e)^{1/2}}{h}s\left(\frac{\phi}{e}\right)^{1/2}\right]$
$\displaystyle s$ $\displaystyle=$
$\displaystyle\max\left[\left(\Lambda_{\mathrm{crit}}-|\Lambda|\right)\frac{s_{\mathrm{max}}}{\Lambda_{\mathrm{crit}}}\;,\;0\right],$
(22)

with effective mass $m_{*}=0.99m_{\mathrm{e}}$ and PVP layer (assumed
homogeneous) thickness $s_{\mathrm{max}}=5\,\text{nm}$. The difference between
Fermi levels of PVP and Ag determines the potential barrier
$\phi=0.82\,\text{eV}$. $A=(0.41\,\text{nm})^{2}=0.17\,\text{nm}^{2}$ is the
area of a face of the silver unit cell. The resistances of nanowires are
considered to be negligible compared to the resistance of the junctions [77].
Consequently, $G_{t}(\lambda)$ exhibits an additional nonlinear dependence on
$V$, through the filament growth parameter $s=s(\Lambda(V))$ which modulates
junction switching due to filament formation.

All junctions start with a high resistance “off” state initially. When the
value of $\lambda$ of a single junction reaches the predefined threshold
$\lambda_{\text{crit}}$, the resistance of the junction switches to “on”
state. The ratio of these resistance states is
$R_{\text{off}}/R_{\text{on}}=10^{3}$, with $R_{\text{on}}=10^{4}$ Ohms.

The evolution of $\lambda(t)$ is described by a polarity-dependent voltage-
threshold model [74, 52, 24, 21]:

$\dfrac{d\lambda}{dt}=\begin{cases}(|V(t)|-V_{\text{set}})\text{sgn}[V(t)],&|V(t)|>V_{\text{set}}\\\
0,&V_{\text{reset}}<|V(t)|<V_{\text{set}}\\\
b(|V(t)|-V_{\text{reset}})\text{sgn}[\lambda(t)],&|V(t)|<V_{\text{reset}}\end{cases}$
(23)

where $V_{\text{set}}$ is the on-threshold and $V_{\text{reset}}$ is the off-
threshold, and $b$ is a positive constant defining the relative rates of decay
of the filament. All simulation results presented in this work are generated
using the following parameters: $V_{\text{set}}=10^{-2}$V,
$V_{\text{reset}}=10^{-3}$V and $b=10$. Experimental validation of this model
is presented elsewhere [22, 21].

More specifically, the junction voltage distribution $V$ across the network is
obtained numerically by solving

${\cal L}^{\dagger}V=I\qquad,$ (24)

in which $I$ is current and $\cal{L}^{\dagger}$ is the expanded graph
Laplacian of the network, expressed in block matrix representation as

$\cal L^{\dagger}=\left[\begin{array}[]{c|c}\cal L&C\\\ \hline\cr C&0\\\
\end{array}\right]\qquad,$ (25)

where $\cal L$ is the graph Laplacian and where the elements of $C$ are either
1 if the nanowire node is connected to an external electrode or 0 otherwise.
The Laplacian $\cal L$ is generated by

${\cal L}=D-W\qquad,$ (26)

where $W$ is the weighted adjacency matrix of the network, with the weights
representing junction conductance distribution:

$W_{ij}=A_{ij}G(i,j)\qquad,$ (27)

where $G(i,j)$ is conductance on the edge connecting nodes $i$ and $j$ and $D$
is the weighted degree matrix generated from $W$:

$D=\textbf{diag}(d_{i})\qquad,\qquad d_{i}=\sum\limits_{k=1}^{N}W_{ik}\qquad.$
(28)

After Kirchoff laws are solved, the parameters $\lambda(t)$ are evolved
accordingly.

##  Appendix B Reservoir computing and multiplexing

Reservoir computing is a popular approach to machine learning that utilizes a
fixed, randomly connected network of nodes, known as a reservoir, to generate
temporal patterns in response to an input signal. The reservoir acts as a
dynamical system that maps input signals to high-dimensional feature spaces,
where a readout layer can then be trained to classify or predict the desired
output. This approach is particularly effective for processing time-varying
data and has been applied to a range of tasks, such as speech recognition,
natural language processing, and control systems. One of the key advantages of
reservoir computing is that the reservoir can be pre-designed and does not
require tuning during the training phase, reducing the complexity of the
learning process. Furthermore, the reservoir can be implemented using a
variety of physical substrates, such as electronic circuits, optical systems,
and biological neural networks, making it a versatile approach that can be
adapted to a wide range of applications.

The training of a reservoir computer involves finding a set of optimal readout
weights that map the reservoir states to the desired output. This can be
formulated as a linear regression problem, where the readout weights
$\mathbf{W}_{out}$ are estimated by minimizing the mean squared error between
the target output $\mathbf{y}(t)$ and the estimated output
$\hat{\mathbf{y}}(t)$:

$\min_{\mathbf{W}{out}}\frac{1}{T}\sum_{i=1}^{T}|\mathbf{y}(t_{i})-\mathbf{W}_{out}\mathbf{x}(t_{i})|^{2}$

where $\mathbf{x}(t_{i})$ is the reservoir state at time step $n$, and $N$ is
the total number of training samples. The solution to this optimization
problem can be obtained using the pseudo-inverse of the reservoir state matrix
$\mathbf{X}$:

$\mathbf{W}_{out}=(\mathbf{X}^{t}\mathbf{X}+\alpha\mathbf{I})^{-1}\mathbf{X}^{t}\mathbf{Y}$

with $\ {}^{t}$ representing the vector or matrix transposition, where
$\mathbf{Y}$ is the target output matrix, and $\alpha$ is a regularization
parameter that prevents overfitting. The reservoir state matrix is defined as:

$\mathbf{X}=[\mathbf{x}(1),\mathbf{x}(2),...,\mathbf{x}(N)]^{t}$

and the target output matrix is defined as:

$\mathbf{Y}=[\mathbf{y}(1),\mathbf{y}(2),...,\mathbf{y}(N)]^{t}$

Once the readout weights are estimated, they can be used to generate the
output for new input signals by computing:

$\hat{\mathbf{y}}(t_{i})=\mathbf{W}_{out}\mathbf{x}(t_{i})$

where $\hat{\mathbf{y}}(t_{i})$ is the estimated output at time step $t$. In
the following, we will assume that $\mathbf{x}(t)$ are either the internal
memory states or the conductance states.

In the case of a single conductance measured, in order to multiply the number
of virtual nodes, we used the multiplexing technique introduced in [60]. The
time series is time multiplexed at certain sampling frequencies. Then, the
time series is parsed into $N$ vectors of length $K$, such that $NK$ is the
length of the time series. Then, elements of such vectors at locations $mN+l$
for integer $l\in[1,K]$ are considered as different internal virtual nodes.
The regression then follows the same scheme as the previous section.

Reservoir computing has been successfully applied in various physical systems,
including opto-electronic and nanoscale systems such as silver nanowire
networks (NWNs) to implement neuromorphic computing. NWNs are particularly
interesting for this application due to their memristive behavior and complex
network topology that resembles that of biological neural networks. By
leveraging the dynamics of NWNs as a reservoir, the input signals can be
mapped to a high-dimensional state space through a random projection, followed
by a simple readout operation to obtain the desired output. This allows for
the implementation of complex computational tasks such as pattern recognition,
time-series prediction, and control in a highly efficient and parallel manner.
In this way, reservoir computing offers a promising approach to harness the
capabilities of NWNs for neuromorphic computing, and silver nanowire learning
represents an exciting area of research at the intersection of nanotechnology
and machine learning.

To be considered a feasible reservoir, a dynamical system must satisfy the
condition that its state approaches a nontrivial function of the input
trajectory in the long time limit. This can be expressed in terms of two
requirements[56]:

_Fading Memory_ : If the system were started from two different initial
conditions while driven with the same input trajectory $u$, the system’s
trajectories should eventually converge to the same state. This implies that
the system has a finite temporal memory;

_State Separation_ : Different input sequences should drive the system into
different trajectories. That is, if the same initial condition were driven
with two different input trajectories, the resulting reservoir trajectories
must be sufficiently different. These properties could be proven analytically
for the case of the memristive reservoir toy model.

These conditions ensure that the state of the reservoir becomes a function of
the input trajectory and that this function carries information about the
input trajectory. In the case of memristive systems, the memristive junctions
act as the dynamic elements of the reservoir, and the network dynamics arise
from the collective behavior of the junctions. The input signal drives the
memristive junctions to switch between high and low resistance states,
resulting in a temporal sequence of network states. The final state of the
network can be used as a high-dimensional feature vector for further
processing. In the case of silver nanowire networks, the nanowires themselves
act as the memristive elements, and their collective behavior gives rise to
the network dynamics. The input signal drives the nanowires to switch between
metallic and insulating states, resulting in a temporal sequence of network
states that can be used for further processing.

As we saw above, reservoir computing requires the fading memory property for
the dynamics of a reservoir, which is equivalent to the system tending to some
degree of effective ergodicity when expressed in physical terms. The scope of
this manuscript is to investigate fading memory and understand how quickly a
reservoir “forgets” its initial conditions, using the Thirumalai-Mountain
metric as our benchmark. The goal is to determine the optimal level of memory
capacity for different system parameters. However, there are limitations to
using ergodicity and the ergodic theorem, as ensemble averages are only
defined for systems at equilibrium. Nonetheless, we can still consider the
case of systems locally at equilibrium.

The original formulation of the TM metric aimed to predict the dynamics of a
system numerically simulated on a computer without having to calculate the
Lyapunov exponents of systems with a high number of variables. This is
particularly relevant when dealing with large-scale systems, where calculating
Lyapunov exponents can be computationally expensive.

Let us now briefly introduce the time series we used for our benchmarks. The
Mackey-Glass time series is a widely used benchmark for testing the
performance of time-series prediction models. It is a nonlinear,
nonautonomous, and nonperiodic time series that exhibits chaotic dynamics. The
series is defined by the following difference equation:

$\frac{dx(t)}{dt}=r\frac{x(t-\tau)}{1+x(t-\tau)^{n}}-\gamma x(t),$ (29)

where $x(t)$ is the value of the time series at time $t$, $\beta$, $\gamma$,
$n$, and $\tau$ are parameters that control the dynamics of the system. The
Mackey-Glass time series exhibits a range of dynamic behaviors, depending on
the values of these parameters. For example, for $r=0.2$, $\gamma=0.1$,
$\tau=17$, and $n=10$, the time series exhibits chaotic dynamics with a
positive Lyapunov exponent.