An analytical approach to engineer multistability in the oscillatory response of a pulse-driven ReRAM

A nonlinear system, exhibiting a unique asymptotic behaviour, while being continuously subject to a stimulus from a certain class, is said to suffer from fading memory. This interesting phenomenon was first uncovered in a non-volatile tantalum oxide-based memristor from Hewlett Packard Labs back in 2016 out of a deep numerical investigation of a predictive mathematical description, known as the Strachan model, later corroborated by experimental validation. It was then found out that fading memory is ubiquitous in non-volatile resistance switching memories. A nonlinear system may however also exhibit a local form of fading memory, in case, under an excitation from a given family, it may approach one of a number of distinct attractors, depending upon the initial condition. A recent bifurcation study of the Strachan model revealed how, under specific train stimuli, composed of two square pulses of opposite polarity per cycle, the simplest form of local fading memory affects the transient dynamics of the aforementioned Resistive Random Access Memory cell, which, would asymptotically act as a bistable oscillator. In this manuscript we propose an analytical methodology, based on the application of analysis tools from Nonlinear System Theory to the Strachan model, to craft the properties of a generalised pulse train stimulus in such a way to induce the emergence of complex local fading memory effects in the nano-device, which would consequently display an interesting tuneable multistable oscillatory response, around desired resistance states. The last part of the manuscript discusses a case study, shedding light on a potential application of the local history erase effects, induced in the device via pulse train stimulation, for compensating the unwanted yet unavoidable drifts in its resistance state under power off conditions.

in such a way to induce monostability or different forms of multistability (see section "Application of the theory to endow the ReRAM cell with three, four, or five oscillatory behaviours") in the device oscillatory behaviour upon request.The capability of the ReRAM cell to act as a monostable or multistable oscillator under suitable periodic pulse train stimulation may be leveraged to develop novel forms of data detection and computation in memory for artificial intelligence applications in the years to come.For example, as revealed in section "Compensating for the drift in the resistance of crosspoint devices under power off conditions", the regular application of a suitable generalised pulse train voltage signal across each crosspoint nanodevice may allow to correct an unwanted drift in its resistance state under power off conditions.Finally, the conclusions, summarising the most significant results of this research study, are drafted in section "Conclusions".

Memristor Model
The Strachan model 4 falls in the class of first-order extended voltage-controlled memristors 13 , defined via the DAE set 14   where the ODE (1), referred to as state equation, dictates the rate of change of the memory state x of the oneport, as an input voltage signal v is let fall between its terminals.The algebraic constraint (2), known as state-and input-dependent Ohm law, defines how state and voltage affect the flow of the output current signal i through the device stack.In (1) ((2)) g(x, v) (G(x, v)) represents the state evolution (memductance) function.Let us assume x to be constrained to lie at all times within a closed set D [x L , x U ] .In the Strachan model the state evolution function reads as where step(•) is the Heaviside function, while the SET g SET (x, v) and RESET g RESET (x, v) components of g(x, v), referred to as SET and RESET state evolution functions, and governing the evolution of the device memory state under positive and negative input voltages, respectively, are in turn defined as in which p = i • v denotes the power dissipated in the memristor, as a voltage signal is applied between its ter- minals.The formula for the memductance function in the Strachan model assumes the form For any given voltage v, the higher is the memory state x, and the larger is the memductance.In the remainder of this paper, the ReRAM cell model, consisting of the ODE (1), with state evolution function (3), where the SET and RESET components are respectively expressed by equations ( 4) and (5), and of the algebraic relation (2), with memductance function (6), is referred for simplicity as Strachan DAE set.Table 1 reports the values assigned to the parameters in Eqs.(4), (5), and (6) so as to allow the resulting model to reproduce experimental data, extracted from a Ta 2 O 5−x physical sample, to within a preliminarily specified degree of accuracy 3 .

Theoretical tools
This section introduces the theoretical concepts applied in the research study discussed later on.
(1) ẋ =g(x, v), and • exp p σ p , and (5) Table 1.Strachan model parameters fitted to a Ta 2 O 5−x physical sample.The lower x L and upper x U bounds in the state existence domain D are respectively equal to 0 and 1.

The time average state dynamic route technique
When a voltage signal v S falls across the ReRAM cell, as illustrated in Fig. 1a, the time average x of its memory state x, referred to as time average state for short, evolves with time according to the formula Taking the time derivative of both sides of Eq. (7) gives where the last step follows from the integration of the state equation (1) for v = v S .In principle Eq. ( 8) could be employed to explore the response of the device to any periodic stimulus, but then recurring to some numerical integration method would be necessary for determining its solutions.However, as explained shortly, for stimuli composed of rectangular pulses of suitable widths and heights, the study of the behaviour of the periodically-forced device may be considerably simplified, and some analytical developments, following from a simplification of Eq. ( 8), are possible.Let us first consider the excitation scenario analysed in the bifurcation study from Pershin and Slipko 8 .With reference to Fig. 1b, the train voltage stimulus v S , applied across the memristor, features here a first τ + -long SET pulse of positive polarity and amplitude V + and a second τ − -long RESET pulse of negative polarity and amplitude V − over each cycle of length T = τ + + τ − .Equation ( 8) reduces then to Assuming the positive (negative) SET (RESET) pulse induces a relatively small increase (decrease) in the device memory state over the first (second) τ + ( τ − )-long part of each cycle, it is possible to substitute the state x with its time average x in each integrand without introducing a large error in the resulting approximation.Equation (9) boils then down to where (7)  RESET pulse per cycle.The RESET pulse of height V − and width τ − follows the series of SET pulses.The i th SET pulse is V +,i high and τ +,i wide, with i ∈ {1, . . ., P} .The ordering of the positive pulses from the lowest to the highest in each input cycle follows the convention adopted in the systematic methodology to engineer multistability in the steady-state oscillatory response of the ReRAM cell to a generalised train stimulus (refer to section "A systematic methodology to craft the pulse stimulus for enabling the ReRAM cell to support multiple oscillations around prescribed resistance levels".However, this has no effect on the simulations.In fact, to facilitate their convergence, in the numerical investigations, discussed in section "Conclusions", the SET pulses were listed from the most narrow to the most wide before being applied in this order one after the other across the device.10) (refer to Fig. 2c), is the so-called time average state dynamic route (TA-SDR) resulting from the earlier arbitrarily specified pulse train stimulation of the nanodevice 16 .A state dynamic route (SDR), namely the ẋ versus x locus, derivable from the state Eq.(1) for a given DC value V assigned to the voltage v, governs the time evolution of the memory state of a first-order memristor under the specified bias stimulus.In this regard, it is worth to observe that a number of research studies have recently www.nature.com/scientificreports/ reported laboratory measurements of SDRs acquired from memristive nanodevices, enabling to establish an important communication channel between theoreticians and experimenters.The interested readers are invited to consult works from Messaris 17 , Maldonado 18 , and Marrone 19 for the details.A TA-SDR can then be interpreted as an extension of the SDR, enabling the investigation of the response of the same device to a particular AC periodic square pulse train.On condition that the choice of the pulse train parameters does not jeopardise the accuracy of the approximation inherent in Eq. ( 10), the analysis of this new graphic tool enables to determine number, mean values, and stability properties of all the admissible asymptotic oscillations in the memory state of the periodically-forced device.Moreover the predictive capability of the TA-SDR technique may be verified by means of another more rigorous system-theoretic methodology, described shortly in section "The state change per cycle map tool".

