{'fraction_non_alphanumeric': 0.05182184090668361, 'fraction_numerical': 0.014934858595487766, 'mean_word_length': 4.260481961275944, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mott Memristors based on Field-Induced Carrier Avalanche Multiplication

We present a theory of Mott memristors whose working principle is the non-linear carrier avalanche multiplication in Mott insulators subject to strong electric fields. The internal state of the memristor, which determines its resistance, is encoded in the density of doublon and hole excitations in the Mott insulator. In the current-voltage characteristic, insulating and conducting states are separated by a negative-differential-resistance region, leading to hysteretic behavior. Under oscillating voltage, the response of a voltage-controlled, non-polar memristive system is obtained, with retarded current and pinched hysteresis loop. As a first step towards neuromorphic applications, we demonstrate self-sustained spiking oscillations in a circuit with a parallel capacitor. Being based on electronic excitations only, this memristor is up to several orders of magnitude faster than previous proposals relying on Joule heating or ionic drift.


INTRODUCTION
In strongly correlated materials, many-body electronic interactions cannot be treated as a weak perturbation. A spectacular consequence is the breakdown of standard band theory in Mott insulators, which display a charge gap despite having nominally partially-filled bands. Even more interesting, from both fundamental and applied points of view, are states of matter obtained from a Mott insulator by applied pressure or chemical doping [1,2], photo-doping [3][4][5][6], or applied electric field [7,8].
The resistance of Mott insulators may vary over several orders of magnitude across different branches of the I-V curve. Owing to this resistive switch, Mott materials are promising candidates for replacing conventional semiconducting transistors in the field of information processing. More specifically, in neuromorphic applications [36] they are proposed to fabricate memristors [37][38][39][40], electronic devices whose resistance depends on the history of the input signal, which are regarded as the building blocks of bio-inspired novel computing architectures [41][42][43][44][45].
From a formal point of view, a voltage-controlled mem-ristive system is defined by its state-dependent resistance, or memristance M (x) (x: state variable) and by the equation of motionẋ = f (x, V ). The instantaneous resistance depends, therefore, on the past voltage. From a more empirical perspective, the fingerprint of a memristor is a pinched hysteresis loop in the I-V plane when the device is subject to a bipolar periodic signal [38,39]. Following semiconducting thin films with intertwined electronic and ionic motion [40], diverse other solid-state platforms are being investigated as physical realizations of memristors; in particular Mott materials, using Joule heating to locally trigger the insulating-to-metal transition [46][47][48][49][50][51]. The time scale of these devices is set by the physical mechanism for the resistance switch and is of the order of milliseconds for ionic drift [40,52] and of nanoto microseconds for Joule heating [47].
In this work we present a theory of a new type of memristor made of a narrow-gap Mott insulator, whose state variable is the density of doublon excitations, which are the charge carriers. In stark contrast with previous proposals, the resistance switch in this memristor is based on a purely electronic mechanism: the field-induced nonlinear carrier avalanche multiplication. This results, in particular, in a time scale set by the doublon decay time which is of the order of picoseconds, namely up to several orders of magnitude faster than in previous proposals.
In the following, we illustrate the microscopic working principle in Sec. I, where we present a phenomenological model for the field-induced non-linear carrier avalanche multiplication. Building on this, in Sec. II we introduce our model of Mott memristor, derive the static currentvoltage curve, and study the d. c. transitions between insulating and conducting states. In Sec. III we study the a. c. response, obtaining the typical behavior of a voltagecontrolled, non-polar memristive system; and derive the steady-state diagram. Finally, in Sec. IV, as a first step towards neuromorphic applications, we study a circuit with a parallel capacitor and demonstrate self-sustained current oscillations, reminiscent of the periodic spiking activity of biological neurons. We start by presenting a phenomenological model of Mott insulator as a material with variable concentration of charge carriers. In this model, similarly to electrons and holes in semiconductors, the carriers are doublons and holes, which are one-particle excitations in upper and lower Hubbard bands, respectively, see Fig. 1(a). Note that here we adopt a simplified description and do not consider the dynamical nature of the Mott gap, which is held fixed. Furthermore, we impose the doublon-hole symmetry, such that these excitations differ only for their charge (±e) and have the same concentration n, which hereafter is simply referred to as doublon density.
The density of doublon excitations n can be considered as a state variable which determines the conductivity of the material. In this phenomenological model, doublons have charge e, effective mass m * and they accelerate in an electric field, before scattering after a typical time τ . This is formalized in the Drude formula for the conductivity, which relates the current density j to the electric field E, Similar forms to Eq. (1) also apply to weakly correlated materials, with for example n representing the density of conduction-band electrons. The key difference with the model at hand is in the rate equation for n, in which the strong correlations typical of Mott materials appear as a non-linear term in the doublon density, Here the source term γ describes excitations of doublonhole pairs induced by thermal fluctuations or quantum tunneling across the gap, see Fig. 1(b). In principle, these depend on temperature and electric field; here we hold γ fixed and concentrate on the field dependence of the other terms. The second term in Eq. (3) describes the decay of doublon excitations with a typical time τ d [53]. The equilibrium density, namely the zero-field stationary solution, is n = γτ d ≡ n 0 . The one-body avalanche term (a 1 nE 2 ), also known as impact ionization, is present in both strongly [25] and weakly correlated materials [31]. It describes a process in which the kinetic energy of a carrier is converted into excitation energy of new carriers via scattering with impurities or phonons [ Fig. 1(c),(d)]. The two-body avalanche term (a 2 n 2 E 2 ), on the other hand, describes many-body scatterings of two excitations kicking out new carriers [ Fig. 1(e)-(g)] and is therefore proportional to the squared carrier density. The last term describes carrier diffusion due to density gradients; hereafter we consider the homogeneous case ∇ 2 n = 0. In nonzero electric field, Eq. (3) yields two stationary doublon densities, that is the solutions ofṅ = 0: (4) Here E 0 ≡ (τ d √ a 2 γ) −1 and A ≡ a 1 (a 2 γτ d ) −1 is the ratio of the one-to the two-body avalanche term for n = n 0 . Imposing the solutions (4) to be real and positive yields the condition E < E th , with the threshold electric field where the approximation is valid for small A, namely for predominant two-body avalanche. At this threshold, the two branches of Eq. (4) merge, the doublon density is and the current density reads In contrast with the threshold electric field and doublon density, the threshold current density does not depend on the one-body constant a 1 , but only on the two-body constant a 2 (through E 0 ) and it diverges for a 2 → 0.
Since the conductivity increases with doublon density, we can interpret the lower branch of Eq. (4) as the slightly perturbed equilibrium insulating state, and the upper branch as a conducting state. The corresponding current density j = σ(n(E))E is plotted in Fig. 2(a). It should be stressed that the two branches correspond to the same microscopic state and differ only in the doublon density; in particular, this theory does not cover the field-induced collapse of the Mott gap. Equation (4) also implies that, within this model, there are no stationary solutions for E > E th , meaning that the material cannot sustain such electric fields. In Fig. 2(b) we plot the conductivity as a function of the current density, whereĒ(j) is the inverse function of j(E) = σ(n(E))E and the approximation is valid for small A. Expressions similar to Eq. (8) have been suggested to explain experiments on a class of charge-transfer insulators [7,32]. The results in Fig. 2 are in qualitative agreement with experiments in which a current is passed through a Mott insulator and the electric field (thus the conductivity) is measured, see e. g. Refs. [34,35]. Indeed, up to this point the treatment is suitable to describe situations in which the current, and not the electric field, is the external parameter. To show this from a formal point of view, we linearize Eq. (3) around the stationary solution (4) at fixed E or at fixed j = σ(n(E))E. In the former case we shows that only the lower branch is stable. If we instead fix j, we get τ d δṅ = −δn[1 + A(j/j 0 ) 2 (n 0 /n) 2 ] which is stable for all current densities. Only in the latter case states with large conductivity are stable and can therefore be observed.
Among the parameters introduced in this section, most relevant are τ d , E 0 , j 0 ; which set the characteristic scales of, respectively, time, electric field, current density. The doublon decay time is typically τ d ∼ 1-10 ps, as measured in ultrafast pump-probe optical spectroscopy [3][4][5][6], while electric fields of the order E 0 ∼ 1-10 kV cm −1 and current densities j 0 ∼ 1-10 mA cm −2 have been measured in Refs. [7,34,35]. Together with the physical dimensions of the memristor, E 0 and j 0 also set the characteristic scales of, respectively, voltage and current.

II. CURRENT-VOLTAGE CHARACTERISTIC AND INSULATING-CONDUCTING TRANSITIONS
We introduce now our model of Mott memristor as a device composed of a Mott insulator connected in series with a conventional resistor. Adopting the description in Sec. I, the resistance of the Mott insulator is a function of carrier density through the conductivity σ(n) [Eq. (1)]: where L and S are length and section area. Instead, the conventional resistor has a fixed resistance R s . The total resistance of the memristor, or memristance, is therefore and the doublon density n is its state variable. Attaching a voltage generator V to the memristor, the electric field internal to the Mott material is with r s = R s /R 0 , R 0 ≡ R(n 0 ). Thus, the electric field does not depend solely on the applied voltage, but also on the doublon density. For small density the resistance of the Mott material is large, R(n) R s , and the field is approximately proportional to the voltage. On the other hand, for large density the resistance of the Mott material drops, R(n) R s , and so does the field. This mechanism is crucial for the stabilization of the conducting state of the memristor, as we discuss in this section.
The state-dependent resistance [Eq. (10)] and the rate equation for the state variable [Eqs. (3) and (11)] define a non-polar voltage-controlled memristive system [38]. In practice, the fixed term in the resistance corresponds to either the contact resistance, often present especially in two-probe measurements (see e. g. Ref. [33]), or a resistor added to obtain a stable conducting state [7,34,35].

A. Stationary doublon density
The stationary condition is obtained plugging Eq. (11) into Eq. (3) and imposingṅ = 0 (we set A = 0 hereafter). We solve the resulting equation for V : where V 0 ≡ LE 0 . This is plotted in Fig. 3(a) as n versus V (n) which allows us to visualize the stationary densitȳ n as a function of voltage. This solution is stable only if dn/dV > 0, namely forn outside a range [n * 1 , n * 2 ], where these values are therefore obtained imposing which yields For small r s we can approximate n * 1 ≈ 2n 0 (1 + 2r s ) and n * 2 ≈ n 0 r −1 s . Therefore, the two stable branches are well separated (n * 2 /n * 1 ≈ 0.5r −1 s ) and we can interpret them as the insulating (n < n * 1 ) and conducting (n > n * 2 ) states of the memristor. Increasing r s the two branches approach each other as n * 2 − n * 1 = n 0 √ 1 − 8r s /r s and eventually merge for r s = 0.125. Beyond this value, we have one continuous stable state with no clear separation between insulating and conducting states. In the opposite limit, r s → 0, the stable conducting branch vanishes (n * 2 → ∞). In the remainder of this work we set r s = 0.01. Between 5V 0 (1 + 2r s ) insulating and conducting states coexist. In particular, to V * 1 correspond the densities n * 1 on the insulating branch and n * 3 on the conducting branch.

B. Current-voltage characteristic
In the stationary state with voltageV (n) and doublon density n, the current through the memristor is where I 0 ≡ V 0 R −1 0 . Plotting Eq. (15) versus Eq. (12) we obtain the current-voltage curve in Fig. 3(b). This has a distinct "S" shape composed of three branches with alternating differential resistance R ≡ dV /dI, which is positive in the stable insulating and conducting branches; and negative in the unstable region in between [negativedifferential-resistance region (NDR)].
A voltage sweep across the range [V * 2 , V * 1 ] results in a current hysteresis, see Fig. 3(b). If the voltage change is adiabatic, meaning so slow that at each moment the memristor is stationary, then from the insulating branch the current follows the I-V curve up to V * 1 , where a jump discontinuity leads from I * 1 =Ī(n * 1 ) ≈ I 0 (1 + 2r s ) to the conducting branch in I * 3 =Ī(n * 3 ). Then, upon decreasing the voltage, the current remains large down to V * 2 where a second discontinuity leads from I * 2 =Ī(n * 2 ) ≈ I 0 / √ r s back to the insulating branch. If the voltage change is non-adiabatic, namely rapidly increasing and decreasing, the current does not follow thoroughly the I-V curve but instead traces a larger hysteresis area. In Fig. 3(c) we plot the same quantities as in Fig. 3(b) versus the electric field internal to the Mott insulator. Since current and current density are proportional, I = jS, the stationary curve is a rescaled copy of Fig. 2(a) with the crucial difference that this is now stable also for I > I * 2 . The trajectories appear different in the I-E plane with respect to the I-V curves; since during the constant-voltage insulating-conducting transitions both current and internal field vary. Also in this case, a nonadiabatic voltage results in a wider trajectory.

C. Delay time and relaxation time
To study the time scales associated with the transitions between insulating and conducting states, we consider a voltage with a ramp function f (t) = [1 + tanh(t − 10)]/2 and we numerically integrate Eqs. (3), (11). From the insulating state, as the voltage increases above V * 1 , the transition takes place in two steps [ Fig. 4(a),(b)]: first, during a delay time τ D the current remains low; then, it rapidly increases above I * 2 , meaning that the memristor has become conducting. Notice that after the transition I ∝ V f since in the conducting state the memristance is approximately constant The delay time is plotted in Fig. 4(c) versus the voltage and for varying initial conditions. While the insulting-toconducting transition naturally starts from the insulating branch [n 0 , n * 1 ], here we consider also initial conditions in the unstable region [n * 1 , n * 2 ] which are relevant in the case the voltage changes while the memristor is not at equilibrium. The delay time decreases with increasing voltage and larger initial density. It diverges in V * 1 if the initial density is below n * 1 ≈ 2.02n 0 , or inV (n(0)) otherwise. This difference can be explained with the aid of the stationary curve in Fig. 3(a), which shows that for n > n * 1 the minimum voltage leading to the conducting branch is indeedV (n). To get analytical insight into the delay time and its dependence on voltage and initial density, we solve Eq.
In this approximation the transition to the conducting state happens where Eq. (16) diverges, giving the delay time which we plot in Fig. 4(c) alongside the numerical result.
In the limit V → (V * 1 ) + we haven av → 2n 0 ≈ n * 1 and ∆ → 0. The behavior of τ D depends on whether n(0) is smaller or larger than n * 1 , in the former case it diverges as τ d ≈ 2τ d π/∆, while in the latter case it stays finite and diverges at a lower voltageV (n(0)).
Also the transition from the conducting state, as the voltage decreases below V * 2 , takes place in various steps [ Fig. 4(d),(e)]: first, the current rapidly decreases; then, it remains high during a relaxation time τ R ; finally, it decreases below I * 1 . The relaxation time is plotted in Fig. 4(f) versus the voltage and for varying initial conditions. Analogously to what discussed for the delay time, we consider initial conditions in the conducting branch [n 1 , ∞] as well as in the unstable region [n * 1 , n * 2 ]. The relaxation time increases with increasing voltage and larger initial density; and diverges in V * 2 if the initial density is above n * 2 ≈ 100 n 0 , or inV (n(0)) otherwise. The results in this section, in particular the currentvoltage characteristic and the delay time, qualitatively agree with various experiments on similar devices [7,8,[32][33][34][35]. Moreover, the analysis of delay and relaxation times sets the stage for the discussion of the a. c. response, a fundamental characteristic of a memristive system.

III. RESPONSE TO ALTERNATING VOLTAGE
We proceed now with the study of the a. c. response of the Mott memristor introduced in Sec. II and defined by its state-dependent resistance and state-variable equation of motion [Eqs.  3) and (11)] for various amplitude and frequency of the voltage V (t) = V a cos(Ωt) and for two different initial conditions. We distinguish four qualitatively different steady states. In Fig. 5(a),(b) the steady state is respectively insulating (n < n * 1 ) or conducting (n > n * 2 ) independently of the initial condition. In contrast, in the case of Fig. 5(c) there are two possible steady states depending on the initial condition. Finally, in Fig. 5(d) the steady state goes back and forth the insulating and conducting states.
The corresponding steady-state current is plotted in Fig. 5  In Fig. 5(i)-(l) we plot the steady-state current versus the voltage. This curve traces a pinched hysteresis loop (so called because it crosses the coordinate axes only in the origin) which is considered the empirical definition of a memristive system [38]. Also here, we have rescaled the axes for the sake of comparing the different steady states. The curve slope is the instantaneous inverse differential resistance R −1 = dI/dV , meaning that the greater the resistance change, the larger the area encircled by the loop. Indeed, this is more evident in the insulating state These correspond to, respectively, insulating state, insulating-to-conducting transition, conducting state; while the conducting-to-insulating transition happens near the origin [cf. arrows in Fig. 5(l)].
The direction of the loop, namely whether it is traced clockwise or anti-clockwise, is related to the polarity of the memristive system. In bipolar memristors, e. g. based on ionic drift [40], the resistance changes depending on the sign of the input. Consequently, it is either maximum or minimum in the origin of the I-V plane, and the loop is anti-clockwise for positive and clockwise for negative input. In contrast, in the present case the memristor is non-polar, meaning the resistance change is independent of the sign of the input, cf. Eq. (3). As a result, the loop is anti-clockwise both for positive and negative inputs [see arrows in Fig. 5(k),(l)]. Moreover, this implies that the slope in the origin, namely the zero-voltage inverse instantaneous resistance, is the same for increasing or decreasing voltage.
Other characteristics of a memristive system are more conveniently discussed in terms of the relation between charge q and flux φ, namely the integrals of, respectively, current and voltage. Indeed, originally the memristance was introduced as the quantity relating flux to charge [dφ = M (q) dq] similarly to how the resistance relates voltage to current [dV = R(I) dI] [37]. The steady-state charge-flux relation is plotted in Fig. 5(m)-(p). The multivaluedness of this relation is the empirical evidence that the memristive system belongs to the class of non-ideal memristors [38]. For ideal memristors, the state-variable equation of motion depends on the input only [ẋ = f (V )] giving a unique relation between charge and flux [37,40]. Instead, in the broader class of non-ideal memristors, the equation of motion depends also on the state variable itself [ẋ = f (x, V )] which yields a multivalued charge-flux relation, as in the this case. On a practical level, an ideal memristor is non-volatile, meaning its state does not change on zero input [f (0) = 0], while the state of a non-ideal memristor typically relaxes [f (x, 0) = 0] which makes it a volatile memory.

C. Steady-state diagram
In Fig. 6 we plot the steady-state diagram as a function of voltage frequency and amplitude. This contains four regions, delimited by the frequency-dependent amplitude thresholds V * 1,2 (Ω), corresponding to each of the steady states discussed above: 1. For amplitude smaller than V * 1,2 (Ω) (blue region in Fig. 6) the steady state is insulating, as in Fig. 5(a).

For frequency not too low and amplitude within
the range [V * 2 (Ω), V * 1 (Ω)] (purple region) the steady state is insulating or conducting depending on the initial condition, as in Fig. 5(c).

For low frequency and amplitude within the range
[V * 1 (Ω), V * 2 (Ω)] (green region) the steady state goes back and forth insulating and conducting, as in Fig. 5(d).
The a. c. thresholds V * 1,2 (Ω) are closely related to the d. c. thresholds V * 1,2 of Sec. II. If we apply a voltage with small amplitude, such that the memristor is insulating, and then gradually increase it, V * 1 (Ω) is the minimum value at which the memristor becomes conducting. Notice the analogy with V * 1 , which is the minimum voltage to trigger the d. c. insulating-to-conducting transition. However, in the a. c. case two scenarios are possible: the memristor either stays conducting indefinitely, or it goes back to insulating at a later point of the voltage period. The two regions above V * 1 (Ω) (respectively red and green in Fig. 6) correspond to these two cases. Analogously, applying a voltage with large amplitude, such that the memristor is conducting, and then gradually decreasing it, V * 2 (Ω) is the amplitude at which the memristor becomes insulating.
To discuss the frequency dependence of V * 1,2 (Ω), it is convenient to separately consider the regimes of low, intermediate, and high frequency.

Low frequency
At low frequency, the a. c. response to a voltage V (t) = V a cos(Ωt) is in a sense singular. On the one hand, at zero frequency the voltage reduces to constant. On the other hand, at non-zero albeit low frequency, it successively assumes all values in [−V a , V a ]. In other words, in the low-frequency limit the a. c. voltage is equivalent to an adiabatic sweep, such as considered in Sec. II. Thus, the steady state is insulating if V a < V * 1 and back and forth insulating and conducting if V a > V * 1 , as for repeated sweeps, cf. Fig. 3(b). Note the absence of conducting steady states in this limit, since no matter how large the amplitude, the memristor invariably turns insulating during the long interval in which the voltage assumes low values. This is reflected in the divergence of V * 2 (Ω), while V * 1 (Ω) is continuous and tends to the d. c. threshold V * 1 .

Intermediate frequency
The intermediate-frequency regime can be understood in terms of a competition of time scales: the half-period τ Ω = π/Ω; and the delay (τ D ) and relaxation (τ R ) times, namely the time scales for, respectively, the insulating-toconducting and the conducting-to-insulating transitions. While these were precisely defined in Sec. II for the d. c. transitions, here the discussion is more qualitative and depends only on τ D,R being, respectively, decreasing and increasing as a function of voltage amplitude.
Since within the range [V * 1 (Ω), V * 2 (Ω)] (green region in Fig. 6) there are one insulating-to-conducting and one conducting-to-insulating transition during each half a period [cf. Fig. 5(d),(h)], this region is characterized by the relation τ D , τ R < τ Ω . Indeed, if either time scale were longer than τ Ω , the corresponding transition could not take place. This suggests the interpretation of V * 1,2 (Ω) as the curves where, respectively, τ Ω = τ D and τ Ω = τ R . Crossing for example V * 1 (Ω), the region with insulating steady states is characterized by τ R < τ Ω < τ D , that is by the inhibition of the insulating-to-conducting transition. Within this perspective, V * 1 (Ω) increases with frequency because -as τ Ω decreases -a larger voltage amplitude is needed to match the condition τ Ω = τ D .
Where V * 1,2 (Ω) intersect each other, the time scales are all equal: τ Ω = τ D = τ R . Crossing this point at constant voltage amplitude, τ Ω becomes at the same time shorter than both τ D,R , meaning that both the insulting-to-conducting and the conducting-to-insulating transitions are inhibited, and the memristor remains in the same state as the initial condition (purple region in Fig. 6).

High frequency
The behavior at high frequency is better illustrated in terms of the infinite-frequency limit, in which the voltage is equivalent to a d. c. V a / √ 2. Indeed, the voltage enters the equation of motion [Eqs. (3), (11)] through the square [V a cos(Ωt)] 2 which at high frequency is equivalent to its average V 2 a /2. The steady state is therefore insulating if Similarly to the d. c. coexistence region, in the range [ √ 2V * 2 , √ 2V * 1 ] the steady state is insulating or conducting depending on the initial condition. Note that in the high-frequency limit V * 1,2 (Ω) tend to constant. To reconcile this with the previous discussion in terms of time scales, we have to consider that at high frequency the transitions can happen across multiple voltage periods.

IV. SELF-SUSTAINED OSCILLATIONS AND SPIKING BEHAVIOR
We study now a first use case of the Mott memristor in electric circuits. In the circuit in Fig. 7(a) the memristor is connected in parallel with a capacitor C and is attached to a voltage generator V through a load resistor R . This setup allows us to study self-sustained current oscillations as observed, e. g., in Refs. [33][34][35]; a phenomenon at the basis of spiking-based computational schemes.

A. Nullclines and fixed point
The equation for the voltage V across the memristor is obtained applying Kirchhoff's law of current conservation at the nodes of the circuit in Fig. 7(a): which are the currents through, respectively, capacitor, memristor and voltage generator. Equation (18) has to be solved together with the rate equation for the doublon density [Eqs.
(3), (11)]. Defining r = R /R 0 , r t = r s + r and the time scale τ c = R C we rewrite these equations as a dynamical system The fixed point of this system is at the intersection of the so-called nullclines, namely the curves along whicḣ V = 0 andṅ = 0, which read respectively V = V l (r s n + n 0 )(r t n + n 0 ) −1 , (20b) We subtract now the nullclines and, similarly to Sec. II, we solve the resulting equation for V , thereby expressing the fixed-point doublon density as the inverse function of Because the current through the capacitor is zero at the fixed point, the resistances R s and R are in series and the circuit reduces to the situation considered in Sec. II with the substitutions V → V and R s → R s + R , which indeed make Eq. (21) identical to Eq. (12). Therefore, the same analysis applies here: if r t > 0.125 the solution is unique, while if r t < 0.125 there is a region with three solutions, see Fig. 7(c). Together with r s = 0.01, we set hereafter r = 0.1 which gives r t = 0.11.
Depending on load voltage and load resistor, the fixed point can be in the NDR region of the Mott memristor, see Fig. 7(b)-(d), which is necessary for having limit-cycle self-sustained oscillations, as we discuss in the following.

