In recent years, a surge of interest has been witnessed in exploiting nonlinear dynamical phenomena in emerging devices for novel circuit applications, such as neuromorphic computing. A subject that has been intensively studied is one-port (two-terminal) passive memristors, which exhibit a pinched hysteresis that always passes through the origin in their current–voltage (I–V) loci, thereby possessing a non-volatile memory effect (Strukov et al., 2008; Dittmann and Strachan, 2019). Passive memristors offer a scalable and energy-efficient approach to emulating biological synapses and implementing computationally efficient neuromorphic learning rules (Kim et al., 2015; Wang et al., 2017; Covi et al., 2018; Brivio et al., 2021).

Although a canonical memristor is a passive circuit element, any one-port device that exhibits a pinched hysteresis is considered an extended memristor (Chua, 2014), which includes a class of one-port devices that exhibit non-monotonicity in their experimental quasi-direct current (quasi-DC) I–V curves—a negative differential resistance (NDR). These devices typically exhibit a pronounced I–V hysteresis when driven by a voltage stimulus; however, the hysteresis collapses at a finite voltage. Therefore, they only have a transient (volatile) memory effect. Importantly, these globally passive one-port devices are locally active within the NDR region. Under proper biasing conditions, they can function as one-port amplifiers that increase the power of an alternating-current (AC) signal applied to the same port. It was recently shown that a locally active memristor (LAM) in a distributed form can act as an axon-like signal-amplifying transmission line (Brown et al., 2024). Signal amplification is an essential capability for information processing and communication, a field that has been dominated by transistors. Figure 1 shows the comparison of typical quasi-DC I–V curves measured from a passive memristor and a locally active memristor. Such a measurement varies V or I stimulus slowly and measures time-averaged device responses, which could capture the resistance states before and after a resistance switching event (see arrows) but not the ultrafast switching dynamics that could occur within femtoseconds. Figure 1A shows a bipolar tantalum oxide ( TaO y - Ta 2 O 5 ) passive memristor, which switches from a low-resistance state (LRS) to a high-resistance state (HRS) if a sufficiently large positive voltage is applied (Reset) and switches from the HRS back to the LRS if a sufficiently large negative voltage is applied (Set). Both Reset and Set operations are nonvolatile, i.e., resistance changes are retained even after power is turned off. In contrast, Figure 1B shows a unipolar vanadium dioxide ( VO 2 ) locally active memristor, which abruptly switches from an HRS to an LRS if a sufficiently large voltage is applied, regardless of its polarity. As the voltage is reduced below a smaller threshold level, the device switches from the LRS back to the HRS—these resistance changes are volatile and are lost once the power is turned off. It should be noted that the exemplary characteristics in Figure 1 are by no means exhaustive. Passive memristors can exhibit either bipolar or unipolar non-volatile resistance switching behaviors, which are determined by the intertwined ionic and electronic transport mechanisms within the nanoscale device volume (Waser et al., 2009; Jeong et al., 2012). They may also exhibit a fading memory effect, where the asymptotic behavior is solely determined by the state dynamics, irrespective of the initial condition (Ascoli et al., 2016; Pershin and Slipko, 2019).

We now focus on one-port devices that possess the peculiar I–V characteristics shown in Figure 1B. One-port devices with such switching characteristics have been extensively studied and implemented in engineering practice. They have been made out of many materials based on a variety of operating mechanisms. A familiar category is electro-thermal threshold switches such as ovonic threshold switches (OTSs), which show rapid changes in resistance due to nonlinear interactions among local temperature, metastable structural change, and electrical conductivity (Ovshinsky, 1968; Noé et al., 2020). Being a one-port device, LAMs and passive memristors share the same level of 4 F 2 (F: half pitch in lithography) scalability in a crossbar device geometry (Amsinck et al., 2005), resolving the trade-off between scalability and biological fidelity.