The state change per cycle map tool
The State Change Per Cycle Map (SCPCM) analysis tool 20 is inspired from the Poincaré map technique 12 , a powerful method from Nonlinear Dynamics Theory, which facilitates the study of a nth-order non-autonomous periodically-forced continuous-time system, equivalent to a (n + 1)th-order autonomous continuous-time system, where the time assumes the role of a state variable, through the analysis of a simpler n-dimensional discrete-time one.For the Strachan model, featuring order n = 1 , and continuously driven by an input signal v, forced to follow a generic periodic voltage stimulus v S , the Poincaré map assumes the one-dimensional form where k ∈ N >0 , and x k stands for the sample x(k • T) of the solution of the ODE (1), with g(x, v) expressed by the formula (3), at the end of the kth input cycle.For k = 1 the map reduces to x 1 = P (x 0 ) , where x 0 denotes the ODE initial condition x(0), and x 1 is the state sample x(T) at the end of the first input cycle.The Poincaré map accurately provides the sequence of values x 0 , x 1 , . . ., also referred to as return points, extracted from the time series of the ODE solution at regular T-long time intervals from the initial instant t = 0 of the simulation.The one-dimensional discrete-time system is said to admit a fixed point x * if it maps such point into itself, which is mathematically formulated via the equality x * = P (x * ) .A fixed point of the map corresponds to a steady-state oscillatory solution for the original non-autonomous continuous-time system.However, the periodic attractor of the non-autonomous ODE (1) is asymptotically stable if and only if the fixed point of the map is also asymptotically stable, which implies the inequality |P ′ (x k )| x k =x * < 1 to hold true. Figure 3a sketches qualitatively how the graph of the map may look like for an exemplary case study, where, similarly as assumed in Fig. 2  , of which the outer ones (inner one) are stable (is unstable).A few coloured zig-zag trajectories, known as cob-web plots 12,20 in Nonlinear Dynamics Theory, are also displayed to show the discrete-time evolution of the map from distinct initial conditions toward one of the two LAS fixed points.In our study a map of this kind can be extracted from the Strachan DAE set, when the input voltage v is enforced to follow a given periodic voltage stimulus v S , e.g. in the form of a rectangular pulse train, by recording samples of the memristor state x at regular T-long time intervals from the beginning of each of a large ensemble of simulations, differing in the initial conditions, and then plotting for each of the resulting time series the kth sample x k = x(k • T) versus the (k − 1) th one x k−1 = x((k − 1) • T) , with k ∈ N >0 .For k = 1 the SCPCM reduces to � 1;0 = x 1 − x 0 = P (x 0 ) − x 0 , providing the change in the memory state over the first input cycle. (b) x k;k−1 = x k − x k−1 versus x k−1 locus, illustrating the SCPCM of the ReRAM cell subject to the periodic stimulus, which induces a state motion resulting in the Poincaré map shown in plot (a).www.nature.com/scientificreports/ i.e. within the kth input cycle, as a function of its value x k−1 at time t = (k − 1) • T , i.e. either at the end of the (k − 1) th input cycle, if k > 1 , or at the beginning of the simulation, if k = 1 .The SCPCM admits a graphical visualisation on the x k;k−1 versus x k−1 plane, as sketched qualitatively in Fig. 3b, corresponding to the P (x x−1 ) versus x k−1 locus in plot (a) of the same figure.For any initial condition x 0 from a set of values uniformly distributed across the state existence domain D , the net change �x(k; k − 1) in the state x over the time interval [(k − 1) • T, k • T] may be marked on this plane at the abscissa corresponding to the state value x k−1 at t = (k − 1) • T for each k ∈ N >0 .A suitable interpolation method can then be employed to derive the curve, which best fits the sequences of return points collected for all the selected initial conditions.Arrows, pointing to the east (west), are then superimposed along the graph of a SCPCM in the upper (lower) half of the x k;k−1 versus x k−1 plane to indicate a progressive step-wise increase (decrease) in the discrete-time evolution of the Poincaré return point when x k;k−1 is positive (negative).For each k value in N >0 , the k th return point x k of the Poincaré map for a given initial condition x 0 may be obtained by adding the abscissa x k−1 , representing either the (k − 1) th return point, if k > 1 or the initial condition, if k = 1 , to the ordinate x k;k−1 of the point of intersec- tion between the graph of the SCPCM and the vertical line passing through the point(x k−1 , 0) .A fixed point x * of the Poincaré map corresponds to the state value, at which the SCPCM crosses the x k−1 axis, as x k;k−1 = 0 therein.The stability of a fixed point of the map may be inferred by monitoring the direction of the arrows in its neighbourhood.Arrows, pointing toward (away from) a fixed point on both its left and right sides, provide clear evidence for its asymptotic stability (instability).Alternatively, the same information can be retrieved by inspecting the slope of the graph of the SCPCM at the fixed point.In fact, the fixed point is asymptotically stable (unstable) if and only if the slope of the x k;k−1 versus x k−1 locus is negative (positive) at its location.
Remark 1 The vector field of the original non-autonomous ODE (1) maps a given state value into some other one over a T-long time span, irrespective of the number of input cycles, elapsed since the beginning of the simulation.Therefore, a more efficient strategy to compute a SCPCM, in comparison to the method described earlier, envisages to test the response of the ReRAM cell to a predefined periodic stimulus across a T-long time span only, for each initial condition x 0 from a set of values uniformly distributed across the state existence domain D .Specifically, in each iteration step, the state value x 1 x(T) at the end of a T-long simulation, and the state change x 1;0 , relative to the initial condition, would be recorded, allowing to identify a particular point on the plane, spanned by x 0 and x 1;0 on the horizontal and vertical axes, respectively.Interpolating the data through some best-fit curve, and renaming the label on the horizontal (vertical) axis as x k−1 ( x k;k−1 ) finally results in the graph of the SCPCM of interest.
All in all, the SCPCM technique enables to explore the response of a first-order nonlinear dynamic system to any periodic stimulus.It thus extends the applicability scope of the TA-SDR tool, which is employable solely in those case studies, where a periodic train of rectangular pulses stimulates a system of this kind.Furthermore, following the steps, described in Remark 1, it should be possible to acquire a SCPCM experimentally, provided access to the device state were possible.On the other hand, the measurement of a TA-SDR seems to pose harder challenges.Moreover, as elucidated in this section, no approximation is involved in the derivation of a SCPCM, which, as a result, may be used to verify the predictions of the TA-SDR investigation technique.Despite its weak points, however, the latter method allows the derivation of an analytical approach to engineer multi-stability in the oscillatory response of the nano-device to a generalised periodic pulse train voltage stimulus, as described in section "An analytical methodology to operate the pulse-driven ReRAM cell as a multimodal device with initial condition-dependent oscillatory behaviour".Moreover, upon availability of a reliable model for the ReRAM cell, and for any given rectangular pulse train stimulus, the computation of the SCPCM takes a much longer time than the determination of the respective TA-SDR.In fact, while the latter task simply requires to plot the right hand side of the TA-SE (10), adapted to the excitation signal of interest, against the time average state, the first one requires the numerical integration of the state equation over one input cycle for an adequate number of initial conditions, as explained in Remark 1.

Insights into the model
Before introducing the analytical framework, allowing to induce mono-or multi-stability in the oscillatory response of the ReRAM cell to a generalised pulse train voltage stimulus, this section discusses numerical investigations, which shed light into the properties of the Strachan model equations as well as into its response to a square wave excitation signal from the class illustrated in Fig. 1b.

Dynamic route map
The Dynamic Route Map (DRM) of the first-order ReRAM cell under focus is a family of SDRs, each of which corresponds to the plot of the state evolution function (3) against the state for a particular DC value V assigned to the voltage v.When V is negative (positive), the resulting g(x, V) versus x locus is referred to as a RESET SDR (SET SDR).A number of RESET (SET) SDRs, obtained by sweeping |V| in 0.2 V-long steps from 0.2 V to 1 V , are shown in plots (a), (c), (e), (g), and (i) ((b), (d), (f), (h), and (l)) of Fig. 4. As may be inferred by inspecting the graphs on the left (right) column of this figure, the choice of the negative (positive) DC value V has no significant (a strong) impact on the shape of the resulting RESET (SET) ẋ versus x locus.We may thus conclude that, toward the development of a strategy to massage the SET ẋ| SET and RESET ẋ| RESET components of the TA-SE (10) in such a way to enable a desired number of intersections between their graphs, the fine control of the position of the gaussian bell-shaped g(x, V) versus x locus across the horizontal axis through smooth changes in the positive DC voltage V is worth of exploitation.τ − , depending upon the selection of the RESET V − and SET V + pulse heights.This is clearly illustrated in Fig. 5a, visualising through a three-dimensional surface each admissible equilibrium xeq = xeq (V + , V − ) , which the TA-SE (10) admits for r = 1 , endowing the pulse train with a 50% duty cycle, when V − and V + are in turn chosen as the abscissa and ordinate of any point of the coloured map in plot (b) of the same figure.The dark blue domain from plot (a) contains the only globally asymptotically stable (GAS) equilibrium xeq , which the TA-SE features, upon selecting (V − , V + ) anywhere within the green region from plot (b).On the other hand, the bottom and top violet domains (the cyan domain) from plot (a) include (includes) the leftmost and rightmost LAS equilibria xeq,1 and xeq,3 (the unstable equilibrium xeq,2 ) of the TA-SE corresponding to any choice for the input parameter pair within the red domain from plot (b).For the sake of completeness, the white area in the coloured map of plot (b) contains input parameter pairs, whereby there exists no state value, at which ẋ| SET = − ẋ| RESET .In a scenario of this kind, if ẋ is found to be strictly negative (strictly positive), the memory state of the periodically-forced ReRAM cell shall progressively decrease (increase) toward the lower (upper) bound x L ( x U ) in its existence domain D .As the operation of the device around its fully-RESET or fully-SET state is not recommendable, the selection of input parameter pairs, belonging to the white region in the map of plot (b), shall be avoided in the analysis of exemplary excitation case studies to follow.