B. Limit-cycle oscillations
Self-sustained oscillations are periodic solutions of a dynamical system, such as Eqs. (19), in absence of any periodic input. In the system configuration space (here the n-V plane) the corresponding trajectories are limit cycles, namely isolated closed trajectories which either attract or repel nearby ones [54]. Simply stated, the conditions for a limit cycle are the non-linearity of the system and the instability of its fixed point. In this case, the former is provided by the non-linear rate equation for the doublon density. The latter is satisfied if the fixed-point doublon density is between the values (see Appendix B) Once the load voltage -thus the fixed-point doublon density -is chosen, Eq. (22) gives an implicit expression for the critical τ * c , which varies with the fixed-point doublon density and whose minimum is obtained setting to zero the argument of the square root in Eq. (22): The region [ñ 1 ,ñ 2 ] is included in the NDR region of the memristor, coinciding with it in the limit of large τ c . As depicted in Fig. 7(d), to enter this region one can either tune the load voltage (thus the doublon density) or the capacitor (thus the characteristic time τ c ). At this point, a supercritical Hopf bifurcation takes place [54], namely the fixed point loses stability and a limit cycle arises.

Tuning the load voltage
At fixed τ c we consider three load voltages such that the fixed-point doublon density is below, inside, or above the unstable region [ñ 1 ,ñ 2 ], see Fig. 7(b)-(d). For each load voltage we numerically integrate Eqs. (19) with varying initial conditions and plot the trajectories in the n-V plane in Fig. 8(a)-(c). Outside the unstable region (V = 0.60, 0.90) all trajectories tend to the fixed point. Notice that this implies the absence of closed trajectories. In stark contrast, inside the unstable region (V = 0.75) there is an isolated closed trajectory (i. e. a limit cycle) which attracts all other trajectories. Notice that the limit cycle is around the unstable fixed point. In this case there is no stationary stable solution and, despite the constant load voltage, density and voltage oscillate indefinitely. In other words, the system undergoes limit-cycle self-sustained (or autonomous) oscillations.
The current profile is markedly different in the three cases. Let us consider [see Fig. 8(d)-(f)] the trajectories with initial condition (n 0 , 0). For V = 0.60 the current increases monotonically to the stable fixed point, which is on the insulating branch, see Fig. 8(d). In contrast, for V = 0.90 the stable fixed point is near the conducting branch and is reached only after a transient, which in the n-V plane takes the form of a spiral around the fixed point [ Fig. 8(c)], and the current profile has a single spike followed by damped oscillations, see Fig. 8(f). Finally, corresponding to the limit cycle, for V = 0.75 the current has periodic spiking, see Fig. 8(e). Each spike consists of a sudden increase and a similarly rapid, but slower, decrease. These are due to repeated transitions between the memristor insulating and conducting states, Note that this is consistent with the spiking behavior of biological neurons, in which the neural-cell membrane also transitions between insulating and conducting in the course of an oscillation [55].

