A mean-field model of memristive circuit interaction

TL;DR

This paper presents an exactly solvable mean-field model of interacting memristors, analyzing its dynamics, fixed points, and asymptotic behavior using a Lyapunov functional analogous to a long-range Ising Hamiltonian.

Contribution

It introduces a novel mean-field approach to model interacting memristors with an exactly solvable circuit and a Lyapunov functional for analyzing non-equilibrium dynamics.

Findings

Exact predictions for decay to lower resistance state

Reasonable predictions for decay to higher resistance state

Lyapunov functional resembles a long-range Ising Hamiltonian

Abstract

We construct an exactly solvable circuit of interacting memristors and study its dynamics and fixed points. This simple circuit model interpolates between decoupled circuits of isolated memristors, and memristors in series, for which exact fixed points can be obtained. We introduce a Lyapunov functional that is found to be minimized along the non-equilibrium dynamics and which resembles a long-range Ising Hamiltonian with non-linear self-interactions. We use the Lyapunov functional as an Hamiltonian to calculate, in the mean field theory approximation, the average asymptotic behavior of the circuit given a random initialization, yielding exact predictions for the case of decay to the lower resistance state, and reasonable predictions for the case of a decay to the higher resistance state.

Equations

$$\frac{d}{dt}w(t)=\alpha w(t)-\frac{R_{on}}{\gamma }I=\alpha w(t)-\frac{R_{on}}{\gamma }\frac{S}{R(w)}$$

$$\frac{R(w^{\ast })}{R_{on}}{w}^{\ast }=\left((1-w^{\ast })+\frac{R_{off}}{R_{on}}{w}^{\ast }\right)w^{\ast }=\frac{S}{\alpha \gamma }.$$

$$n{R}_{int}{i}_0(t)+\sum _{k=1}^n{R}_{ref}(i_0(t)-i_k(t))$$

$$R\left(w_k(t)\right)i_k+R_{ref}(i_k(t)-i_0(t))$$

$$\alpha w_k(t)-\frac{R_{on}}{\gamma }{i}_k(t)$$

$$\frac{d{w}_k}{dt}=\alpha w_k-\sum _{j=1}^n\frac{R_{on}}{\gamma }(\mathcal{I}-\frac{1}{n}M\mathcal{J}{)}_{kj}^{-1}\frac{S_j}{R_{ref}+R(w_j)},$$

$$\frac{d}{dt}{w}_i=\alpha w_i-f_i(\vec{w}),$$

$$i_k$$

$$i_0$$

$$\frac{S_k}{\alpha \gamma }$$

$$S_c=\alpha \gamma \left(\frac{R_{off}+R_{ref}}{R_{on}}-\frac{R_{ref}^2}{R_{on}(R_{int}+R_{ref})}\right).$$

$$\frac{S}{\alpha \gamma }$$

$${w^{\ast }}_{\pm }=-\frac{\chi }{2\xi }\pm \sqrt{\frac{{\chi }^2}{4{\xi }^2}+\frac{S}{\xi \alpha \gamma }}=-\frac{1}{2\xi }\left(1\pm \sqrt{1+4\frac{S\xi }{\alpha \gamma }}\right)$$

$$H(w_i)=\frac{\rho }{2n}\sum _i{w}_i^2-\frac{\rho }{n}\sum _{i,j}{w}_i{w}_j-\sum _i\frac{S_i}{\alpha \gamma }{w}_i+\sum _{i}E(w_i)$$

$$P(w(t=\infty )=1)$$

$$P(w(t=\infty )=0)$$

$$P(w_{in}<w^{\ast })$$

$$\psi ={\text{arg sup}}_{w\in [0,1]}\left(\left(\rho \psi +\frac{S}{\alpha \gamma }\right)w-E(w)\right).$$