Monostability
Taking V − and V + in turn as the abscissa and ordinate of the point (−0.4V, +0.46 V) , indicated as a black cross marker, and belonging to the green region in the map of Fig. 5b, as may be inferred from Fig. 6a, showing the | ẋ| SET | and | ẋ| RESET | versus x loci for the earlier specified (V − , V + ) pair, the TA-SDR analysis predicts a monostable oscillatory behaviour for the periodically-forced ReRAM cell.Its memory state x is expected to experience a steady-state oscillation around the TA-SE equilibrium xeq = 0.308 .With reference to plot (a) of Fig. 5, the left vertical black dashed line crosses the blue domain of the three-dimensional surface in a single point, specifically (V − , V + , xeq ) = (−0.4V, +0.46 V, 0.308) , as indicated through a green filled circle.Choosing sufficiently small values for the RESET τ − and SET τ + pulse widths is instrumental to prevent the error in the approximation inherent in Eq. ( 10) to jeopardise the accuracy of its predictions.Fixing both τ − and τ + to 1µs , the SCPCM of the ReRAM cell, shown in Fig. 6(b), confirms the conclusions drawn via TA-SDR analysis from plot (a) of the same figure.Figure 6c shows the progressive approach of the solution x to the state Eq.(1) toward the only possible asymptotic periodic waveform, revolving approximately around xeq = 0.308 , from either of two initial conditions, lying one well below the minimum and the other well above the maximum of the steady-state oscillation.Plot (d) of the same figure visualises both the periodic pulse train voltage stimulus v S (in blue) and the steady-state oscillation x ss in the memristor state (in green) together with its mean value xss , its prediction,  10), associated to a train voltage stimulus, featuring two pulses of opposite polarity per cycle, may possibly admit, when the SET τ + and RESET τ − pulse widths are identical, as a function of the SET V + and RESET V − pulse heights, swept across the ranges [−2, 0] V and [0, 1.2] V , respectively.The dark blue surface includes all the GAS equilibria of the TA-SE in the monostable oscillatory operating mode of the ReRAM cell.The cyan (magenta) surface contains all the unstable (all the LAS) equilibria of the TA-SE in the bistable oscillatory operating mode of the ReRAM cell.(b) Projection of the surface from (a) onto the V + versus V − plane.Choosing the pulses' heights of the pulse train voltage stimulus, featuring a 50% duty cycle, according to the coordinates of any point in the green (red) region, the TA-SE features a single GAS equilibrium (two LAS equilibria) for r = 1 .The black cross marker (black plus sign) identifies the input parameter pair (V − , V + ) , inducing the particular monostable (bistable) oscillatory response, illustrated in Fig. 6 (Fig. 7), in the nanodevice.
Vol:.( 1234567890 www.nature.com/scientificreports/namely the TA-SE equilibrium xeq , as well as the map fixed point x * , corresponding to the minimum value it assumes over each cycle.

Bistability
Setting the RESET V − and SET V + pulse heights to −0.6V and +0.54V , respectively, which identifies a point, indicated through a black plus sign, and lying within the red region in the map of Fig. 5b, the TA-SDR analysis predicts the coexistence of two LAS steady-state oscillatory solutions for the memory state x of the periodicallyexcited ReRAM cell, as may be inferred from Fig. 7a, revealing the existence of a triplet of crossings between the loci of the moduli of the SET and RESET TA-SE components.The abscissa of each of the two outer intersections xeq,1 = 0.106 and xeq,3 = 0.370 (of the inner intersection xeq,2 = 0.237 ) denotes a LAS (an unstable) equilib- rium for Eq.(10).The right vertical black dashed line in Fig. 5a intersects the bottom and top violet domains in the green-filled points   cell under the application of a two-pulse-per-cycle pulse train voltage stimulus v S , when its SET V + and RESET V − pulse heights are in turn set to +0.46 V and −0.4 V , and for r = 1 , irrespective of the choice of its SET τ + and RESET τ − pulse widths.Note that scaling the widths of the 2 pulses in the train per cycle by the same factor does not affect the TA-SDR prediction.The only GAS equilibrium xeq of the TA-SE lies at 0.308, which is the abscissa of the black-filled circle.A marker, indicating the zero of the RESET component at x = 0 , is omitted from the graph, so as to avoid clutter.(b) SCPCM of the ReRAM cell subject to a particular pulse train voltage stimulus v S , belonging to the class considered in (a), and characterised by parameters ) (refer to the blue signal of period T = τ + + τ − = 2 µs in plot (d)).The Poincaré map, from which it is extracted, features a GAS fixed point x * (see the black-filled circle).Differently from what is the case for the TA-SDR, scaling the widths of the 2 pulses in the train per cycle by the same factor may affect the SCPCM.(c) Brown (Green) trace: progressive approach of the solution x to the Strachan DAE set, when v is forced to follow the particular excitation voltage signal v S , employed for the derivation of the SCPCM, from the initial condition x 0 = x 0,1 = 0.15 ( x 0 = x 0,2 = 0.85 ) toward a unique steady-state oscillation.(d) Green trace: steady-state time series x ss of the memristor state x, as extracted from the solution featuring the same colour in plot (c).Horizontal lines mark the locations of the map fixed point x * , of the TA-SE equilibrium xeq , and of the time average xss of the steady-state time series.As the RESET pulse follows the SET pulse over each cycle of the input train, x ss attains its minimum value at the end of any period.Therefore x * directly reveals the minimum of x ss across one input cycle.
).In each of the two cases the choice of the initial condition ensures that no transients appear in the device response.The time average x1 ( x3 ) of the solution x 1 ( x 3 ), as well as the corresponding LAS TA-SE equilibrium xeq,1 ( xeq,3 ) and LAS map fixed point x * 1 ( x * 3 ) are also marked in plot (e, f).
Vol:.( 1234567890) www.nature.com/scientificreports/respectively, and the cyan domain in the red-filled point (V − , V + , xeq ) = (−0.6V, +0.54 V, 0.237) .Setting the RESET τ − and SET τ + pulse widths to a relatively small value, specifically 40 ps , as shown in plot (b) of Fig. 7, visualising the time waveform of the resulting train voltage stimulus v S , allows to limit the change in the memory state over each cycle, which endows the TA-SDR graphic tool with predictive capability.In fact, as may be inferred from plot (c) of the same figure, the SCPCM of the ReRAM cell, subject to the stimulus from plot (b), validates the conclusions drawn through the analysis of the TA-SE (10).The cyan (violet) trace in Fig. 7d depicts the transient behaviour of a solution to the ODE (1), as it approaches the LAS oscillatory waveform revolving approximately around the leftmost (rightmost) TA-SE equilibrium xeq,1 ( xeq,3 ).Due to the slow/fast dynamical effects, emerging in the nanodevice, it takes a rather long (rather short) time for the first (second) solution to attain the steady state.However, an ad hoc choice of the initial condition may allow to retrieve the asymptotic behaviour of the state of the periodically-driven ReRAM cell without the need to wait for transients to vanish.In fact, plot (e) ((f)) of the same figure illustrates the transient-free solution x 1 ( x 3 ) to the ODE (1), initiated from the leftmost (righmost) map fixed point x * 1 ( x * 3 ), corresponding to the minimum value the state assumes over each cycle, together with its mean value x1 ( x3 ), and the respective approximation xeq,1 ( xeq,3 ).

On the crossings between one scaled SET SDR and one scaled RESET SDR
In general, a single gaussian bell-shaped SET SDR may cross a single RESET SDR nowhere, which is of no practical interest, as elucidated in section "Numerical investigation of the ReRAM response to the basic pulse train stimulus", in one point, endowing the resulting TA-SE with a GAS equilibrium xeq , or in three locations, specifically xeq,1 , xeq,2 , and xeq,3 , the outer of which denote LAS equilibria for the corresponding TA-SE.For a fixed choice of the negative pulse height V − , this depends upon the amplitude V + of the positive pulse as well as upon the SET-to-RESET pulse width ratio r, as may be inferred from the coloured map, shown in Fig. 8a, which was derived by means of a numerical procedure for V − = −0.5 V , and depicts through a white, a green, and a red hue the regions of the r versus V + plane, where, according to the TA-SDR analysis, the ReRAM cell is expected to admit no, a monostable, and a bistable oscillatory behaviour at steady state, respectively.Fig. 8b, c) illustrates the loci of the SET and RESET TA-SE components for a choice of the input parameter pair (V + , r + ) , specifically (+0.50 V, 1 × 10 8 ) ((+0.75 V, 1 × 10 −30 )) (see the black cross marker (black plus sign) within the green (red) domain in plot (a) of the same figure), which is expected to trigger a monostable (bistable) oscillatory response in the ReRAM cell, when V − is fixed to −0.5 V.

An analytical methodology to operate the pulse-driven ReRAM cell as a multimodal device with initial condition-dependent oscillatory behaviour
Despite an in-depth numerical investigation may allow to explore the response of the ReRAM cell to a rectangular pulse train stimulus, the availability of an analytical strategy to craft the excitation signal so as to endow the memory state of the ReRAM cell with a prescribed number of steady-state oscillatory solutions, revolving around predefined levels, would be of greater interest for circuit designers.In order to address this point, this section first presents a thorough analytical investigation of the Strachan model, and then employs its findings to propose a systematic approach to engineer multistability in the oscillatory response of the ReRAM cell to a generalised pulse train voltage stimulus from the class defined in the section to follow.

Adaptation of the TA-SE to a generalised pulse train stimulus
In fact, here the periodic voltage source v S in the test circuit of Fig. 1a is assumed to emit a generalised pulse train from the class illustrated in Fig. 1c.In each cycle the generalised train is composed of a tunable number P ∈ N >0 of positive SET pulses followed by a single RESET pulse.Let the ith SET pulse feature a height V +,i and a width τ +,i , with i ∈ {1, 2, . . ., P} .The height and width of the RESET pulse are indicated as V − and τ − , respectively.The input period is thus computable as T = τ +,1 + τ +,2 + . . .+ τ +,P + τ − .Under this hypothesis, Eq. ( 8) may be expanded as Let us further assume each pulse in the train of Fig. 1c to induce a negligible change in the device memory state.This allows to approximate the state x, appearing in each integrand function from Eq. ( 13) with its time average state x , allowing to derive the TA-SE of the ReRAM cell, subject to the generalised train voltage stimulus.In its approximate formula, still provided by Eq. ( 10), the RESET component keeps the expression reported in (12), while the SET component reads as with τ+,i τ +,i /T.

