Modeling networks of probabilistic memristors in SPICE

Efficient simulation of probabilistic memristors and their networks requires novel modeling approaches. One major departure from the conventional memristor modeling is based on a master equation for the occupation probabilities of network states [arXiv:2003.11011 (2020)]. In the present article, we show how to implement such master equations in SPICE - a general-purpose circuit simulation program. In the case studies, we simulate the dynamics of ac-driven probabilistic binary and multi-state memristors, and dc-driven networks of probabilistic binary and multi-state memristors. Our SPICE results are in perfect agreement with known analytical solutions. Examples of LTspice codes are included.

However, there is a strong indication that the deterministic description fails when applied at least to certain realizations of resistors with memory [18][19][20].In particular, it was shown experimentally that when a constant voltage is applied to such devices, their state changes in a step-like fashion at random times.In one group of devices, a Poisson distribution of switching times was observed [18][19][20].Furthermore, another group of devices is characterized by a log-normal distribution [21].Several theoretical models were pushed forward to account for the randmoness in the memristor switching [22,23].
The dynamics of networks with discrete-state memristors can be imagined as a sequence of transitions between network states.Recently, we have introduced a master equation approach for the occupation probabilities of the network states [24] that can be used to describe circuits that include binary and multi-state memristors, resistors, voltage and current sources, and possibly some other components 1.In this previous work [24], the solution of the master equation was found analytically for networks of N in-series/in-parallel connected binary memristors driven by a constant voltage source.It has been demonstrated in Ref. [24] that the master equation solution allows to calculate many quantities of interest including various mean switching times, mean current, resistance, etc.
There are two major advantages of the master equation compared to stochastic/Monte Carlo simulations: i) in principle, the master equation can be solved analytically (see Ref. [24] for examples), and ii) using the master equation, many network characteristics can be found in a single calculation without the need for averaging.In the case of symmetries in the circuit, the additional benefit of the master equation is its compactness.This means that a single degree of freedom is required to describe equivalent circuit configurations.
In this article (which is our second work in a series dedicated to probabilistic memristive networks), we introduce a methodology to simulate the probabilistic memristive networks in SPICE.The paper is organized as follows.We start with an overview of the master equation in relation to binary and multi-state probabilistic memristor networks.This is followed by a description of the SPICE implementation scheme supplemented by several examples.In particular, we consider individual probabilistic binary and tri-state memristors driven by ac-voltage, and dc-driven networks thereof.LTspice codes for some of our examples are provided in the Appendix.
The approach presented in this work is relatively general and can be used to model networks combining resistors, probabilistic memristors, constant and time-dependent voltage and current sources.The application of the master equation to probabilistic memristor networks is a paradigm change in the probabilistic memristor modeling, and its SPICE im-1A generalized approach is needed for circuits combining probabilistic memristors and capacitors/inductors. plementation makes it affordable to students and researchers working in the field.

Probabilistic memristors and master equation 2.1 Binary memristors
Binary memristors are characterized by two resistance states, R on and R o f f (with R on < R o f f ) corresponding to the states 1 (on) and 0 (off).The switching between these states is defined by a probabilistic law with voltage-dependent switching rates (inverses of the mean switching times) given by [18][19][20] Here, τ 01 (10) and V 01 (10) are constants and V is the voltage across the device.For a memristor in state 0, the probability to switch to state 1 within small time interval ∆t is γ 0→1 (V)∆t.The probability of swithching from 1 to 0 is defined similarly.
The master equation is written with regard to the occupation probabilities of network states.The network state is defined by a specific combination of the off-and on-states of memristors.For a system containing N binary memristors, there exists 2 N such states.The network evolution consists of a chain of consecutive switchings of memristors (simultaneous switchings can be neglected).On average, such a process is described by the master equation with form where p Θ (t) is the occupation probability of state Θ, Θ m is the network state obtained from Θ by flipping the state of m-th memristor, γ m Θ is the switching rate for m-th memristor in the configuration Θ, and γ m Θ m is defined similarly.The switching rate γ m Θ equals the switching probability (Eqs.(1) or ( 2)) for m-th memristor in the state Θ.
To demonstrate Eq. ( 3), consider two in-series connected identical memristors subjected to a voltage waveform V a (t).There are 4 possible network states that we denote as 00, 01, 10, and 11.In 00, both memristors are in the off-state, in 01, the first is in the off-, while the second is in the on-state, etc. Eq. ( 3) has the form The similarity of memristors is taken into account by relations like γ 1 00 = γ 2 00 , γ 2 01 = γ 1 10 , p 01 (t) = p 10 (t), etc.Therefore, Eqs. ( 5) and ( 6) are the same and the total number of equations that need to be solved reduces by one.In our notation, γ 1 00 describes the switching rate from state 00 with the flipping of the 1-st memristor.The corresponding switching rate is given by Eq. ( 1) with V = V a (t)/2, etc. Importantly, the computation of the switching rate involves the voltage across the switching memristor in the given configuration at the time moment t.