LAMs may generally be classified into two types, namely, current-controlled (S-type) and voltage-controlled ( λ -type), where the letters “S” and “ λ ” resemble the characteristic shape of NDR in the I–V curve plotted with voltage as the x -axis (Ridley, 1963). S-type LAMs are “normally off” devices with an HRS when the power is turned off. λ -type LAMs such as resonant tunneling diodes are “normally on” devices with an LRS when the power is turned off (Esaki, 1958). Therefore, S-type LAMs are superior to λ -type LAMs in terms of standby power consumption. Hereafter, we focus our discussion on current-controlled LAMs.

A particularly interesting class of current-controlled LAMs is based on the insulator-to-metal phase transition (IMT) phenomena in strongly correlated materials that arise from a coupling between structural distortions (Peierls transition) and electronic instabilities (Mott transition) (Andrews et al., 2019). They possess several attractive features for circuit applications, such as ultralow (femtojoule) switching energy (Prinz et al., 2020), ultra-fast (tens of femtosecond) switching speed (Jager et al., 2017), and electroforming-free operations (Yi et al., 2018). We term all these IMT-based LAMs “Mott memristors” without discerning the subtle differences in their phase transition mechanisms. Vanadium dioxide ( VO 2 ) and niobium dioxide ( NbO 2 ) are two intensively studied Mott memristor materials among many others (Andrews et al., 2019).

For neuromorphic computing applications, circuits of self-excited oscillators and spiking neuron emulators have been built with one or more LAMs that are coupled with reactive elements (capacitors) (Farhat and Eldefrawy, 1993; Moon et al., 2015; Ignatov et al., 2015; Stoliar et al., 2017; Wang et al., 2018). An illustrative example is a scalable spiking neuron, which constitutes two oppositely energized (“polarized” in neuroscience glossary) LAMs to mimic the voltage-gated potassium and sodium cell membrane ion channels. When coupled with parallel membrane capacitors and series load resistors, the composite circuit emulates a single-compartment nerve cell, initiating all-or-nothing action potentials upon a suprathreshold stimulus (Yi et al., 2018; Pickett et al., 2013) or acting as a delayed buffer, which allows bidirectional, distortion-free propagation of action potentials when daisy-chained (Pickett and Williams, 2013; Yi, 2022). Such a circuit topology bears some resemblance to a biologically plausible Hodgkin–Huxley (HH) axon model (Chua et al., 2012) and the early 1960s proposals of “neuristor” axons utilizing non-scalable components such as inductors (Crane, 1960; Crane, 1962; Nagumo et al., 1962). Experimental spiking neurons built with VO 2 Mott memristors have shown about two dozen biological neuronal temporal dynamics, including all three classes of neuronal excitability (Yi et al., 2018). Arguably, LAM-based neurons and passive memristor-based synapses form a self-sufficient basis to construct a transistorless neural network (Yi and Cruz-Albrecht, 2019).

Despite numerous experimental demonstrations, predictive modeling and analysis of LAM elements and circuits remain challenging and hinder technology development. These difficulties are, in part, due to the fundamental mathematical challenges associated with nonlinear differential systems. One illustrious example is the second part of Hilbert’s 1 6 th problem that questions whether there exists a finite upper bound for the number of limit cycles of planar polynomial differential systems. It remains unsolved today, even for quadratic polynomials (degree n = 2 ) (Ilyashenko, 2002). Qualitative local analyses, on the contrary, are facilitated by small-signal linearization techniques, which allow linear analysis to be applied to a nonlinear system near a hyperbolic fixed point with all eigenvalues of the linearization having non-zero real parts (Perko, 1991). A key theoretical contribution made by Chua is the local activity (LA) theorem, which provides a rigorous mathematical definition of the LA as a necessary prerequisite for the emergence of complexity in nonlinear systems (Chua, 2005). Moreover, Chua provided a set of explicit and computable criteria in the parameter space, which allows for identifying the edge-of-chaos (EOC) region that is both locally active and stable, where most of the complexity phenomena emerge. In recent years, Chua’s LA principle has been applied to clarify several long unsolved fundamental problems about dissipative systems, including Prigogine’s symmetry-breaking instability in homogeneous cellular media (Prigogine and Nicolis, 1967); the emergence of Turing instability (Turing, 1952) and its higher-order analog, Smale’s paradox (Smale, 1976), in reaction-diffusion systems; and Hodgkin–Huxley all-or-nothing excitability of nerve cells (Hodgkin and Huxley, 1952). All these complex phenomena are associated with the EOC domain within a system’s parameter space (Ascoli et al., 2022a; Ascoli et al., 2022b; Ascoli et al., 2022c; Chua, 2022).