Extraction of key geometrical features from a gaussian bell-shaped SET state evolution function
The proposed strategy to endow multistability in the oscillatory response of the ReRAM cell to a generalised pulse train envisages an ad hoc choice for the P + 1 stimulus parameters (V +,1 , τ +,1 , V +,2 , τ +,2 , . . ., V +,P , τ +,P , V − , τ − ) , for P ∈ N >0 , so as to shape the locus of the SET TA-SE component ẋSET in such a way to let it intersect the locus of the RESET TA-SE component ẋSET , while keeping above (below) it to the left (right) of the crossing, in as many locations as specified in the design requirements.In fact, it is fundamental to devise an ad hoc linear combination between scaled gaussian bell-shaped SET state evolution functions for massaging the SET TA-SE component according to the design specifications.The derivation of a few key geometrical properties of a generic gaussian bell-shaped SET SDR, i.e. the g SET (x, v) versus x locus, resulting from assigning an arbitrary positive DC value V + to the voltage v (recall the graphs along the right column of Fig. 4), is instrumental to accomplish this goal.

State value at the peak of a SET SDR
This section derives an exact closed-form expression for the state value, at which a generic gaussian bell-shaped SET SDR attains its peak level.Employing Eqs.(1), (3), and (4), the rate of change ẋ of the memory state x under the application of a positive bias voltage V + across the device may be cast as where ( 13) with G(x, V + ) expressing the dependence of the device conductance upon its state under the prescribed posi- tive DC stimulus, according to Eq. ( 6).The abscissa of the maximum of the ẋ versus x locus for v = V + may be analytically computed by employing Eq. ( 15), and finding the state value, at which ∂g SET (x, V + )/∂x vanishes.As the exponential function on the right hand side of Eq. ( 15) is a monotonically increasing function of its argument, it is sufficient to find the state value, at which ∂α(x, V + )/∂x vanishes.After some algebraic calculation, the formula for x max (V + ) is found to read as where Figure 9a shows the x max versus V + locus, extracted from the formula (17) (red trace) together with its numerical approximation (blue trace).