Multi-state memristors
It is assumed that in a K-state memristor the switching between its boundary states (R on and R o f f ) occurs consecutively through K − 2 intermediate resistance states.The master equation (3) preserves its form for multi-state memristor networks, but the network configuration space becomes more complex.Now the indices i, j, k, and so on, in the set Θ = (. . .k ji) denoting the states of the first memristor, the second one, and so on, in the network can have more than two values.Generally, this leads to the exponential growth of the number of network states and, correspondingly, the number of independent equations for occupation probabilities p Θ (t) when N, the number of memristors, increases.Luckily, the number of nonzero switching rates γ, corresponding to the nonzero terms in the right hand side of the master equation (3) for a given network configuration Θ, does not typically grow as fast.
In order to account for potential change in parameter values between resistance states, Eq. ( 1) and Eq. ( 2) are modified to with τ i j(ji) and V i j(ji) being the constant values describing the resistance switching from i( j)-th to j(i)-th memristor state, and i changes from 0 to K − 1.
It is convenient to represent the interdependencies between different occupation probabilities in the master equation using transition schemes.As an example, Figure (1) shows the transition schemes for a single three-state memristor (a) and two such memristors connected into network.An important feature of these schemes is the sequential change in the state of multi-state memristors that approximates the sequential growth of filaments in physical devices.Additionally, we emphasize that the transition schemes in general do not depend on the specific connections of memristors in the network.Such information is contained in the transition rates.In practice some of the transitions may be almost or entirely forbidden.For instance, when a positive voltage is applied to memristor described by Eqs. ( 1) and ( 2), the transition 1 → 0 is forbidden as it occurs at negative voltages.If one neglects low rate and/or forbidder transitions, we obtain the reduced transition scheme, which simplifies the solution of the master equation (3) (see Ref. [24] for some examples).

SPICE modeling approach
Let M be the number of non-equivalent equations for the occupation probabilites (like the set of Eqs. ( 4), (5), and ( 7)).The supremum of M is K N , where K is the number of memristor states, and N is the number of memristors in the network.However, in practical cases M can be much smaller than K N .For instance, if there are N binary (K = 2) identical memristors connected in series, M = K + 1 (see Ref. [24]).
In the SPICE environment, we model each differential equation (such as Eq. ( 4)) by a 1 Farad capacitor charged by a voltage-controlled current source.The occupation probabilities are represented by capacitor voltages.Each source current depends on the voltage across some of the capacitors which forms the right-hand side of the master equation.These circuits are shown in the top rows of SPICE models in Figs. 2, 3, 5 and 6.
To account for the voltage-dependent switching rates (Eqs.( 1)-( 2)), M copies of the network with memristors in non-equivalent combinations of states are utilized.These circuits (shown in the bottom row in Figs. 2, 3, 5 and 6) are connected to the input voltage.The voltages across memristors in these circuits are utilized to calculate the transition rates between the states.
To calculate the mean current, we use a voltagecontrolled current source connected by a resistor to ground to provide a current path.For instance, in the case of in-series connected binary memristors, the current source output is defined by where

Simulation examples 4.1 AC-driven binary memristor
In this simulation, a single binary memristor driven by an ac source is considered as seen in Fig. 2(a).Fig. 2(b) contains the schematic for the SPICE implementation and the corresponding SPICE code can be found in appendix A.1.The memristor has two possible states, R on and R o f f , with resistance values of 1k and 10k Ohms respectively.We used the model parameter values τ 01 = τ 10 = 3 • 10 5 s and V 01 = V 10 = 0.05 V.The ac source,V a (t), has a peak voltage of 1V and is driven at various frequencies.The memristor is initialized in the off-state and will continue switching between the resistance states until the simulation has ended.The current is calculated using B4 and R4 components in Fig. 2(b).The current-voltage curves generated through SPICE simulation can be seen in Fig. 2(c) and they show the frequency behavior typical to deterministic memristive devices [16,17].We verified that Fig. 2(b) SPICE model reproduces some previous results found through Monte Carlo simulations [24].