Tuning the capacitor
In Fig. 9 we plot the trajectories obtained by numerical solution of Eqs. (19) with initial condition (n 0 , 0), fixed load voltage and varying τ c . With its location unaltered, the fixed point loses stability across a critical τ * c ≈ 6.06 τ d (for V = 0.75), see Fig. 7(d). It is stable for τ c < τ * c and reached after a number of oscillations which become more dense as τ * c is approached. As soon as τ c > τ * c , the fixed point becomes unstable and a small limit cycle appears. Increasing τ c further, the limit cycle grows and tends to a loop with segments at constant voltage connecting lower and upper branches of theṅ = 0 nullcline, see Fig. 9(d). Since this nullcline is nothing but the stationary doublon densityn versus the voltage [cf. Fig. 3(a)] this limit cycle is equivalent to the hysteresis loop in adiabatic voltage considered in Sec. II. In other words, in this limit the circuit behaves like a relaxation oscillator [54].
The limit-cycle current spikes can be characterized by height (difference between maximum and minimum) and period, see Fig. 10. Evidently, these quantities are only  defined for τ c ≥ τ * c . At τ * c we have the typical behavior for a supercritical Hopf bifurcation [54]: the height grows from zero (the limit cycle has vanishing amplitude) while the period is finite and equal to 2π[det(J)] −1/2 ≈ 12 τ d , with J the Jacobian of the dynamical system (19) at the fixed point (see Appendix B). Increasing τ c , the height first rapidly increases, then it slowly saturates to a value close to I * 3 which is, together with I * 1 , the stationary current at the threshold voltage V * 1 , cf. Fig. 3(b). At the same time, already at τ c ≈ 7 τ d , the period is linear in τ c , showing a decoupling of time scales for doublon density and voltage, as expected for a relaxation oscillator.