Positive DC voltage for programming the abscissa of the peak of a SET SDR
This section derives an approximate analytical formula for the positive bias voltage V + to be assigned to the volt- age v for the respective ẋ = g SET (x, v) versus x locus to feature the peak level at a preliminarily specified state value x max .The function in Eq. ( 17) may not be inverted analytically, which explains the reason for the search of a suitable approximate formula.The solid blue trace in Fig. 9b shows the dependence of the function γ (V + ) upon V + , as it descends from the exact closed-form expression (18).Let us approximate the expression for γ (V + ) with a quadratic polynomial of the form where a 0 (V +,1 , V +,2 ) and a 1 (V +,1 , V +,2 ) are functions of two voltage parameters, specifically V +,1 and V +,2 , allowed to range between 0 V and 1 V , and shortly subject to an optimisation procedure.Importantly, a 0 and a 1 are strictly positive-and negative-valued, respectively, since, as may be evinced by inspecting the blue trace in Fig. 9b, the original function γ (V + ) features a positive polarity for V + = 0 V , and admits a downward concav- ity, respectively.Replacing γ (V + ) with γ (V + , V +,1 , V +,2 ) into the formula (17) for x max delivers an approximate analytical formula for the abscissa of the peak of the gaussian bell-shaped SET SDR, reading as Inserting now the second-order polynomial (19) in place for γ (V + , V +,1 , V +,2 ) into this equation yields the biquadratic equation (17) The exact analytical solution, descending from the formula (17), is illustrated in red.The numerical solution, depicted in blue, saturates abruptly to the unitary value at the first positive DC voltage V + , specifically 0.957 V , where x max exceeds the upper bound x U of the state existence domain D , keeping unchanged for any larger V + value.(b) Blue trace: Graph of γ as a function of V + , according to the exact analytical formula (18).Red trace: approximation of the γ versus V + locus via the analytical function γ (V + , V +,1 , V +,2 ) from Eq. ( 19) for +,2 ) = (0.662, 0.923)V .(c) Positive value V + to be assigned to the DC voltage V in order for the abscissa of the peak of the resulting SET SDR to lie at a pre-specified location x max .The blue curve shows the V + versus x max locus determined numerically from the blue-coloured numerical solution in (a) by exchanging the data series reported along horizontal and vertical axes.At x max = 1 the blue trace abruptly turns into a vertical segment stretching from V + = 0.957 V to V + = 1 V .The red curve is the Ṽ+ versus x max locus, extracted from the analytical formula (22), proposed to approximate the inverse of the function (17), for V +,1 = V (opt) +,1 , and V +,2 = V (opt) +,2 .(d) Blue trace: graphical illustration of the exact analytical formula (17) for x max .Red trace: xmax versus V + locus, obtained from the approximate closed-form expression (20) for . (e) Peak value g SET, max (V + ) of a SET SDR as a function of the positive DC voltage V + across the ReRAM cell.The red trace shows the exact analytical solution, derived from the closed-form expression (27), while the blue trace depicts its numerical counterpart.(f) Impact of the positive DC voltage V + on the width w k (V + ) of the respective SET SDR, measured as the distance between the state values x +,k and x −,k , at which g SET (x, V + ) appears to be scaled down by a factor k as compared to its peak value g SET (x max , V + ) , for each k value from the set {1.5, 2, 3} .The exact analytical solution, descending from the formula (35), (The numerical solution) is illustrated through a dashed (solid) trace with red (blue), magenta (black), and green (brown) hue for the first, second, and third k value from the triplet.When 1.5, 2, and 3 is assigned to k, the numerical solution deviates from the corresponding analytical one as soon as V descends below +0.184 , +0.211 , +0.237 V (increases above +0.937 ,+0.932 , and +0.925 V ), since then x − ( x + ) descends below (rises above) the lower (upper) bound x L ( x U ) in the state existence domain D. www.nature.com/scientificreports/which can be solved for V + , resulting in an approximate analytical formula for V + (x max ) , featuring the form in which the positive (negative) sign in front of the first (second) square root sign descends from the polarity of V + (from the monotonic increase of V + with x max , as inferable from the graphs along the right column of Fig. 4), and where the functions a 0 (V +,1 , V +,2 ) and a 1 (V +,1 , V +,2 ) are in turn defined as Let us define the error e between x max and its approximation x max ( Ṽ+ (x max ), V +,1 , V +,2 ) as in which the second addend on the right hand side calls for the use of the approximate formula Ṽ+ (x max , V +,1 , V +,2 ) for V + (x max ) , as given in Eq. ( 22), into the closed-form expression for x max (V + ) , reported in Eq. ( 17), with γ (•) denoting the function in (18).For each voltage parameter pair (V +,1 , V +,2 ) , we first com- puted the maximum of the squared error e 2 , as x max was swept in D , and then found the minimum of the result- ing list of numbers via According to this optimisation procedure, assuming V +,2 > V +,1 , the best voltage parameter pair (V +,1 , V +,2 ) , let us denote it as (V (opt) +,2 ) , was found to be equal to (0.662, 0.923) V , which delivered the lowest possible maximum squared error, amounting to 5.635 × 10 −8 .(see the rightmost red-filled circle in Fig. 10).The comple- mentary hypothesis V +,2 < V +,1 results in an optimal voltage parameter pair (V (opt) +,2 ) = (0.923, 0.662) V , whereby once again the maximum squared error descends to its global minimum (see the leftmost red-filled circle in Fig. 10).Using these optimal values for V +,1 and for V +,2 , the formulas for γ (V + , V +,1 , V +,2 ) and for Ṽ+ (x max , V +,1 , V +,2 ) are respectively plotted through red traces in Fig. 9b and c.In the latter plot the blue curve , and Figure 10.Surface of the maximum squared error max x max ∈D {e 2 (x max , V +,1 , V +,2 )} as a function of the voltage parameters V +,1 and V +,2 under optimisation.At each of the points (V +,1 , V +,2 ) = (0.662, 0.923) V and (V +,1 , V +,2 ) = (0.923, 0.662) V , marked as red circles, and symmetrically located relative to the plane V +,2 = V +,1 , the surface assumes the minimum possible value, specifically 5.635 × 10 −8 .Without loss of generality, in the remainder of this paper V +,2 is assumed to be larger than V +,1 .As a result the optimal parameter pair is chosen as (V www.nature.com/scientificreports/illustrates the numerical approximation for the inverse of the function in Eq. (17). Figure 9d depicts the approximate analytical formula (20) for +,1 , and V +,2 = V (opt) +,2 (red trace) together with the exact analytical closed-form expression for x max from (17) (blue trace).

Peak value of a SET SDR
The maximum value attained by the ẋ versus x locus for a given positive bias level V + assigned to the voltage v may be easily derived by inserting the analytical closed-form expression (17), with γ (V + ) expressed via Eq.( 18), in place for x max on the right hand side of (15), which employs (16) for α(x, V + ) .Algebraic manipulations allow to express the maximum g SET, max (V + ) of the SET state evolution function g SET (x, V + ) for a given choice of V + as where Figure 9e shows the g SET, max (V + ) versus V + locus, as extracted from the exact analytical formula (27) (red trace) and by means of a numerical procedure (blue trace).

Width of the Gaussian bell-shaped SET state evolution function
For any positive DC value V + assigned to the voltage v, the SET state evolution function g SET (x, V + ) features a gaussian shape on the ẋ versus x plane.Let x − (V + ) and x + (V + ) denote two state values, lying to the left and to the right of the abscissa of the maximum x max (V + ) of the gaussian function, respectively.Assume x −,k (V + ) and x +,k (V + ) to hold the same distance from x max (V + ) , and the common value of the SET state evolution function at each of these points to appear scaled down by a factor k relative to the maximum level g SET,max (V + ) , which may be expressed in mathematical terms as For a given choice of k ∈ R , let us now define the kth-scale width w k (V + ) of the gaussian function as the distance between x +,k (V + ) and x −,k (V + ) , i.e.Using Eqs. ( 15) and ( 27), the condition (29) at either state value x ∓,k (V + ) ∈ {x −,k (V + ), x +,k (V + )} can be recast as Employing Eqs. ( 16) and (28), defining and recalling the earlier specified formula (18) for γ (V + ) , the constraint (31) reduces to the second-order polynomial which can be easily solved for x ∓,k (V + ) , yielding Inserting (34) into Eq.( 30), the kth-scale width w k (V + ) of the gaussian function g SET (x, V + ) is then comput- able via which, as indicated, is surprisingly found to be independent of V + .Fig. 9f depicts the kth-scale width w k (V + ) of the gaussian bell-shaped SET state evolution function g SET (x, V + ) against V + for the first, second, and third k value from the set {1.5, 2, 3} , as computed through the exact analytical formula (35), delivering in turn the con- stant values 0.076, 0.100, and 0.126 (red, magenta, and green dashed traces, respectively) as well as by numerical means (blue, black, and brown solid traces, respectively).Having acquired key geometrical properties of a gaussian bell-shaped SET state evolution function g SET (x, V ) , and, particularly, the positive value V + to be assigned to the DC voltage V to program its peak level at a prescribed state value x max , as reported in section "Positive DC voltage for programming the abscissa of the peak of a SET SDR", and its kth-scale width w k , as described in section "Width of the Gaussian bell-shaped SET state evolution function", we are now in a position to present a systematic technique for choosing 2 • P − 1 parameters of a gen- eralised pulse train from the class illustrated in Fig. 1c, specifically V +,1 , V +,2 , . . ., V +,P , τ +,1 , τ +,2 , . . ., τ +,P , τ − , with V +,1 < V +,2 < . . .< V +,P , so as to induce the coexistence of a predefined number of asymptotic oscillatory solutions for the memory state of the ReRAM cell around prescribed levels ( P ∈ N >0 ), for a given RESET pulse height V − , chosen beforehand.
Our methodology envisages to endow the TA-SE (10) with as many stable equilibria as the number P of positive pulses over each cycle of the pulse train voltage signal, falling across the ReRAM cell.In particular, the proposed systematic procedure is calibrated so as to ensure that the graph of the leftmost scaled gaussian bellshaped state evolution function τ+,1 • g SET (x, V +,1 ) , corresponding to the first positive input pulse, which features the smallest height V +,1 , creates one and only one stable TA-SE equilibrium xeq,1 with the graph of the modulus of the scaled RESET state evolution function τ− • g RESET (x, V − ) , appearing on the right hand side of Eq. (12).It also guarantees that the graph of the scaled gaussian bell-shaped state evolution function τ+,j • g SET (x, V +,j ) , cor- responding to the j th positive input pulse, which exhibits height V +,j , forms a pair of TA-SE equilibria, specifically xeq,2•j−2 and xeq,2•j−1 , featuring an unstable and stable nature, respectively, with the graph of τ− • |g RESET (x, V − )| , for j ∈ {2, . . ., P} .The target of the methodology is to massage the aforementioned 2 • P − 1 parameters of the generalised pulse train so as to endow the resulting TA-SE (10) According to the TA-SDR analysis, the stable ones, endowed with odd labels, and referred to as xeq,1 , xeq,3 , . . ., xeq,2•P−1 , are expected to denote the mean values of the P admissible stable oscillatory solutions x 1,ss (t) , x 3,ss (t) , . . ., and x 2•P−1,ss (t) for the memory state x of the periodically-forced ReRAM cell.
In order for the ith scaled SET state evolution function, with i ∈ {1, 2, •, P} , to dominate over the other P − 1 terms in the sum, composing the SET TA-SE component ẋSET from Eq. ( 14), locally, around the respective maxi- mum, which is a necessary critical measure to ensure that the existence of TA-SE equilibria in the region around x max,i is determined mainly by the interaction between the τ+,i • g SET (x, V +,i ) and the τ− • |g RESET (x, V − )| versus x loci, the P stable equilibria to be provided as a design specification to the input of the systematic procedure need to hold a suitable distance, which is at least one kth-scale width w k , one from any adjacent other.Moreover, for each i value from the set {1, 2, . . ., P} , the abscissa x max,i of the maximum of the gaussian SET state evolution function g SET (x, V ) , sampled at the DC voltage V +,i , which corresponds to the height of the ith positive pulse within each period of the train stimulus, is placed to the left of the prescribed (2 • i − 1) th stable TA-SE equi- librium xeq,2•i−1 , at an appropriate distance, amounting to one quarter of the kth-scale bell width w k , from its location.This step ensures that for each i ∈ {1, 2, . . ., P} the SET (RESET) forces win over the RESET (SET) ones to the left (right) of the (2 • i − 1) th TA-SE equilibrium xeq,2•i−1 , which, as a result, acquires a stability nature, as explained in section "The time average state dynamic route technique".Having computed the state value, at which the peak of the ith gaussian bell should appear, via x max,i = xeq,2•i−1 − w k /4 , for i ∈ {1, 2, . . ., P} , the approximate closed-form expression (22), with +,2 , is then employed to compute the positive value V +,i to be assigned to the DC voltage V in the expression for the SET state evolution function g SET (x, V ) , appearing in the ith addend of the sum to the right hand side of Eq. ( 14), i.e. the height of the ith positive pulse over each cycle of the train excitation signal v S .P algebraic equations are then written down to enforce the TA-SE (10) to feature equilibria at xeq,1 , xeq,3 , . . ., and xeq,2•P−1 .More specifically, these constraints, imposing an equality between the moduli of the SET ẋSET and RESET ẋRESET TA-SE components at each of the stable equilibria, which Eq. ( 10) is expected to admit, read as where r +,1 τ +,1 /τ − , r +,2 τ +,2 /τ − , . . ., and r +,P τ +,P /τ − express the first, second, . . ., and Pth SET to RESET pulse width ratio, respectively.This set of P equations is then solved for the unknowns r +,1 , r +,2 , . . ., and r +,P .The last unknown parameter, specifically the RESET pulse width τ − , which automatically fixes all the SET pulse widths τ +,1 , τ +,2 , . . ., and τ +,P , is finally chosen so small as to guarantee the accuracy of the predictions, drawn from the TA-SDR analysis, as verifiable through the investigation of the SCPCM of the periodically-forced memristive system, as well as via numerical simulations.

Remark 2
Even though an adequate distance between adjacent TA-SE equilibria is preliminarily observed in their prescription, as specified above, for an arbitrary choice of the negative input pulse height V − , out of the theoretic methodology, presented in this section, the leftmost scaled gaussian bell-shaped τ+,1 • g SET (x, V +,1 ) versus x locus may cross the graph of τ− • |g SET (x, V − )| as a function of x a couple of additional times, for some choices of the (36) of the first positive pulse and of the first SET-to-RESET pulse width ratio r +,1 .However, ad hoc control measures can be set in place to ensure that the interaction between the first scaled SET SDR and the modulus of the only scaled RESET SDR forms one and only one GAS equilibrium at the prescribed location xeq,1 for Eq. (10).For example, with reference to the numerical study, discussed in section "On the crossings between one scaled SET SDR and one scaled RESET SDR", and referring to the simple case, where a two-pulse-per-cycle pulse train stimulus is let fall across the ReRAM cell, choosing V − = −0.5 V , the resulting TA-SE equation admits one and only one GAS equilibrium xeq , irrespective of the pulse width ratio r, if V + is set to a value lower than the abscissa V+ = 0.494 V of the cusp in Fig. 8a.The state value x max , at which the g SET (x, V + ) versus x locus features a peak Figure 11.Illustrations elucidating how to choose the design parameter k for a case study, where it is requested for the ReRAM cell to act as a bistable oscillator under the application of a three-pulse-per-cycle pulse train voltage stimulus between its terminals.Let the ith positive pulse in the input sequence over each cycle have amplitude V +,i and width τ +,i , for i ∈ {1, 2} .The negative pulse, following the two positive ones in each input cycle, is assumed to feature a fixed amplitude V − of −0.5V , while its width τ − is to be determined.It is further required for the left LAS TA-SE equilibrium xeq,1 to be located at 0.280.The right LAS TA-SE equilibrium xeq,3 should be apart from the left one by one bell width w k .When k is set to 1.5, 2, and 3, xeq,3 is expected to lie at 0.356, 0.380, and 0.406, respectively.(a) For k = 1.5 the application of the design methodology first employs the approximate analytical formula (22), with x max set to x max,1 = xeq,1 − w 1.5 /4 ( x max,2 = xeq,3 − w /4 ), +,1 , and V +,2 = V (opt) +,2 , to fix the amplitude V +,1 ( V +,2 ) of the first (second) SET pulse to 0.483 V ( 0.550 V ).It then specifies the values 0.815 and 2.687 × 10 −5 for r +,1 and r +,2 , respectively, by solving the linear system of equations (36)-(37).Regardless of the choice for the RESET pulse width τ − , which automatically fixes the values for the SET pulse widths τ +,1 and τ +,2 , the TA-SE is found to admit the triplet of equilibria (x eq,1 , xeq,2 , xeq,3 ) = (0.132, 0.28, 0.356) .Clearly, the design specifications are not satisfied here.(b) For k = 2 , applying the proposed methodology delivers first the SET pulse heights V +,1 = 0.478 V , and V +,2 = 0.564 V , and then the SET-to-RESET pulse width ratios r +,1 = 10.866 and r +,2 = 8.974 × 10 −7 .The TA-SE equilibria are then found to lie at xeq,1 = 0.251 , xeq,2 = 0.28 , and xeq,3 = 0.38 .Also in this case the systematic procedure, introduced in this paper, fails to fulfil the design tasks.(c) Recurring to the proposed design methodology with k = 3 , the pulse train voltage stimulus is crafted as specified by the parameters V +,1 = 0.472 V , V +,2 = 0.580 V , r +,1 = 54.759, and r +,2 = 1.715 × 10 −8 .The TA-SE admits here the equilibria xeq,1 = 0.280 , xeq,2 = 0.309 , and xeq,3 = 0.406 .Therefore, choosing k = 3 , the combination between the two gaussian bells and the red curve, increasing monotonically with the time average state, allows to endow the TA-SE with two LAS equilibria at the desired locations, meeting the design requirements.Graphs revealing the instrumental role of the TA-SDR analysis tool to guide the circuit designer toward an appropriate choice for the parameter k for a case study, where a pulse train voltage stimulus, composed of one negative and three positive pulses per cycle, is expected to induce tristability in the oscillatory response of the ReRAM cell.Let V +,i ( τ +,i ) indicate the pulse amplitude (width) of the i th SET pulse, for i ∈ {1, 2, 3} .The pulse amplitude V − of the RESET pulse is fixed to −0.5 V , while its width τ − is an unknown variable.The leftmost LAS TA-SE equilibrium xeq,1 should lie at 0.275.The jth equilibrium xeq,j should appear to the right of the (j − 1) th equilibrium xeq,j−1 by as much as one bell width w k , for j ∈ {2, 3} .For k equal to 1.5, 2, and 3, the inner (rightmost) LAS TA-SE equilibrium xeq,3 ( xeq,5 ) is expected to lie at 0.351 (0.428), 0.375 (0.475), and 0.401 (0.527), respectively.(a) Choosing k = 1.5 , the proposed systematic design procedure first specifies the values 0.478 V , 0.546 V , and 0.606 V for the SET pulse amplitudes V +,1 , V +,2 , and V +,3 , respectively, via the approximate analytical formula (22), for V +,1 = V (opt) +,1 , and V +,2 = V (opt) +,2 , and fixing x max in turn to x max,1 = xeq,1 − w 1.5 /4 , x max,2 = xeq,3 − w 1.5 /4 , and x max,3 = xeq,5 − w 1.5 /4 .It then solves the system of linear Eqs.(36)-(38) with P = 3 for r +,1 , r +,2 , and r +,3 , in turn found to equal 13.228, 2.375 × 10 −5 , and 5.399 × 10 −12 .Irrespective of the choice for τ − , which directly sets values for τ +,i , with i ∈ {1, 2, 3} , the intersections between the loci of the moduli of the SET and RESET TA-SE components, identifying the equilibria xeq,1 , xeq,2 , and xeq,3 , the outer (the inner) of which are LAS (is unstable), for Eq. ( 10), are found to lie at 0.275, 0.351, and 0.428, respectively.As the TA-SDR analysis predicts bistability in the memristor steady-state oscillatory behaviour, assigning 1.5 to k is not an appropriate design choice.(b) For k = 2 , out of the proposed design procedure, the input parameters V +,1 , V +,2 , V +,3 , r +,1 , r +,2 , and r +,3 , are respectively set to 0.473 V , 0.560 V , 0.636 V , 31.913, 1.578 × 10 −6 , and 2.016 × 10 −16 .Correspondingly, the TA-SE admits the five equilibria xeq,1 = 0.275 , xeq,2 = 0.319 , and xeq,3 = 0.370 , xeq,4 = 0.375 , and xeq,5 = 0.475 , of which those labelled with odd numbers are LAS.Here the systematic parameter tuning procedure meets the design specifications.However the robustness of the design is questionable, given the non-ideal proximity between the TA-SE equilibria xeq,3 and xeq,4 . (c) With k = 3 , the application of the design procedure allows to choose the input parameters V +,1 = 0.467 V , V +,2 = 0.576 V , V +,3 = 0.668 V , r +,1 = 1.115 × for V + = V+ is 0.266.This directly sets the maximum value, which may be prescribed for the TA-SE equilibrium xeq , so as to ensure its GAS property, irrespective of r, to x max + w k /4 , that equals 0.285, 0.291, and 0.297, for the first, second, and third k value from the set {1.5, 2, 3} .In principle, as is the case for the examples illustrated in Figs. 13, 14, and 15, assuming a P-pulse-per-cycle pulse train voltage stimulus were let fall across the ReRAM cell, it is also possible to set the first TA-SE equilibrium xeq,1 to a value larger than this upper bound, but then, after solving the system of linear Eqs. ( 36)-(38), it would be necessary to check that the selection of values for the first pulse height and for the first SET-to-RESET pulse width ratio would fall in the monostability green region of the coloured r versus V + map, with r = r +,1 , under the specified value for V − .Finally, it is worth pointing out that a suitable change in the value, assigned to V − , may allow to move the abscissa V+ of the cusp, indicating the left bound of the red bistability domain, to the right, relative to its location along the horizontal axis in the coloured map of Fig. 8a.With reference to the proposed methodology, this would result in a corresponding increase in the maximum value, which may be prescribed for the stable equilibrium xeq,1 , that the leftmost scaled gaussian bell-shaped SET SDR would form with the graph of the modulus of the scaled RESET state evolution function over the time average state, irrespective of the first SET-to-RESET pulse width ratio r +,1 .

Remark 3
The selection of the real-valued parameter k is a critical design choice.In order to gain insights into this important aspect, Figs.11 and 12 illustrate two examples, where the methodological approach, presented in this section, is applied for different k values in the attempt to endow the TA-SE with two or three prescribed equilibria, respectively.In each of the two figures, plots (a), (b), and (c) show the loci of the moduli of the SET and RESET components of the corresponding TA-SE for the first, second, and third k value in the set {1.5, 2, 3} , revealing how only assigning the largest value in this triplet to the parameter under discussion allows to satisfy the design specifications (see the respective captions for more detail).
Importantly, under a proper selection for k, constraining the SET and RESET TA-SE components to comply with the set of P constraints (36)-( 38), together with the imposition of a minimum distance between adjacent stable equilibria, prescribed for the TA-SE, as well as with a sufficient leftward shift of each SET SDR relative to the respective stable TA-SE equilibrium, the scaled gaussian bells, resulting from the application of the proposed algorithm, gracefully pass over the locus of the modulus of the scaled RESET state evolution function versus the time average state in the regions of the respective peaks only, as may be inferred from either of Figs.11c and  12c, which refer to a particular case study for P = 2 and for P = 3 , respectively, and where, as a result, the blue- coloured single-valued curve, illustrating the TA-SE component, is found to oscillate around the graph of the modulus of the RESET TA-SE component as a function of the time average state, creating 2 • P − 1 equilibria, of which P stable, as prescribed, for the TA-SE.In each of the case studies, illustrating the application of the theory in section 6, keeping such a value for k, which implies a minimum distance between adjacent prescribed stable TA-SE equilibria of w 3 = 0.126 , and a spacing between each prescribed stable TA-SE equilibrium and the abscissa of the peak of the respective gaussian bell of w 3 /4 = 0.031 , proves to be a suitable choice to accomplish a robust design.

Discussion
The first part of this section applies the rigorous system-theoretic methodology, presented in section "A systematic methodology to craft the pulse stimulus for enabling the ReRAM cell to support multiple oscillations around prescribed resistance levels", to the Strachan model 4 for the determination of heights and widths of all the pulses, appearing cyclically across the ReRAM cell, so as to endow it with three, four, or five coexisting oscillatory operating modes around prescribed resistance levels.The second part of this section is devoted to show an interesting potential application, where the local fading memory effects, emerging across the nonvolatile resistance switching memory under periodic pulse train stimulation, could be leveraged to counteract certain non-idealities, which may be responsible for the corruption of the synaptic weights stored in a crossbar array.

Application of the theory to endow the ReRAM cell with three, four, or five oscillatory behaviours
The first, second, and third examples, illustrated in turn in Figs. 13, 14, and 15, result from the application of the theoretical method, presented in section "A systematic methodology to craft the pulse stimulus for enabling the ReRAM cell to support multiple oscillations around prescribed resistance levels", to the Strachan model 4 for the specification of suitable values for the 2 • P − 1 tuneable parameters of a generalised pulse train volt- age stimulus, belonging to the class, visualised in Fig. 1c, namely V +,1 , V +,2 , . . ., V +,P , τ +,1 , τ +,2 , . . ., τ +,P , τ − , with V +,1 < V +,2 < . . .< V +,P , when V − is preliminarily set to −0.5 V , so as to induce the coexistence of P stable asymptotic oscillations with prescribed mean values xeq,1 , xeq,3 , . . ., xeq,2•P−1 in the memory state of the periodically-forced ReRAM cell, with P set to 3, 4, and 5, respectively.

Remark 4
The theoretical framework, presented in this manuscript, provides evidence for the support, which nonlinear system theory may provide to experimenters and circuit design engineers.The experimental validation of the theory is the aim of our future research efforts.There are several challenges to tackle in order to achieve this goal.
1. Memristor devices available today can have limited endurance and their electrical behaviour may be subject to subtle drifts under operation, requiring much care and numerous repetitions to acquire convincing results.2. Intrinsic variability in memristors requires the procurement of significant statistics, regarding device-todevice and cycle-to-cycle variability effects, for the provision of convincing experimental results.3. The input pulse sequences, required in our programming schemes, are rather complex, requiring finelyprogrammable high-frequency pulse generators to support the experimental validation activities.This calls for the need to adapt existing measurement routines, available in house, or to acquire new experimental setups.
In regard to the third challenge from the above list, it might finally turn out to be less problematic than it seems, as explained next.The application of the rigorous system-theoretic methodology, presented in this section, to the Strachan model results in the specification of input pulse widths, decreasing at exponential rate with increases in their heights.While this issue does not undermine the significance of the theoretical work, which is applicable mutatis mutandis to any other memristor model, it originates here as the Strachan mathematical description 4  was not optimised for regions of the state-voltage space, where the ReRAM cell undergoes local fading memory effects, supporting multistable oscillatory operating modes.In fact, in these regions-refer to the order of magnitude of the bell peak value in either of plots (f), (h), and (l) of Fig. 4, extracted from the DRM upon assigning the first, second, and third positive value V + from the set {0.6 V, 0.8 V, 1.0 V} to the DC voltage V-the Strachan model may overestimate the speed of the oxygen vacancies as they move across the longitudinal extension of the nanodevice during a SET resistance switching process.This issue points to the necessity to retune the Strachan model so as to reproduce more accurately the behaviour of the nanodevice in the state-voltage space domain, where it is subject to local fading memory effects.Importantly, research investigations, applying the proposed systematic methodology to a recent reformulation of the Strachan model 21 , which employs ad hoc functions to limit to some extent the maximum admissible velocity, attainable by the ions under positive voltages, and was introduced to resolve some numerical issues, the original DAE set may suffer from, resulted in a dramatic increase in the minimum pulse width by several orders of magnitude relative to the case, where no upper bound was enforced on the rate of change of the memory state, in various scenarios, where the amplitudes assigned to the SET pulses were found to trigger local fading memory effects in the ReRAM cell.This provides proof-of-concept evidence that the application of our theory to a properly-optimised variant of the Strachan model might lead to the specification of widths and heights for the pulses, composing cyclically the train stimulus, which would be programmable in the control settings of existing physical AC voltage waveform generators.
(a) Decomposition of the TA-SDR into its SET (blue trace) and RESET (red trace) contributions, here plotted together on the | ẋ| versus x plane to visualise each possible equilibrium xeq of equation (10), where ẋSET = − ẋRESET , for a case study, where the proposed methodology from section 5.3 set the values for the parameters V +,1 , V +,2 , V +,3 , r +,1 , r +,2 , and r +,3 of a four-pulse-per-cycle pulse train voltage stimulus, with V − preliminarily chosen as −0.5 V , to +0.490 V , +0.649 V , +0.778 V , 4.594, 1.489 × 10 −18 , and 1.361 × 10 −47 , respectively, so as to endow the TA-SE with the 3 stable equilibria xeq,1 = 0.3 , xeq,3 = 0.5 , and xeq,7 = 0.7 , which in turn place the maxima of the first, second, and third gaussian bells at x max,1 = 0.269 , x max,2 = 0.469 , and x max,3 = 0.669 .The TA-SE equilibria are found to lie at xeq,1 = 0. ◂ of CMOS circuitry, promise to overcome the performance limitations of traditional technical systems, opening a wide spectrum of opportunities for electronics in the post-Moore era.Due to the strong nonlinearity, characterising the operating principles of these nanodevices, recurring to powerful concepts from Nonlinear Circuit and System Theory 2 is a necessary step for drawing a full picture of their dynamics.In fact the common approach of electrical engineers to linearise the model of a nonlinear device before commencing the investigation of its dynamics is insufficient to explore their global behaviour.As an example of the significant impact that this theory may have on the progress of memristor research, this paper reveals how the application of some of its powerful techniques to a predictive model 4 of a Ta 2 O 5−x Resistive Random Access Memory cell from Hewlett Packard Labs may allow the development of a systematic strategy, supported by a rigorous analytical framework, to craft a generalised rectangular pulse train voltage stimulus, composed of P ∈ N >0 SET positive pulses and of a single RESET negative pulse, so as to endow the memory state of the nano-device with P of coexisting oscillatory solutions, revolving around mean values, prescribed as design specification, and observable at steady state for all initial conditions drawn from their basins of attraction.The availability of an algorithm, which, evaluating analytical formulas, and solving a linear system of equations, automatically massages the properties of a generalised pulse train stimulus for triggering a monostable (multistable) periodic response in a Resistive Random Access Memory cell, triggering the emergence of global 3 (local 9,10 ) fading memory effects across its physical medium, and forcing it to oscillate around a specific resistance level, for any initial condition from the state existence domain (from a certain basin of attraction), may inspire the development and circuit implementation of novel in-memory sensing and computing paradigms in the years to come.As an example of a potential application of the theory, the local fading memory effects, emerging in the ReRAM cell according to the Strachan model, have been leveraged to propose a novel scheme to compensate for the unavoidable drift in the resistance of a crosspoint nanodevice under power off conditions.A similar theoretical approach, as the one, presented in this paper for the Strachan model, may be developed to investigate the response of the mathematical description of any other non-volatile 22 or volatile 23 resistance switching memory to periodic pulse train stimuli.to follow the voltage stimulus v S from (c) at all times, and for each initial condition x 0 from the set {x 0,1 , x 0,2 , x 0,3 , x 0,4 , x 0,5 , x 0,6 , x 0,7 , x 0,8 , x 0,9 , x 0,10 } = {0.15,0.34, 0.35, 0.485, 0.49, 0.62, 0.625, 0.745, 0.75, 0.9} .When initiated from either initial condition in the first, second, third, fourth, and fifth pair, the memristor state x converges progressively toward the steady-state waveforms x ss,1 , x ss,3 , x ss,5 , x ss,7 , and x ss,9 , respectively, as illustrated in plots (f), (g), (h), (i), and (l), which further visualise in turn the mean values x1,ss , x3,ss , x5,ss , x7,ss , and x9,ss of the asymptotic oscillations, together with the corresponding stable map fixed points x * 1 , x * 3 , x * 5 , x * 7 , and x * 9 .

Figure 1 .
Figure 1.(a) Circuit employed to investigate the response of the Ta 2 O 5−x ReRAM cell 15 to periodic square pulse-based voltage excitation signals.(b) Time course of a two-pulse-per-cycle train voltage stimulus.(c) Time waveform of a generalised pulse train voltage stimulus v S , including P positive SET pulses and one negative RESET pulse per cycle.The RESET pulse of height V − and width τ − follows the series of SET pulses.The i th SET pulse is V +,i high and τ +,i wide, with i ∈ {1, . . ., P} .The ordering of the positive pulses from the lowest to the highest in each input cycle follows the convention adopted in the systematic methodology to engineer multistability in the steady-state oscillatory response of the ReRAM cell to a generalised train stimulus (refer to section "A systematic methodology to craft the pulse stimulus for enabling the ReRAM cell to support multiple oscillations around prescribed resistance levels".However, this has no effect on the simulations.In fact, to facilitate their convergence, in the numerical investigations, discussed in section "Conclusions", the SET pulses were listed from the most narrow to the most wide before being applied in this order one after the other across the device.
https://doi.org/10.1038/s41598-024-55255-7www.nature.com/scientificreports/with τ+ τ + /T and τ− τ − /T .This ODE, referred to as time average state equation (TA-SE) 8 , governs the time evolution of the time average state x of the memristor when a voltage source, generating a specific square pulse train voltage stimulus v S , belonging to the class illustrated in Fig.1b, and characterised by the parameter quartet (V + , τ + , V − , τ − ) , is connected between its terminals, as shown in plot (a) of the same figure.Equations (11) and (12) are respectively referred to as SET and RESET TA-SE components.The blue (red) trace in Fig. 2a illustrates qualitatively a ẋSET ( ẋRESET ) versus x locus of the ReRAM cell subject to an arbitrary pulse train volt- age stimulus.The SET (RESET) resistance switching process tends to increase (decrease) the time average state over each input cycle, as indicated by the arrows on the first (latter) single-valued curve.Plotting the RESET TA-SE component in modulus, as depicted in plot (b) of the same figure, allows to visualise clearly each point, at which the SET and RESET forces balance out.A point of this kind denotes an equilibrium x = xeq for the TA-SE, as | ẋSET | x=x eq = | ẋRESET | x=x eq implies ẋ = 0 .A TA-SE equilibrium is asymptotically stable if and only if | ẋSET | > (<)| ẋRESET | locally to the left (right) of its location, and unstable otherwise.A blue filled circle (red hollow circle) is employed to indicate the location of a stable (an unstable) TA-SE equilibrium.The ẋ versus x locus, derivable by summing the ordinates of the vertically-aligned points, sitting along the SET and RESET traces from plot (a), for each x , as dictated by Eq. (

Figure 2 .
Figure 2. (a) Blue (Red) trace: SET ẋSET (RESET ẋRESET ) component of the TA-SE (10) of the ReRAM cell under an arbitrary pulse train stimulation.(b) Moduli of the SET and RESET TA-SE components.Their intersections identify the TA-SE equilibria.(c) TA-SDR of the ReRAM cell subject to the arbitrarily chosen pulse train stimulus.Arrows, pointing to the east (west), are superimposed along any TA-SDR branch, which visits the upper (lower) half plane, so as to indicate a progressive increase (decrease) in the time average state x when ẋ is positive (negative).An equilibrium for the TA-SE exists at the abscissa x = xeq of any point, at which the TA-SDR crosses the horizontal axis, as ẋ = 0 therein.The equilibrium is asymptotically stable (unstable), as indicated via a black filled (red hollow) circle, if and only if the slope ∂ ẋ/∂ x of the ẋ versus x locus is negative (positive) at its location.According to the TA-SDR analysis, the ReRAM cell is expected to act as a bistable oscillator under the given periodic excitation.
, a periodic stimulus endows the memory state of the ReRAM cell with two locally asymptotically stable (LAS) steady-state oscillatory solutions.A black filled (red hollow) circle is employed to indicate the location of a stable (an unstable) fixed point of the map.For each k value in N >0 the SCPCM expresses the net change �x k;k−1 x k − x k−1 = P (x k−1 ) − x k−1 , which the ODE solution x undergoes, over the time interval

Figure 3 .
Figure 3. (a)Exemplary illustration of a one-dimensional discrete-time system x k = P (x k−1 ) , referred to as Poincaré map, which admits three intersections with the identity map x k = P I (x k−1 ) = x k−1 , representing its fixed points, specifically x * 1 , x * 2 , and x * 3 , of which the outer ones (inner one) are stable (is unstable).A few coloured zig-zag trajectories, known as cob-web plots12,20 in Nonlinear Dynamics Theory, are also displayed to show the discrete-time evolution of the map from distinct initial conditions toward one of the two LAS fixed points.In our study a map of this kind can be extracted from the Strachan DAE set, when the input voltage v is enforced to follow a given periodic voltage stimulus v S , e.g. in the form of a rectangular pulse train, by recording samples of the memristor state x at regular T-long time intervals from the beginning of each of a large ensemble of simulations, differing in the initial conditions, and then plotting for each of the resulting time series the kth sample x k = x(k • T) versus the (k − 1) th one x k−1 = x((k − 1) • T) , with k ∈ N >0 .For k = 1 the SCPCM reduces to � 1;0 = x 1 − x 0 = P (x 0 ) − x 0 , providing the change in the memory state over the first input cycle. (b) x k;k−1 = x k − x k−1 versus x k−1 locus, illustrating the SCPCM of the ReRAM cell subject to the periodic stimulus, which induces a state motion resulting in the Poincaré map shown in plot (a).

Figure 4 .
Figure 4. (a), (c), (e), (g), (i) ((b), (d), (f), (h), (l)) g RESET (x, V ) ( g SET (x, V ) ) versus x locus, denoting the RESET (SET) SDR of the ReRAM cell15 when V is chosen as the first, second, third, fourth, and fifth value from the set {−(+)0.2, −(+)0.4,−(+)0.6,−(+)0.8,−(+)1} V.Over a RESET (SET) resistance switching transition the device state undergoes a progressive decrease (increase), as the arrows, superimposed on top of the respective SDR, clearly indicate through their westward (eastward) direction.With reference to each graph along the first column, the red filled circle shows the location of the stable equilibrium x eq = x L , which the ODE (1) admits for any negative bias value V assigned to the input voltage v. On the other hand, the state equation admits no equilibrium under any positive DC stimulus.

Figure 5 .
Figure 5. Three-dimensional illustration, showing each admissible stable or unstable equilibrium xeq = xeq (V + , V − ) , which the TA-SE(10), associated to a train voltage stimulus, featuring two pulses of opposite polarity per cycle, may possibly admit, when the SET τ + and RESET τ − pulse widths are identical, as a function of the SET V + and RESET V − pulse heights, swept across the ranges [−2, 0] V and [0, 1.2] V , respectively.The dark blue surface includes all the GAS equilibria of the TA-SE in the monostable oscillatory operating mode of the ReRAM cell.The cyan (magenta) surface contains all the unstable (all the LAS) equilibria of the TA-SE in the bistable oscillatory operating mode of the ReRAM cell.(b) Projection of the surface from (a) onto the V + versus V − plane.Choosing the pulses' heights of the pulse train voltage stimulus, featuring a 50% duty cycle, according to the coordinates of any point in the green (red) region, the TA-SE features a single GAS equilibrium (two LAS equilibria) for r = 1 .The black cross marker (black plus sign) identifies the input parameter pair (V − , V + ) , inducing the particular monostable (bistable) oscillatory response, illustrated in Fig.6(Fig.7), in the nanodevice.

Figure 6 .
Figure 6.SET | ẋSET | (blue trace) and RESET | ẋRESET | (red trace) components of the TA-SDR of the ReRAM cell under the application of a two-pulse-per-cycle pulse train voltage stimulus v S , when its SET V + and RESET V − pulse heights are in turn set to +0.46 V and −0.4 V , and for r = 1 , irrespective of the choice of its SET τ + and RESET τ − pulse widths.Note that scaling the widths of the 2 pulses in the train per cycle by the same factor does not affect the TA-SDR prediction.The only GAS equilibrium xeq of the TA-SE lies at 0.308, which is the abscissa of the black-filled circle.A marker, indicating the zero of the RESET component at x = 0 , is omitted from the graph, so as to avoid clutter.(b) SCPCM of the ReRAM cell subject to a particular pulse train voltage stimulus v S , belonging to the class considered in (a), and characterised by parameters (V + , τ + , V − , τ − ) = (+0.46V, 1µs, −0.4 V, 1µs) (refer to the blue signal of period T = τ + + τ − = 2 µs in plot (d)).The Poincaré map, from which it is extracted, features a GAS fixed point x * (see the black-filled circle).Differently from what is the case for the TA-SDR, scaling the widths of the 2 pulses in the train per cycle by the same factor may affect the SCPCM.(c) Brown (Green) trace: progressive approach of the solution x to the Strachan DAE set, when v is forced to follow the particular excitation voltage signal v S , employed for the derivation of the SCPCM, from the initial condition x 0 = x 0,1 = 0.15 ( x 0 = x 0,2 = 0.85 ) toward a unique steady-state oscillation.(d) Green trace: steady-state time series x ss of the memristor state x, as extracted from the solution featuring the same colour in plot (c).Horizontal lines mark the locations of the map fixed point x * , of the TA-SE equilibrium xeq , and of the time average xss of the steady-state time series.As the RESET pulse follows the SET pulse over each cycle of the input train, x ss attains its minimum value at the end of any period.Therefore x * directly reveals the minimum of x ss across one input cycle.

Figure 7 .
Figure 7. Decomposition of the TA-SDR into its SET | ẋSET | (blue trace) and RESET | ẋRESET | (red trace) components-plot (a)-for the ReRAM cell subject to a two pulse-per-cycle pulse train voltage stimulus v S , composed of one SET (RESET) pulse of positive (negative) amplitude V + = +0.54V (V − = −0.6V) over the first (second) τ + (τ − )-long half of each period of duration T = τ + + τ − , irrespective of the common value assigned to τ + and τ − .The TA-SE admits a triplet of equilibria, namely xeq,1 = 0.106 , xeq,2 = 0.237 , and xeq,3 = 0.370 .Each of the outer ones (The inner one), indicated via a black-filled (red hollow) circle, is LAS (unstable).(b) Time waveform of a particular pulse train voltage stimulus, belonging to the class assumed in (a), and identified via the parameter quartet (V + , τ + , V − , τ − ) = (+0.54V, 20 ps, −0.6 V, 20 ps) .(c) SCPCM of the ReRAM cell in the case, where the excitation voltage signal v S from (b) is let fall continuously between its terminals.A black-filled (red hollow) circle denotes a locally-stable (an unstable) fixed point for the associated Poincaré map.(d) Cyan (Violet) trace: time course of the memory state x of the ReRAM cell, with voltage v forced to follow v S from (b) at all times, from the initial condition x = x 0,1 = 0.2 ( x = x 0,2 = 0.3 ).Unlike the latter solution, the first one takes a very long time to attain the steady state.(e, f) Locally-stable oscillatory solution x 1 ( x 3 ) for the state x of the ReRAM cell, as recorded in a numerical simulation of the Strachan DAE set under v = v S from (b) for x 0 = x * 1 ( x 0 = x * 3).In each of the two cases the choice of the initial condition ensures that no transients appear in the device response.The time average x1 ( x3 ) of the solution x 1 ( x 3 ), as well as the corresponding LAS TA-SE equilibrium xeq,1 ( xeq,3 ) and LAS map fixed point x * 1 ( x * 3 ) are also marked in plot (e, f).

Figure 8 .
Figure 8.(a) Coloured map, depicting how the number of admissible stable or unstable equilibria for the TA-SE of the ReRAM cell, subject to a two-pulse-per-cycle pulse train voltage stimulus from the class illustrated in Fig.1bis influenced by the SET pulse amplitude V + as well as by the ratio r between the SET and RESET pulse widths, given a RESET pulse amplitude V − of −0.5 V .The green and red regions respectively enclose input parameter pairs, which endow the TA-SE with one and only one GAS equilibrium (three equilibria, of which the outer ones are LAS).(b, c) Graphical illustration, showing the decomposition of the TA-SDR into its SET and RESET components for a scenario, where the input pair (V + , r) , lying at (+0.50 V, 1 × 10 8 ) ((+0.75 V, 1 × 10 −30 )) (see the black cross marker (black plus sign) within the green (red) region of the map in (a)), determines the existence of one and only one GAS equilibrium xeq = 0.314 (three equilibria xeq,1 = 0.042 , xeq,2 = 0.516 , and xeq,3 = 0.725 , of which the outer ones are LAS) for the respective TA-SE.