DC-driven binary memristor network
For this next simulation, we consider a network of binary memristors connected in-series as shown in Fig. 3(a).The network is composed of five memristors driven by a dc source with a voltage of 5V.Fig. 3(b) contains the schematic for the SPICE implementation.Each memristor is identical to one another, meaning the model parameters and the two states are equivalent from memristor to memristor.The memristors have two possible states, R on and R o f f , with resistance values of 1k and 10k Ohms respectively.We used the model parameter values τ 01 = τ 10 = 3•10 5 s and V 01 = V 10 = 0.05 V. T i m e ( m s ) Each memristor starts in the off-state and as time progresses each will switch to the on-state.When a memristor switches to the on-state, the drop in resistance causes an increase in the voltage across the off-state memristors increasing the probability of switching for the off-state memristor.
According to the analytical theory [24], the network mean switching time can be calculated as For the parameters of simulations in Figs. 3 and 4, the above equation gives T 5 = 126 µs.Numerically, the same quantity can be evaluated using the following integral Technically, the integration is done by the components B8 and C7 in Fig. 3, so that the averaged switching time corresponds to the saturation limit of V(Vt) curve in Fig. 4. We emphasize that the analytical and numerical (SPICE) values for T 5 are in full agreement.

Multi-state memristors
The first multi-state simulation considered is a single tri-state memristor driven by an ac source.The ac source has a peak voltage of 1.5 V and is driven at various frequencies.Fig. 5(a) contains the schematic for the SPICE implementation and the corresponding SPICE code can be found in appendix A.2.The memristor now has three possible states, off-, intermediate, and on-state.To account for the added resistance state, a new copy of the memristor network is necessarily added to the SPICE implementation.These states have resistance values of 10k, 3k, and 1k Ohm respectively.The model parameters, τ i j and V i j , are as specified in the SPICE model schematics (Fig 5(a)).The memristor is initialized in the off-state and will continue switching between (c) This next simulation is a network of two tri-state identical memristors driven by a 1.5V dc source shown in Fig. 6(a).The resistance states and model parameters are identical to the memristor used in the previous configuration.Fig. 6(b), the SPICE schematic used for this simulation is shown.The SPICE model is designed according to the transition scheme in Fig. 1(b).The memristors are initialized in the off-state and will switch to the intermediate state before switching to the on-state during the simulation.The evolution of resistance state probabilities for this network is shown in Fig. 6(c) and the mean current as a function of time for this SPICE simulation is shown in Fig. 6(d).The mean current increases in two steps because of the different time scales for the 0 → 1 and 1 → 2 memristor switchings.

Summary
In summary, the use of the master equation in probabilistic circuit modeling [24] offers significant benefits compared to the routine Monte Carlo/stochastic simulations.Many circuit characteristics can be found on average in a single run and the master equation can be, in principle, solved analytically, with several analytical solutions already known [24].In this work, we have shown how to implement the master equation in SPICE.Our examples include simulations of binary and multi-state probabilistic memristors and their circuits subjected to ac-and dc-voltages.We expect that our approach will be useful to a broad range of researchers working in the area of emerging memory devices.

Fig. 2 .Fig. 3 .
Fig. 2. Ac-driven probabilistic binary memristor: (a) simulated circuit, (b) schematics of SPICE model, and (c) example of current-voltage curves found with SPICE simulations.The listing of SPICE model is given in Table A.1.
n c u r r e n t ( m A )

Fig. 4 .
Fig. 4. Current as a function of time (black solid line), and calculation of the network switching time (red dashed line) in the dc-driven network of five probabilistic binary memristors.

Fig. 5 .
Fig. 5. Ac-driven probabilistic three-state memristor: (a) schematics of SPICE model, and (c) example of current-voltage curves found with SPICE simulations.The listing of SPICE model is given in Table A.2.The simulated circuit is the same as in Fig. 2(a) with the difference of different memristor type used.

Fig. 6 .
Fig. 6.Dc-driven network of two three-state memristors: (a) simulated circuit, (b) schematics of SPICE model, (c) time-evolution of occupation probabilities, and (d) current as a function of time.
the number of states with the same number of memristors in the on-state is taken into account by the binomial