$${\partial }_S\left(\text{max}\left(0,w^{\ast }(S)\right)\right)=\{\begin{array}{ll}\frac{1}{\xi \alpha \gamma }\frac{1}{2\sqrt{\frac{(\chi -\rho {)}^2}{4{\xi }^2}+\frac{S}{\xi \alpha \gamma }}}& S<0\\ 0& S>0\end{array}$$

$${\partial }_{S}w(t\gg 1)\approx \frac{R_{on}}{2\gamma \sqrt{\alpha (R_{off}-R_{on})}}{e}^{-\alpha t},$$

$$i_k=\frac{S_k}{R_{ref}+R(w_k)}+\frac{R_{ref}}{R_{ref}+R(w_k)}{i}_0$$

$$i_0=\frac{R_{ref}}{R_{ref}+R_{int}}\frac{1}{n}\sum _k^n{i}_k$$

$$i_k=\frac{S_k}{R_{ref}+R(w_k)}+\frac{R_{ref}}{R_{ref}+R(w_k)}\frac{R_{ref}}{R_{ref}+R_{int}}\frac{1}{n}\sum _j^n{i}_j$$

$$i_k=\sum _{j=1}^n(\mathcal{I}-\frac{1}{n}M\mathcal{J}{)}_{kj}^{-1}\frac{S_j}{R_{ref}+R(w_j)}$$

$$\frac{d{w}_k}{dt}=\alpha w_k-\sum _{j=1}^n\frac{R_{on}}{\gamma }(\mathcal{I}-\frac{1}{n}M\mathcal{J}{)}_{kj}^{-1}\frac{S_j}{R_{ref}+R(w_j)}$$

$$\left(\begin{array}[]{ccccc}n(R_{int}+R_{ref})&-R_{ref}&-R_{ref}&\cdots&-R_{ref}\\ -R_{ref}&R(w_{1})+R_{ref}&0&\cdots&0\\ -R_{ref}&0&R(w_{2})+R_{ref}&0&\vdots\\ \vdots&\vdots&\ddots&\ddots&0\\ -R_{ref}&0&\cdots&0&R(w_{n})+R_{ref}\end{array}\right)\left(\begin{array}[]{c}i_{0}\\ i_{1}\\ \vdots\\ \vdots\\ i_{n}\end{array}\right)=\left(\begin{array}[]{c}0\\ S_{1}\\ \vdots\\ \vdots\\ S_{n}\end{array}\right)$$

$$M=\left(\begin{array}[]{ccccc}a_{0}&-b&-b&\cdots&-b\\ -b&a_{1}&0&\cdots&0\\ -b&0&a_{2}&0&\vdots\\ \vdots&\vdots&\ddots&\ddots&0\\ -b&0&\cdots&0&a_{n}\end{array}\right)$$

$$(M^{-1}{)}_{ij}=\frac{C_{ij}}{D}$$

$$\displaystyle C_{12}=C_{21}=-det\left(\begin{array}[]{cccccc}a_{0}&-b&-b&-b&-b&-b\\ -b&0&0&0&0&0\\ -b&0&a_{3}&0&0&0\\ -b&0&0&a_{4}&0&0\\ -b&0&0&0&a_{5}&0\\ -b&0&0&0&0&a_{6}\\ \end{array}\right)$$

$$\left(\begin{array}[]{c}i_{0}\\ i_{1}\\ \vdots\\ \vdots\\ i_{n}\end{array}\right)=\frac{1}{D}C\left(\begin{array}[]{c}0\\ S_{1}\\ \vdots\\ \vdots\\ S_{n}\end{array}\right)\rightarrow\left(\begin{array}[]{c}i_{1}\\ \vdots\\ \vdots\\ i_{n}\end{array}\right)=\frac{1}{D}\left(\begin{array}[]{ccccc}q_{1}&b^{2}\ c_{12}&\ b^{2}\ c_{13}&\cdots&\ b^{2}\ c_{1n}\\ b^{2}\ c_{12}&q_{2}&b^{2}\ c_{23}&&\vdots\\ b^{2}\ c_{13}&b^{2}\ c_{23}&\ddots&\ddots&\\ \vdots&&\ddots&q_{n-1}&b^{2}\ c_{1n-1}\\ b^{2}\ c_{1n}&\cdots&&b^{2}\ c_{1n-1}&q_{n}\\ \end{array}\right)\left(\begin{array}[]{c}S_{1}\\ \vdots\\ \vdots\\ S_{n}\end{array}\right)$$

$$D$$

$$c_{ij}$$

Full text

Structural transitions in nanoscale materials Statistical mechanics of model systems Dynamic critical phenomena

1 Introduction

Neuromorphic circuits are a promising technology to implement at the hardware level the computational power of analog computation inspired by the mammal brain. The type of computation performed by memristors requires a general theoretical understanding of the dynamics, in particular to allow controllability and interpretability of the results. Memristors are becoming the most promising technology for the analogue implementation of artificial intelligence, and their dynamics is known to display memory effects [1, 2, 3, 4], being these very sensible to initial conditions [5, 6, 7]. Memristive circuits are also a new direction of study [8, 9] from a statistical physics standpoint, as these show critical behavior [10, 11] and can be connected to the solution of optimization problems [12, 13, 14]. In its simplest description, a memristor is a 2-ports device behaving as a resistance which changes its value as a function of the flowing current. In this paper we restrict to ideal memristors with zero-crossing in the Voltage-Current diagram [15, 16, 17], though more recently ReRAM (Resistive RAM) devices have further generalized this type of behavior [18].

The analysis presented in this paper could constitute a baseline for modelling and analyzing more complex interacting memristive circuits. We introduce a simple circuit whose asymptotic dynamics we show to be governed by a Lyapunov functional. As we will argue, such a functional can be casted into a spin-like model with long range interactions but with non-linear self-energy. The model we introduce interpolates in fact between a set of non-interacting memristive circuits and a single mesh of memristors. In a recent paper [19], it has been shown however that the interaction strength between memristors for generic circuits is controlled by the Hamming distance on the dual graph of the circuit. That paper however does consider memristors on a mesh, but between meshes. Specifically, we consider the case in which the Hamming distance between each pair of memristors is one, hence representing a fully-connected model in which mean field techniques can be used. The source of our motivation for studying mean-field models of memristors is due to the difficulty of analyzing and understanding the behavior of general memristive circuits. It is the purpose of this paper to show that the model we introduce can be regarded as the analogue of the system studied in the well-known Curie-Weißmodel for spin-spin magnetic interactions. The key difference is that in our case the single interacting element is not the spin of a particle but the internal memory of a memristors, characterized by a different functional form for the pairwise interaction and, most importantly, by a non-Hamiltonian dynamics.

The present paper is structured as follows: in Section 1, we revise the standard model of a simple circuit with one memristor, in Section 2, we define the model of interacting memristors, in which many a memristor are coupled with a central mesh characterized by a given conductance regulating the strengh of the interaction. In Section 3, we analytically and numerically characterize the circuit both in the case of deterministic and random initializations. Finally, in Section 4 we discuss the results and their implications both on real implementations of memristive circuits and on their theoretical modelling.

2 Mean field memristive circuit

We begin by introducing the memristors under study. As first observed in [20], physical memristors slowly relax to a limiting resistance even when a voltage is not applied. This observation implies that there is a competition between a phenomenon of decay and one of reinforcing, which is one of the key mechanisms for the learning ability of biological systems. We consider the simplest dynamical equation which captures such behavior. The time evolution of a simple Ag+ memristor (atomic switch)[21]:

[TABLE]

where $0\leq w(t)\leq 1$ is the internal memory parameter of the memristor, $R(w)=R_{on}(1-w)+R_{off}w$ is the resistance and $I$ and $S$ are the current and applied voltage respectively. Using this parametrization, $R_{on}$ and $R_{off}$ are the limiting resistances for $w=0$ and $w=1$ respectively ($R_{off}>R_{on}>0$). The constant parameters $\alpha$ and $\gamma$ set the timescales for the relaxation and excitation of the memristor respectively111In particular, meanwhile $\alpha$ has the dimension of an inverse time, $\gamma$ has the dimension of time and voltage.. The fixed points $w^{*}$ can be obtained by setting $\frac{d}{dt}w=0$, from which we find the equation

[TABLE]

We immediately observe that this equation is quadratic in $w^{*}$ and that none, one, or two solutions can be obtained depending on the values of the parameters. Physical memristors relax to the state of highest resistance $R_{off}$ at zero voltage, i.e. $\alpha>0$.

The analysis above shows the main characteristics of memristive circuit dynamics: multiple fixed points depending on the value of the external control. With this feature in mind, we now generalize the dynamical model to the case of a circuit composed of $N$ memristors in series to a voltage source, as in Fig. 1: this is a simple modification of the one studied in [22] for machine learning purposes. It was proven in [19] that the interaction strength between memristors decays exponentially with the Hamming distance on the graph. Thus, in order to obtain long range interaction it is necessary to minimize the Hamming distance. The easiest way to obtain this is by arranging each memristor on a single mesh. However, since the memristor internal dynamics depends on the current, when these are arranged on a single mesh the dynamics of the circuit is trivial. The simplest non-trivial dynamics with long range interaction is obtained by inserting a current-divider on a central mesh, as in Fig. 1. Hence, we introduce control resistances in order to have long range interactions between memristors which are tunable. We consider $n$ memristors.

For the derivation of the equations for the circuit of Fig. 1, we first apply the mesh current method [23], in which we assign to each memristive mesh a current $i_{k}$, $k=1,\cdots,n$, and we define the central mesh current $i_{0}$. By construction, in the central mesh there are no memristors but only resistances. Once we have assigned an orientation to the mesh current, the current on each resistance is clear: on the $R_{ref}$ in parallel to the $k$th memristor, we will have a current $i_{0}-i_{k}$, meanwhile in any $R_{int}$ resistance will flow a current $i_{0}$.

On each $k$th memristor the current is given simply by $i_{k}$. The final step is to write the Kirchhoff voltage conservation law for each mesh. Since we have $n+1$ currents and $n+1$ meshes, we have a complete set of equations given by:

[TABLE]

where $S_{k}(t)$ are voltage sources. The system of equations (5) completely determines the circuit dynamics and yields the following memristors’s dynamics, and can be derived via the Sherman-Morrison relation[24]:

[TABLE]

where $M_{ij}=\delta_{ij}\frac{R_{ref}}{R_{ref}+R(w_{i})}\frac{R_{ref}}{R_{ref}+R_{int}}$ and $\mathcal{J}$ is an all-ones matrix (see the Appendix of [25]) and $\mathcal{I}$ the identity matrix. In addition to the Kirchhoff laws above, we need to include the memory dynamics:

[TABLE]

where, more specifically, $f_{i}(\vec{w})=\sum_{j}g_{i}(w_{i},w_{j})$, which reflects the fact that the dynamics is non-linear, pairwise, and all-to-all. Each memristor interacts with every other, through an interaction term $g_{i}(w_{i},w_{j})=\frac{R_{on}}{\gamma}(\mathcal{I}-\frac{1}{n}M\mathcal{J})_{ij}^{-1}\frac{S_{j}}{R_{ref}+R(w_{j})}$ depending both on its own state $w_{i}$ at time $t$ and on the state of memristor $j$. Via a direct computation of the partial derivatives, it is easy to see that the dynamics does not derive from a potential, as $\partial_{w_{i}}f_{j}(\vec{w})\neq\partial_{w_{j}}f_{i}(\vec{w})$. This means that we cannot strictly interpret the behavior of memristive circuits as an Hamiltonian dynamics. In general, the statistical description of such non-Hamiltonian systems cannot restrict itself to the calculation of the equilibrium distribution. In the following subsection, we identify a general approximation to map initial states into asymptotic states, depending on the stability of the fixed points of equations (6).

3 Asymptotic dynamics

In this Section we identify the general fixed point equations of the circuit’s dynamics. The fixed points of this circuit can be obtained from eqns. (5) if we set $\frac{d}{dt}w_{k}=0$. These equations generalize the fixed-point equation obtained for the single memristor. We derive in fact a direct relationship between equilibrium currents and the internal memory $w_{k}$:

[TABLE]

and then find a fixed point equation for the internal memory parameters, which are the solution of the following fixed point equations:

[TABLE]

Considering that for an ideal memristor $R(w_{k})=R_{off}w_{k}+(1-w_{k})R_{on}=R_{on}+(R_{off}-R_{on})w_{k}$, we can rewrite this equation in terms of adimensional quantities only:

[TABLE]

where we defined $\xi=\frac{R_{off}-R_{on}}{R_{on}}$, $\chi=\frac{R_{on}+R_{ref}}{R_{on}}$ and $\rho=\frac{R_{ref}^{2}}{R_{on}(R_{int}+R_{ref})}$. We note that the mean internal memory $\langle w\rangle$ acts as an effective voltage source for the circuit. We see that already for this rather simple circuit, all-to-all interactions can affect the position of the fixed points, as these are modified by a mean field term $\langle w\rangle$ directly proportional to $\rho$.

3.1 Asymptotic distributions for the homogeneous case

Starting from the general dynamics and fixed-point equations, we consider the case of homogenous applied voltage across all memristors, $S_{k}=S$, and consider two different initial conditions. The non-equilibrium dynamics of the system will determine two radically different asymptotic states for these two cases. In the first case, all memristors start from a switched-off condition, $w_{k}(0)=1$ for all $k$. Such uniform condition combined with the homogenous applied voltage results in an exact mean-field dynamics. This implies that all memristors behave alike. This type of dynamics allows us to exactly compute the asymptotic state varying the applied voltage $S$. The critical value for the voltage $S$ can be worked out from eq. (9) and reads:

[TABLE]

Fig. 2 shows the critical line as a function of the parameter $R_{int}/R_{off}$. Above that line, memristors are activated and all converge to the state of lower resistance $w=0$, corresponding to $R(w)=R_{on}$, while below the same line memristors will remain on the state of higher resistance $w=1$, corresponding to $R(w)=R_{off}$. As we will see, this simple analysis shows the existence of a phase transition, that we will study using mean-field techniques.

3.2 Non-interacting case

The case of a random initial condition is radically different as a homogeneous dynamics is not established by the system. In this case, we need to introduce a mean-field approximation of the heterogenouos dynamics of the memristors. This approximation provides an estimate of the asymptotic state of activation. Let us consider first the case in which the memristors are non-interacting, i.e. we assume that $\rho=0$. Using this parametrization, the non-interacting estimate becomes:

[TABLE]

where $\xi$ and $\chi$ are the adimensional parameters defined in eqn. (9), $\frac{S}{\alpha\gamma}$ can be tuned using the external voltage sources, and $\chi$ depends on the interaction between the memristive meshs. We note that $\xi$ is typically positive as $R_{off}\gg R_{on}$, and that $\chi$ cannot be negative for any positive values of the resistances. In the case $R_{ref}=0$, implying $\rho=0$ and $\chi=1$. The solution of equation (11) is:

[TABLE]

We observe that necessarily one root of eqn. (12) falls below zero, and thus only one solution is feasible. In Fig. 3 (top) we plot the numerical solutions obtained for $\alpha<0$ and $\alpha>0$. One important fact that we need to stress is that the dynamics of the circuits greatly depends on the signs of $S$ and $\alpha$. The case in which $S_{i}$’s and $\alpha$ have identical signs is, as the fixed points will not fall in the interval $[0,1]$, and thus the asymptotic state for memristors is binary, either $\{0,\,1\}$. We observe that meanwhile for $\alpha<0$ the asymptotic fixed point is stable, in the case $\alpha>0$ (which is the physical case) it is unstable. A simple calculation of the Jacobian confirms this fact. Fig. 3 (bottom) shows the position of the fixed points as a function of $\alpha$ and $S$ for both values of $\alpha$.

3.3 General case: use of Lyapunov functional

In the case $\alpha<0,S>0$, the asymptotic fixed point can be described as the minimum of a functional. In fact, eqn. (9) can be obtained from $\partial_{w_{i}}H=0$, where $H$ is given by:

[TABLE]

and where $E_{i}(w_{i})=\frac{\xi}{3}w_{i}^{3}+\frac{\chi}{2}w_{i}^{2}$. As we will see, $H(w_{i})$ will be used as an approximate Hamiltonian for the Curie-Wei$\ss$ mean field theory.

It is important to note that the functional in eqn. (13) is a Lyapunov function, as we show in the Appendix of [25], implying that $\frac{d}{dt}H(w_{i})<0$ for $\alpha<0$. Thus, for $\alpha<0$ the asymptotic states $\langle w(t=\infty)\rangle$ are directly connected to fixed points of a functional which can serve as a Hamiltonian. For $\alpha>0$ the Lyapunov functional is given by $L=-H$ in eqn. (13); in which the evolution of the Lyapunov functional is shown in Fig. 4. Although the Hamiltonian of eqn. (13) is reminiscent of an Ising model with long-range interactions but with non-linear self-energy [26, 27], it is worth to mention that here the parameters $w_{i}$ take values in $[0,1]$. To begin with, we anticipate that we will perform all the statistical mechanical calculations at a non-zero temperature, and then take the limit $T\rightarrow 0$ at the end.

The first question we aim to answer is whether we can use the (unphysical) situation $S>0,\alpha<0$ to make any statement regarding the behavior of the system for $S<0,\alpha>0$. As mentioned before, $\alpha>0$ corresponds to a relaxation into an insulating phase, which is the physical case observed in $Ag+$ memristors [20, 21]. We can thus take advantage of a heuristic observation: the fact that location of the unstable fixed determines which fixed point each memristor will reach. For instance, in Fig. 3 (top), if $w_{in}>w^{*}$, since the time derivative is positive, we would expect $w(t=\infty)=1$; if on the other hand $w_{in}<w^{*}$, then we would expect $w(t=\infty)=0$, up to a set of measure zero, $w_{in}=w^{*}$. This also shows that there is a duality between the cases $S<0,\ \alpha>0$, and $S>0,\ \alpha<0$, which is also evident by the fact that since the fixed points depend only on the ratio $S/\alpha$, the position of the fixed point will be unaffected. We can then use the following rule of thumb which connects the probability on the asymptotic states for $\alpha<0$ to the ones of the initial states for $\alpha>0$:

[TABLE]

With this in mind, we use the Lyapunov functional of eqn. (13) as a Hamiltonian in a statistical mechanics setting. If we assume random initial conditions, we can try to predict $\psi=\langle w(t=\infty)\rangle$ using an equilibrium approach in a canonical setting. Given the similarity of the functional to an Ising model, implies a Curie-Weiß approach for the average magnetization $\langle w\rangle$. Using standard mean field theory techniques and after some straightforward calculations (which are provided in the Appendix in [25]), we find the following mean field theory equation at zero temperature [29]:

[TABLE]

Eqn. (15) can be exactly inverted as a function of $\psi=\langle w\rangle$. This gives the same result as the one we would obtain if we substituted $w_{k}\rightarrow\langle w\rangle$ in eqn. (9). This is simply a correction to the physical parameter $\chi$, as in fact we obtain the same effective equation as the non-interacting approximation with the substitution $\chi\rightarrow\chi-\rho$, which is the correction due to the interaction between the memristors. Since the memristor memory is bounded between [math] and $1$, we consider the function $\langle w\rangle=\text{max}\left(0,\text{max}\left(\psi(S),1\right)\right)$.

In Fig. 5 we plot the numerical results on the mean field $\langle w(t=\infty)\rangle$ for $\alpha>0$ and compare these, as a function of $S$, to the non-interacting estimate obtained using eqn. (14). We observe that such approximation fails for larger values of $S$, but yet it provides nonetheless a good estimate for the asymptotic dynamics.

Few comments are in order. First we note that for $\alpha=-1$ the mean field theory calculation exactly reproduces the behavior of $\langle w\rangle$. This approximation suggest a second order phase transition at $S=0$, as in fact one has $\partial_{S}\langle w\rangle$ given by

[TABLE]

which is not a differentiable function.

For $\alpha>0$ we use the approximation of eqns. (14) to calculate the behavior of the system. The validity of this approximation is shown on the right hand side of Fig. 5. We observe that the heuristic assumption of eqn. (14) for $\alpha<0$ is closer to the numerically obtained curve for $S\approx 0$. For larger values of $S$ however, such approximation is less valid. We also observe stronger fluctuations around the mean field theory calculation, which we attribute to the effective instability of the fixed point. Nonetheless, such simple approximation provides a good estimate for the behavior of the internal memory also for larger values of the external applied voltage. The discrepancy is due to the fact that for larger values of $S$, because of the instability aforementioned, there can be trajectories which can invert their course. Since we observe that the real curve lies below the one obtained from the mean field theory, this implies that some memristors whose initial condition lies above the fixed point can invert and reach the asymptotic state $w=0$, rather than $w=1$. We observe that eqn. (14) is valid for $S\approx 0$, this implies that the discontinuity of eqn. (16) is valid also for $\alpha>0$. This seems evident from Fig. 5, where the discontinuity in the first derivative of $\langle w\rangle$ between $\alpha<0$ and $\alpha>0$ is opposite in the sign of the derivative. Also, we note that possibly a divergence would occur if $\chi=\rho$ in eqn. (16), but this does not happen for any positive values of the resistances. From eqn. (16) we can also promptly infer the critical exponent of the transition point $S=0$ is $\frac{1}{2}$.

As a last comment, we note that we can derive the asymptotic susceptibility from the mean field exact solution in the approximation that $S(t)$ is stepwise constant. The result is given by:

[TABLE]

and the derivation provided in the Appendix of [25].

4 Conclusions

Memristors are becoming of interest for their use in on-chip implementations of neural networks [30, 31]. Yet, only a few works focus on their collective dynamics. The results presented in this paper connect the Hamiltonian of interacting spin systems with the Lyapunov function of memristive circuits. This implies a connection between the equilibrium states of a statistical system of spins and the asymptotic states of memristors in a circuit. This statistical mechanics interpretation of memristors, is also analogous to the case of neural networks with temporal delays [28], and provides a solid basis for further generalizations. It has been insofar hard to obtain an analytical control of the dynamics of memristive circuits, in this paper we introduced a mean field theory via a mapping between the equilibrium states into a suitable minima of Lyapunov function. Thus, we have shown, both analytically and numerically, that a first order phase transition occurs for positive $\alpha$, when memristors are initialized to the high-resistance state, and that a second order phase transition occurs when the initial conditions are chosen uniformly at random between the high and low resistance states. Our results on the mapping between non-equilibrium dynamics and equilibrium states could be extended to the more general case of purely memristive circuits. This work gives predictions on the convergence to boundary values of DC-controlled memristors, something which has already being introduced in [32] In future works we will consider the case of noisy memristive dynamics, heterogeneous applied voltages and a generic network topology [9], where complex non-equilibrium glass-like behaviour is expected [33].

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. PB acknowledges support from FET Project DOLFINS nr. 640772 and FET IP Project MULTIPLEX nr. 317532; PB and FC would like to thank the London Institute for Mathematical Sciences where part of this work was done.

5 Single mesh toy model

5.1 Simple derivation

From the mesh circuit equations we have:

[TABLE]

And,

[TABLE]

So that we have:

[TABLE]

That leads to:

[TABLE]

where $\mathcal{J}$ is the all-ones matrix and $M_{ij}=\delta_{ij}\frac{R_{ref}}{R_{ref}+R(w_{i})}\frac{R_{ref}}{R_{ref}+R_{int}}$. And to the dynamics:

[TABLE]

5.2 Alternative derivation and formula

An alternative derivation of the equation above can be obtained via a direct calculation of the inverse. The relation between currents and voltages can be can be written as: {widetext}

[TABLE]

which means we need to invert a matrix of the form:

[TABLE]

which a special case of an arrowhead matrix. The inverse of this matrix is rather complicated, but can be easily obtained by means of a cofactor formula: $(A^{-1})_{ij}=\frac{1}{\text{det}(A)}C_{ji}$ where $C_{ij}=(-1)^{i+j}\text{det}(A^{\tilde{i}\tilde{j}})$ is the determinant of the matrix $A$ where the row $i$ and the column $j$ has been removed. We note that because of the properties of $A$, $C$ is a symmetric matrix. Let us also note that the determinant of matrices of the form as in eqn. (24), $D=\text{det}(M)=\prod_{k=0}^{n}a_{k}-b^{2}\sum_{k=1}^{n}\prod_{j\neq k,j>0}a_{j}$.

Thus:

[TABLE]

We that the cofactor of the matrix in the case in which $i=j>1$ has the same form. We can thus already say that $C_{ij}$ for $i=j>1$ is of the form: $C_{ii}=\text{det}(M^{\tilde{i}\tilde{i}})=\prod_{k=0,k\neq i}^{n}a_{k}-b^{2}\sum_{k=1,k\neq i}^{n}\prod_{j\neq k,j>0}a_{j}$. The special cases $i=1,j>1$ and $j=1,i>1$, take the form $C_{1j}=C_{j1}=-b^{2}\prod_{k\neq 1,k\neq j}a_{k}$. The case $C_{ij}$ with $i\neq j,i,j>1$ has to be calculated on its own. We note that for instance the matrix $C^{12}$ is of the form:

[TABLE]

and thus introduces a zero on the diagonal. For this reason, the determinant of these matrices are of the form $C_{ij}=-b^{2}\prod_{k>0,k\neq i,k\neq j}a_{k}$. We thus have: {widetext}

[TABLE]

where we used the fact that because $S_{0}=0$ and we are also interested in $i_{k\geq 1}$, we can write the equation directly for the side loop. Also, we have implicitly defined:

[TABLE]

We are in particular interested in the inverse of the submatrix which acts only on the memristor currents. This can be easily obtained, and is given by $-\frac{R_{ref}^{2}}{n(R_{int}+R_{ref})}J+(R(\vec{w})+R_{ref})I$, where $J$ is the matrix made of ones and $R(\vec{w})$ is the diagonal matrix with the resistances of each memristor. This implies the following dynamics:

[TABLE]

here $\text{diag}(\vec{x})=\delta_{ij}x_{j}$, and which is the equation we find in the paper.

5.3 Formal mean field dynamics solution and perturbative expansion

In the mean field approximation, eqn. (LABEL:eq:formaleq) can be solved. If we use the Sherman-Morrison formula assuming that all memristors are equal, in such a case the equation becomes:

[TABLE]

with $a=\frac{R_{on}}{\gamma}$, $s=\langle\vec{S}\rangle=\frac{1}{n}\sum_{i=1}^{n}S_{i}$, $c=R_{ref}+R_{on}-\frac{R_{ref}^{2}}{R_{int}+R_{ref}}=R_{on}+\frac{R_{ref}R_{int}}{R_{ref}+R_{int}}$ and $b=R_{off}-R_{on}$, and where we are assuming that $S_{i}=\langle S_{i}\rangle$.

An analytic solution for such an equation can be found in terms of an inverse. Let us define:

[TABLE]

for an arbitrary integration constant $c_{1}$ due to time invariance symmetry. Then, the solution of eqn. (39) is given by the inverse function of $Q(t)$:

[TABLE]

which is not analytical. In order to solve this equation, we use the a perturbative method in $\epsilon=c/b$, assuming that $R_{int}\gg R_{on}\gg 1$. In fact, $\frac{1}{2}<\frac{R_{ref}R_{int}}{R_{ref}+R_{int}}<1$ for positive resistances. In this case, the differential equation becomes:

[TABLE]

We thus search for perturbative solutions $w(t)=w_{0}(t)+\frac{c}{b}\ w_{1}(t)+\cdots$ up to the first order. We have the two differential equations up to the first perturbative order in $\epsilon$, which are:

[TABLE]

The zeroth perturbative order equation has solutions:

[TABLE]

with $z_{0}$ associated with the initial condition, and of which we take the positive sign only. For $t>\frac{1}{2\alpha}$ we have an exponential function of the form:

[TABLE]

The first perturbative order is of the form:

[TABLE]

and has a general solution of the form:

[TABLE]

where we identify $h(t)=\alpha+\frac{as}{bw_{0}^{2}(t)}$ and $g(t)=\frac{a}{bw_{0}^{2}(t)}$. For large times $h(t\gg 2\alpha)\approx\alpha$ and $g(t\gg 2\alpha)\approx\frac{a}{b}\sqrt{\alpha b}^{2}e^{-2\alpha t}=a\alpha e^{-2\alpha t}$. Thus we have:

[TABLE]

We thus obtain the dynamic susceptibility for long times, which is given by:

[TABLE]

which falls off exponentially in time.

5.4 Fixed point structure

Since the interaction matrix is always invertible, we can study the equivalent equation: {widetext}

[TABLE]

which we note we can rewrite the fixed point equation as:

[TABLE]

Fixed points then require that the following two equations are satisfied at the same time:

[TABLE]

for an arbitrary constant $c$ and with $w_{i}\in[0,1]$.

5.5 Lyapunov function

Let us now consider the function

[TABLE]

we have

[TABLE]

and using the equations of motion we obtain

[TABLE]

Let us define

[TABLE]

We have:

[TABLE]

This quantity is positive or negative depending on the sign of $\alpha$ and the eigenvalues of $Q$. If $Q$ is positive definite, then the sign of $\frac{d}{dt}H$ depends only on $\alpha$. We note that $Q$ is the sum of two Hermitean matrices. Thus the minimum eigenvalue of $Q$ satisfies the bound $\lambda_{min}(A+B)\geq\lambda_{min}(A)+\lambda_{min}(B)$ for $A$ and $B$ Hermitean. Thus:

[TABLE]

Since $\frac{1}{n}J_{ij}$ has maximum eigenvalue $1$, we immediately observe that, since $R_{on}$ is positive by construction:

[TABLE]

This implies that $Q$ is positive definite. The function $H$ is thus a decreasing function of the dynamics when $\alpha<0$. Since the function is a weighted negative norm of the derivative of the memristor memories, then it is also zero at the fixed point. For $\alpha>0$, it is sufficient to define $-H$ as a Lyapunov function. In this case, the fixed points become $w=1$ and $w=0$.

6 Mean field theory for random inizialization of homogeneous memristors

We are interested in the low temperature regime of this model. Its partition function can be written as:

[TABLE]

with $H(w)$ from eqn. (13), $\beta=1/T$ and we have implicitly defined the trace:

[TABLE]

We now use the Hubbard-Stratonovich identity, with $m=\frac{\sum_{i}\sigma_{i}}{n}$,

[TABLE]

with $b=\frac{n\rho\beta}{2}$. Let us define $\tilde{E}(w_{i})=E(w_{i})+\frac{\rho}{2n}\sum_{i}w_{i}^{2}$. We write:

[TABLE]

where $Q(\beta,S,\psi)=\text{Tr}_{w}e^{\beta\left((\rho\psi+\frac{S}{\alpha\gamma})w-\tilde{E}(w)\right)}$.If we take the limits $n\rightarrow\infty$ first, for which $\lim_{n\rightarrow\infty}\tilde{E}(w_{i})=E(w_{i})$, which gives

[TABLE]

with $f(\beta)=\text{arg min}_{\psi}\left(\frac{1}{2}\rho\psi^{2}-\frac{1}{\beta}\log\left(Q(\psi)\right)\right)$.

In turn, $f(\beta)$ is given by $\psi$ solution of

[TABLE]

Now we have

[TABLE]

which, in the limit $\beta\rightarrow\infty$ is given by the following mean field equation:

[TABLE]

which is the result presented in the paper. The mean field susceptibility at equilibrium can be calculated from this equation. We have:

[TABLE]

which is the result shown in the main text.