Mathematically rigorous yet unphysical toy models of nonlinear dynamical elements were frequently used in the LA analysis procedure (Mannan et al., 2016; Mannan et al., 2017). For engineering practice, such toy models fall short of establishing a connection between circuit- or network-level dynamics and the measurable physical properties of constituent components. A recent review thoroughly elaborated on the importance of applying appropriate device physics into the mathematical memristor framework and defining physically relevant model parameters to control the circuit dynamic behavior (Brown et al., 2022b).

The main objective of this manuscript is to apply relevant theoretical techniques to understand the dynamics and stability of nonlinear circuits that involve locally active Mott memristors and map the conditions for the LA regime within the design parameter space (Messaris et al., 2021; Ascoli et al., 2021). These theoretical techniques include essential local analysis methods such as the small-signal linearization and the LA theorem and global techniques such as the nullclines and phase portraits. For engineering relevance, we base our analyses on an analytical one-dimensional (1D) Mott memristor compact model that is built on the laws of thermodynamics and only contains physically relevant device parameters. The model was developed by Pickett and Williams for NbO 2 Mott memristors (Pickett and Williams, 2012). Previously, we have verified that it is also applicable to VO 2 Mott memristors after replacing the material properties (Oh et al., 2010; Berglund and Guggenheim, 1969), and our SPICE simulations based on this model faithfully reproduced most of the measured neuronal dynamics in neuron circuits built with VO 2 memristors (Yi et al., 2018). In this study, we demonstrate that this physics-based compact model is mathematically tractable for applying local and global analysis techniques, with closed-form expressions for all the important quantities involved in the analyses. It enables a connection between the system dynamics and component physical parameters to guide circuit designs and process development. The algorithmic analysis procedure we present using a VO 2 Mott memristor model is general in nature and suitable for analyzing other Mott memristor materials. Qualitatively, the predictions regarding the dynamics and stability in this work are similar to those made by compact models based on different choices of state variables and kinetic functions (Brown et al., 2022b).

We focus on theoretical analyses and only included a cursory comparison between the model simulated and experimental characteristics of a VO 2 nano-device relaxation oscillator near the end. More detailed comparisons in the context of VO 2 Mott memristor neurons can be found in supplementary materials of Yi et al. (2018). It is understood that the compact model presented in this study has some simplifications and limitations. It is a nontrivial task to construct a computationally efficient compact model for locally active memristors with an appropriate balance between the dynamical fidelity and the computational complexity of solving the model equations. This is especially important for digital computer simulations of a scaled network that contains many instances of memristor elements, which could be costly in time and energy consumption.

The remainder of this paper is organized as follows: The first three sections (Sections 2–4) provide the analyses of an isolated 1D Mott memristor. In Section 2, we introduce the physics-based compact model and analyze the stability of an isolated 1D Mott memristor by examining its power-off plot (POP) and dynamic route map under constant input currents or voltages. This exercise confirms that the metallic state of a Mott memristor is unstable without power and is asymptotically stable with a finite input current. It also reveals that varying the voltage as the bifurcation parameter leads to a supercritical saddle-node bifurcation. In Section 3, we solve its locus of steady states (fixed points) in the three-dimensional (3D) state space and their two-dimensional (2D) projections. Note that we use both fixed point and steady state for the same concept in an interchangeable manner but avoid the term equilibrium unless the input is 0. See Subsection 2.2 for an elaboration on this topic.