CONCLUSIONS
We have proposed the narrow-gap Mott insulator as a compact realization of a new type of memristor based on the field-induced carrier avalanche multiplication. Due to this purely electronic mechanism for the resistive switch, this Mott memristor has a characteristic time scale set by the doublon-excitation decay time τ d ∼ 1-10 ps, which is up to several orders of magnitude faster than in devices based on Joule heating or ionic drift.
As a first step we have put forward a phenomenological description of the field-induced carrier avalanche in Mott insulators, in which the conductivity depends on the carrier density, whose rate equation contains the non-linear scattering terms induced by strong correlations. Building on this, we have introduced the Mott memristor as a device made of a Mott material in series with a conventional resistor; and we have derived its current-voltage curve, as well as the transitions between conducting and insulating states. While the very definition qualifies the model as a non-polar, voltage-controlled memristive system, we have analyzed in detail its a. c. response, in particular the pinched hysteresis loop and the steady-state diagram as a function of amplitude and frequency. Finally, we have considered a circuit with a capacitor in parallel with the Mott memristor, and demonstrated self-sustained current oscillations and periodic spiking behavior, consistent with the periodic activity of biological neurons.
While similar devices have been subject of intensive experimental study [7,8,[32][33][34][35], this is the first time (to the best of our knowledge) they are proposed as memristors. Moreover, our work provides a comprehensive theory of the key features of those prior studies: threshold electric field, negative differential resistance (NDR), multivalued current-voltage characteristic, delay time, and current oscillations. At the same time, our proposal consists of a tractable set of equations; which stands in contrast with previous more complicated models, see e. g. Ref. [47], and results in two valuable features. First, it allowed us to derive analytical expressions, such as the boundaries of the NDR region and the conditions for limit-cycle oscillations. Second, and perhaps more importantly, it makes promising to include the model into the description of circuits of growing complexity, in the quest for bio-inspired novel computing architectures. real only if E < E th = 0.5 [cf. Eq. (5)]. Since during the delay time the doublon density does not change much, we approximate the field in the Mott insulator as constant, E ≈ V L, which is equivalent to approximating r s ≈ 0, yielding V * 1 ≈ 0.5V 0 and n * 1 ≈ 0.5n 0 . Equation (A1) can then be solved with a variable change: x +ẋ + E 2 x = 0.
Appendix B: Derivation of Eq. (22) In this appendix we derive Eq. (22) for the region with limit cycle in Sec. IV. A limit cycle is guaranteed to exist by the Poincaré-Bendixson theorem when the system is confined in a region with no stable fixed point therein [54]. Such a trapping region is (with V r t r −1 s > V * 1 ) {(n, V ) ∈ [0,n(V r t r −1 s )]×[0, V r t r −1 s ]}. The fixed point turns from stable to unstable (Hopf bifurcation) when, with positive determinant, the trace of the Jacobian becomes positive.

(B4)
The sign of the determinant does not depend on τ c and is positive for n outside the range [n 1 ,n 2 ] witĥ n 1,2 = n 0 (1 ± √ 1 − 8r t )(2r t ) −1 . The sign of the trace depends on τ c . Notice that a necessary condition for the trace to vanish is (−r s n 2 + n 0 n − 2n 2 0 ) > 0 which is the same condition for the NDR region of the memristor, cf. Eq. (13), demonstrating that the region with limit-cycle oscillations is a subset of the NDR region, as depicted in Fig. 7(b)-(d). Imposing the trace to be positive we get the condition that n should be outside the range [ñ 1 ,ñ 2 ] withñ 1,2 given in Eq. (22).