In Section 4, we apply local analysis techniques on an isolated Mott memristor, including linearization and small-signal analysis, pole–zero diagram, Chua’s LA theorem, and frequency response. Its complex domain ( s -domain) equivalent circuit derived by the Laplace transform contains three virtual elements—a negative nonlinear capacitor in parallel with a negative nonlinear resistor, both in series with a positive nonlinear resistor. The negative s -domain capacitance gives rise to an apparent inductive response, similar to the memristive models of potassium and sodium ion channels (Chua et al., 2012). We found that an isolated Mott memristor near a fixed point dwells either in the locally passive (LP) or the EOC region. The EOC region coincides with the NDR region in its steady state or DC I–V locus. Brown et al. (2022b) explained that for an extended electro-thermal memristor, the coincidence between NDR and EOC or LA regions is not guaranteed. Therefore, NDR shall not be used as a sole signature for EOC. In our case, the crossover between the LP and EOC regions also manifests itself in the small-signal frequency response, which shows a sign reversal in the real part of the impedance (complexity) function Re Z ( s ; Q ) , as predicted by the fourth LA criterion. In the frequency domain, an isolated Mott memristor is equivalent to a positive inductor in series with a resistor that is positive in the LP region and negative in the EOC domain. We derived the parametric Nyquist plot of the LP ↔ EOC crossover at a single current level and then extended it to a 2D color-scale map of Re Z ( s ; Q ) to visualize the LP and EOC regions in the parameter space spanned by frequency and current, which is effectively a phase diagram for complexity. We also examined the scaling trend of the EOC region versus the device size, which shows that the VO 2 conduction channel radius is the relevant dimension for device miniaturization to enhance the EOC frequency regime.

Although an EOC region exists in an isolated 1D Mott memristor, the topological constraint limits the dynamics it can possess, making it impossible to exhibit persistent oscillations. In Sections 5–6, we remove the topological constraint for an isolated 1D Mott memristor by coupling it with one or more reactive elements, increasing the system’s dimensionality and dynamical complexity. For simplicity, we choose a DC voltage ( V dc ) -biased Pearson–Anson relaxation oscillator, formed using a Mott memristor coupled with a parallel capacitor C p and a series resistor R s , as an example of a 2D nonlinear system for our analysis (Pearson and Anson, 1921). The same analysis procedure can be applied to higher-dimensional systems, such as spiking neuron circuits consisting of two or more Mott memristors coupled with passive and reactive elements.

In Section 5 we first apply Chua’s LA criteria and local linearization techniques to this example system, including the element combination approach, the Jacobian matrix method, and the trace-determinant plane classification to study the stability and qualitative behaviors of its hyperbolic fixed points. The element combination approach considers a Mott memristor in parallel with a capacitor to be a composite second-order nonlinear element. The small-signal transfer function of the element-combined system has a pair of complex conjugate or real poles. We derived the Nyquist plot and a 2D phase diagram of the system’s poles. The pole phase diagram, combined with Chua’s LA criteria, allows a visualization of the LP, EOC and locally active and unstable (LA\EOC) regions in the circuit parameter space. These results are corroborated by the trace-determinant plane analysis of the Jacobian linearized 2D system, which reveals a stability-change bifurcation as the parametric (trace and determinant) locus crosses the zero-trace axis as one of the three circuit parameters is varied ( R s , C p , and V dc ). However, analyzing the stability behavior of a non-hyperbolic center requires additional theoretical tools since the Hartman–Grobman theorem is not applicable due to the loss of hyperbolicity (Hartman, 1960; Grobman, 1959).

Finally, in Section 6 we apply several global methods, such as the nullclines and numerical phase portrait analyses to understand qualitative behaviors of the non-hyperbolic centers in this example 2D nonlinear system. We found that each of the three circuit parameters ( R s , C p , and V dc ) acts as a bifurcation parameter that switches the stability of a fixed point as the parametric (trace and determinant) locus crosses a center. We verified that there exists a supercritical 2D Hopf-like bifurcation, i.e., the creation of a stable limit cycle encircling an unstable spiral as the fixed point switches its stability from stable to unstable. We also noticed that the limit cycle emerges abruptly over an extremely narrow bifurcation parameter interval, a phenomenon known as “canard explosion” in relaxation oscillations within chemical and biological systems (Krupa and Szmolyan, 2001; Rotstein et al., 2012). This is a prominent distinction from the classical Hopf bifurcation, which predicts a gradual growth proportional to the square root of the bifurcation parameter. Each bifurcation parameter has different bifurcation growth characteristics. We conclude the section with a comparison between the experimental limit cycle characteristics of a VO 2 relaxation oscillator and SPICE simulations based on the Mott memristor model, showing excellent agreements between them.

We conclude the manuscript with brief remarks on the application implications of locally active memristors and scalable neuromorphic dynamic neurons with a high degree of complexity.

In the present approach, the internal temperature is embedded in the biphasic model and is not considered a state variable. The set of all fixed points ( x Q , i Q , v Q ) in the 3D ( x , i , v ) state space that satisfy the instantaneous relationship v Q = R ch ( x Q ) i Q and ( d x / d t ) | Q = 0 is defined as the steady-state or DC locus of a Mott memristor. Solving the steady-state locus of an isolated Mott memristor is among the first steps for the local linearization analysis. Henceforth, both the ( x Q , i Q , v Q ) locus and its 2D projections are called the steady-state loci without discerning the dimensional difference.

To obtain the steady-state ( x Q , i Q , v Q ) locus, one can first define a sequence of i Q ∈ R and then find the solutions of the state variable x = x Q ( i Q ) numerically. This is achieved by setting the numerator in Equation 7 to be 0, which provides an equation C A ( 1 + B x Q 2 ) = − i Q 2 ⁡ ln x Q that can be solved numerically. After solving x Q ( i Q ) , voltage v Q can be calculated using the Ohm’s law relationship v Q ( i Q ) = R ch ( x Q ( i Q ) ) i Q .

However, there is a much easier way to obtain the steady-state ( x Q , i Q , v Q ) locus. Instead of numerically solving the value of x Q from a given i Q , one can first define a sequence of x Q ∈ ( 0,1 ) and then calculate i Q ( x Q ) analytically using the following formula: i Q x Q = − C A 1 + B x Q 2 ln x Q . (9)

Voltage v Q is then calculated using the Ohm’s law v Q ( x Q ) = R ch ( x Q ) i Q ( x Q ) . The sequence of x Q can be chosen to be evenly spaced on a linear or logarithmic scale, depending on how fast these functions vary with x Q . We verified that steady states calculated by both methods are consistent with each other. The analytical method is used for discussions hereafter.

Figure 7A shows the steady-state loci of ( x Q , i Q ) calculated using Equation 9 for three different VO 2 device sizes, plotted as x Q ( i Q ) , since Mott memristors are current-controlled devices. At small currents, the fraction of the metallic phase x Q remains negligibly small. x Q starts to increase with current in a sublinear fashion once i Q exceeds a size-dependent threshold level. The current threshold increases with the device size and is at μ A level for the shown device sizes.

Figure 7B shows the loci of the memristance function R ch ( x Q ) vs i Q , which reveal that R ch ( x Q ) has similar crossover characteristics at the same i Q thresholds. At small currents, R ch ( x Q ) remains elevated with negligible current dependence. Once i Q exceeds a size-dependent threshold, R ch ( x Q ) decreases rapidly with current in a nonlinear fashion. For the midsize VO 2 device ( r ch = 36 nm and L ch = 50 nm), R ch ( x Q ) decreases by more than three orders of magnitude from 122.8 k Ω to 97 Ω as i Q increases from 0 to 1 mA.

Figure 7C shows the steady-state loci of ( x Q , v Q ) plotted as v Q ( x Q ) , which resemble the shape of a left handled cup. The open left handle is nearly vertical. In other words, at very small x Q levels, a tiny change in x Q will cause a large change in v Q . Figure 7C (left inset) shows the plots of the extremely small x Q region of the ( x Q , v Q ) loci on a log–log scale, which reveals that at a given device size, there is a corresponding asymptotic lower bound of steady-state v Q as x Q approaches 0. For the midsize VO 2 device ( r ch = 36 nm and L ch = 50 nm), the v Q lower bound turns out to be 0.0973 V (dashed line). Figure 7C (right inset) shows the plot of the halfway x Q region in linear scale, illustrating that the v Q = 0.0973 V horizontal line is tangent to the ( x Q , v Q ) locus at its trough, located at Q 0 = ( 0.60628 , 0.0973  V ) (marked as × ); this corresponds to the same semi-stable fixed point Q 0 identified in the DR analysis. A slight increase in v Q would bifurcate Q 0 into a pair of fixed points on its left and right. The left inset also indicates that in this case, another fixed point would emerge at an extremely small x Q level (at v Q = 0.2 V, x Q is only 1 0 − 63 ), i.e., an insulating steady state exists at a finite voltage. These observations corroborate our previous DR analysis shown in Figure 6. All three ( x Q , v Q ) loci have a sharp peak at x Q = 0.00567 and a rounded trough at x Q = 0.60628 , resembling the shape of a cup. Notably, the x Q coordinates of these two extrema are size-independent.

Figure 7D shows the steady-state loci of ( i Q , v Q ) plotted as v Q ( i Q ) . As current-controlled memristors, the ( i Q , v Q ) loci are “N”-shaped when plotted with current as the x -axis. They are symmetric with respect to the origin in the first and third quadrants. Therefore, one only needs to analyze the first-quadrant halves. Each ( i Q , v Q ) locus has three distinctive regions: a lower positive differential resistance (PDR) region from 0 to a critical current i c 1 ; an NDR region between i c 1 and the second critical current i c 2 (see inset); and an upper PDR region for even higher currents. Therefore, i c 1 and i c 2 produce a local maximum and minimum in the ( i Q , v Q ) loci. For the shown device sizes, values of i c 1 ( i c 2 ) are 2.522 μ A (269.77 μ A), 9.077 μ A (971.18 μ A), and 14.122 μ A (1510.73 μ A), respectively. Figure 7D also shows that the steady state or DC loci of ( i Q , v Q ) always pass through the origin (0, 0), satisfying the zero-crossing property of memristors.

It should be noted that the volumetric enthalpy change in IMT Δ h tr appears only in the denominator of the kinetic function f x ( x , i ) via the coefficient E . Therefore, it has no effect in determining the steady-state ( i Q , v Q ) loci. The main effect of IMT on the shape of steady-state ( i Q , v Q ) is applied via coefficient B —the coefficient of the quadratic nonlinearity in the memristance function. Coefficient B is approximately the electrical resistivity ratio ρ ins ρ met between the insulating and metallic phases.

The sets of loci plotted in Figures 7A, C, D are 2D projections of the steady-state loci ( x Q , i Q , v Q ) in the 3D state space. Figure 8 shows the locus of ( x Q , i Q , v Q ) calculated for the midsize VO 2 device ( r ch = 36 nm and L ch = 50 nm). It resembles a twisted handle of a binder clip. The two open legs of the clip are rotated out of the plane defined by the looped clip head. Figure 9 provides a zoomed-in view of Figure 8 to visualize the low-current region of the same ( x Q , i Q , v Q ) locus, allowing its 2D projections onto the ( i , x ) , ( v , x ) , and ( i , v ) planes to be directly compared with the loci shown in Figures 7A, C, D, respectively.