Stochastic behaviour of an interface-based memristive device

A large number of simulation models have been proposed over the years to mimic the electrical behaviour of memristive devices. The models are based either on sophisticated mathematical formulations that do not account for physical and chemical processes responsible for the actual switching dynamics or on multi-physical spatially resolved approaches that include the inherent stochastic behaviour of real-world memristive devices but are computationally very expensive. In contrast to the available models, we present a computationally inexpensive and robust spatially 1D model for simulating interface-type memristive devices. The model efficiently incorporates the stochastic behaviour observed in experiments and can be easily transferred to circuit simulation frameworks. The ion transport, responsible for the resistive switching behaviour, is modelled using the kinetic Cloud-In-a-Cell scheme. The calculated current-voltage characteristics obtained using the proposed model show excellent agreement with the experimental findings.


INTRODUCTION
For over 50 years, Moore's law has been the guiding principle of the electronics world. However, recent times have witnessed a breakdown of this law, with modern devices driving the future of computing beyond Moore's law [1]. Memristive devices (or memristors) are a kind of 'More-than-Moore' devices that have become a hot topic of interest for researchers in the field of solid-state physics, materials science, electronics and computer engineering [2,3]. In its generic form, a memristive device is a two-terminal device with an active (often electrically insulating) layer sandwiched between two metal electrodes. In simple terms, it is nothing but a nonlinear resistor with memory, in the sense that the Ohmic resistance depends on the history of the current that has flown through it driven by an applied voltage [4,5]. The switching from a low resistance state (LRS) to a high resistance state (HRS) and vice versa is initiated once the driving voltage reaches a specific value. With features such as low power consumption, passivity, and scalability down to the nanometer range, these devices have emerged out as potential candidates in the field of hardware security [6,7], non-volatile memory applications (e.g., resistive random access memories, RRAMs) [8][9][10], and neuromorphic systems [11][12][13][14].
State-of-the-art memristive devices based on fielddriven nano-ionic transport can be either interface-type or filamentary-type. The latter relies on the formation and dissolution of conducting filaments. These conducting filaments may consist of metal atoms or cations in so-called electrochemical metallization (ECM) cells [15][16][17][18] or oxygen vacancies in, e.g., hafnium oxide valence change mechanism (VCM) cells [19]. The switching behaviour of such filamentary devices is often digital and fast, and the states are stable. Therefore, the applications are predominantly in the field of RRAMs. In contrast, interface type devices are often driven by the movement of ions or charged defects within the active layer that modify the electrical characteristics of the metal/insulator interfaces (Schottky contacts or tunneling contacts), which act as boundaries for the active layer [20][21][22]. The switching dynamics is therefore analogue. The states are less stable, so the application of these devices is more in the field of neuro-inspired electronics. An important feature common in almost all nano-ionic memristive devices is that they show intrinsic stochastic behaviour due to the involved ion movement.
A number of different models have been proposed to explain the operation and switching dynamics of the devices on the experimental time scale of seconds to days. The actual switching mechanism on the atomic time scale is seldom addressed [23]. Although the models are able to mimic the electrical behaviour of the devices, they all have certain shortcomings. Concentrated models are based on complex mathematical formulations of the switching dynamics [24][25][26]. Despite being fast, the physical and chemical processes responsible for the actual switching dynamics are immensely neglected, and thus, they mostly show deterministic behaviour. Besides, more elaborate stochastic models are based on multi-physical 2D or 3D approaches that include more or less realistic physics. For example, 3D kinetic Monte-Carlo (kMC) based simulation models produce characteristics similar to real-world memristive devices, successfully capturing their intrinsic stochastic behaviour [19,[27][28][29]. Nevertheless, these models are computationally expensive, and therefore, it is not feasible to couple them with circuit simulators to address memristive systems. The drawbacks of the models mentioned above form the prime motivation for our work.
In the framework of the present work, a spatially onedimensional physics-inspired model that considers the stochastic nature of resistive switching is proposed. It is focused on the random behaviour observed in a realworld interface type device, showing a comparable cycleto-cycle and device-to-device variation while conserving the physics. The model is based on Newton's laws of motion for the ion motion in the active layer combined with the kinetic Cloud-In-a-Cell (CIC) scheme [30] (also known as the Particle-In-Cell scheme, which is often used in plasma physics) to couple the ion transport mechanism with the electric field within the device. A feature of this approach is that it is deterministic. However, the stochastic behaviour comes into play by randomly perturbed charge positions within the CIC scheme. It is important to notice that the proposed model is computationally less expensive than the usual kMC-based models. Therefore, it is compatible with SPICE-level or Verilog-A circuit simulations. The proposed stochastic model and its effectiveness are discussed in detail, taking the double barrier memristive device (DBMD) [20] as the subject device. To understand the proposed model, we start with the introduction of the DBMD device in Section II. The description of the simulation model is then given in Section III, including the implementation of the numerical scheme, stochasticity modelling and the current transport mechanisms. Finally, the results are discussed in Section IV, and conclusions are drawn in Section V.

II. DEVICE CONFIGURATION
The interface type double barrier memristive device, shown in Fig. 1, was introduced by Hansen et al. [20]. The device mainly consists of three interfaces in the order Au/Nb x O y /Al 2 O 3 /Nb. The Au/Nb x O y interface acts as a Schottky contact, which is assumed to be inert to prevent oxidation. The Al 2 O 3 tunnel barrier is an electrically high-quality barrier that allows for elastic electron tunnelling. The Nb x O y layer in the centre represents a solid-state electrolyte that acts both as an electronic and ionic conductor. The mobile negatively charged oxygen ions (blue circles) in this layer are primarily responsible for changing the resistance of the device, while the stationary positively charged oxygen ions (red circles) maintain charge neutrality. Under the influence of an externally applied electric field, the negatively charged oxygen ions drift within the Nb x O y layer (the charge arrangement in HRS and LRS is indicated in Fig. 1). This oxygen ion drift eventually influences the interface properties of both the Schottky-contact and the tunnel barrier, leading to an overall change in device resistance. While the field-induced ion movement can model the switching mechanism very well, as it is decribed in detail in [20,27], the charging and discharging of interface traps cannot be completely disregarded as switching mechanism [20]. It should be noted that the device stack shown in [20,27] differs slightly from the here presented one. It has been shown by electron energy loss spectroscopy, that the few nanometer thin Al bottom electrode is fully oxidized and even the Nb layer underneath is slightly oxidized to a sub-stoichiometric NbO z at the interface [31]. However, this does neither influence the model shown in [20,27] nor the 1D model presented here, because the NbO z is believed to be a metallic resistor and tunneling through the alumina is still reasonable [31]. The device parameters used in the 1D model in this work are extracted from [20] and are summarized in Table I.

III. MODEL DESCRIPTION
The equivalent circuit of DBMD is used as the starting point for the electrical representation of the compact model. In contrast to the equivalent circuit of DBMD by Hansen et al. [20], a simplified version of the equivalent circuit is adopted. The equivalent circuit diagram of the The equivalent circuit model of a DBMD with a Schottky diode, solid-state electrolyte as an Ohmic region and the tunnel barrier as a voltage-controlled current source. Fig. 2. The Schottky contact is represented as a Schottky diode and the tunnel barrier as a voltage-controlled current source. The ion distribution in the solid-state electrolyte is expressed as a resistance (R SE ). By applying Kirchhoff's voltage law (KVL) and Kirchhoff's current law (KCL) to the equivalent circuit, we obtain the following equations,

DBMD is shown in
V Device , V SC , V SE , V TB are the voltage drops across the device, Schottky contact, solid-state electrolyte and tunnel barrier, respectively. Similarly, I SC , I SE , I TB are the currents across the Schottky barrier, solid-state electrolyte and tunnel barrier, respectively. Based on the equivalent circuit diagram and Eqs. (1)-(2), the ion transport and the current mechanisms are now implemented and discussed in the following subsections.

A. CIC inspired ion transport
As mentioned in Section II, the overall change in device resistance is mainly due to ion transport. Therefore, it is necessary to choose a suitable transport mechanism that produces characteristics equivalent to a physical device. Unlike the available ion transport models, such as the kMC model by Dirkmann et al. [27], a much simpler transport method inspired by the Particle-In-Cell (PIC) scheme [32,33] is employed in the present work. This method is well known as the Cloud-In-a-Cell (CIC) scheme in the semiconductor community [30] .
The computational domain consists of a solid-state electrolyte (l SE ) modelled on a 1D Eulerian grid. Initially, an equal number of mobile and fixed ions are randomly distributed within the solid-state electrolyte considering a Lagrangian grid. The spatial arrangement of ions in this form is regarded as the HRS depicted in Fig. 1. For simplicity, from here, we consider two loops, (a) the inner loop for solving Newton's laws and (b) the outer loop for ion transport.
Inner loop: A basic mathematical optimisation is performed to obtain V TB by minimising the error between I TB and I SC . For this, the initial value of V TB is assumed to be V Device . Then, the other circuit parameters (mentioned in Fig. 2) are calculated and reiterated until Eq. (2) is satisfied. At the end of each loop, V TB is modified according to the relative error (ξ) between I TB and I SC . Finally, using the voltage drops, the potentials at the Au/Nb x O y interface and Nb x O y /Al 2 O 3 interface are respectively given by and Outer loop: Followed by the inner loop, the electric potential is calculated using the 1D Poisson's equation, where ε is the permittivity of Nb x O y and ρ is the charge density. The potential equation is subjected to Dirichlet boundary conditions at the simulation domain boundaries, given by Eqs. (3) and (4). The electric field is then obtained from the electric potential by This electric field E is used to push the mobile ions with a certain drift velocity assuming the Lagrangian grid. The spatial arrangement of drifted ions in LRS is illustrated in Fig. 1. In solid-state physics, ion transport on an atomic level is characterised by jump attempts over a potential barrier. Based on this theory, the drift velocity of an ion is given by [34,35] v D = dp( where d is the jump distance (lattice constant), z is the charge number of the ion, k B is the Boltzmann constant, T is the temperature, and e is the elementary charge. The transition probability, p(E A ), for an ion to jump to the neighbouring site is where n c = 0.5 exp −∆Hg 2kBT is the probability that a site is occupied, then (1 − n c ) gives the probability that the next neighbour site is empty. ν 0 is the phonon frequency, ∆H g is the intrinsic energy gap, E A is the activation energy, N is the number of these empty sites, and f is a geometrical factor of order unity. For a 1D transport, the term N (1 − n c )f is approximately equal to unity. It must be taken care that the Lagrangian grid is only assumed for distributing the ions and calculating drift velocity, and the rest of the simulation is performed on an Eulerian grid. As further detailed in the subsequent sections, the mobile ion movement consistently changes the internal state of the device. The average relative distance between the current position of mobile ions and their inertial position determines the internal state in the proposed model. The internal state is given by Here,d is the average distance between the mobile ions and Schottky contact, andd r is the initial average distance.x i is the position of i th mobile ion,x SC is the position of Schottky contact and N ions is the number of mobile ions. The flowchart of the overall implementation of the simulation approach is shown in Fig. 3.

B. Stochasticity modelling
An essential property of the considered oxide-based DBMD is its intrinsic stochastic nature. The stochastic behaviour is usually measured in terms of cycle-tocycle (C2C) variability and device-to-device (D2D) variability. It has been experimentally observed for a DBMD as shown in Figs. 6(a) and 7(a). So far, the stochasticity in the model is only induced by random initial positioning of the mobile and fixed ions in the solid-state electrolyte. However, this is insufficient to replicate the stochastic behaviour observed in a physical device since it only mimics the D2D variability. As the drift theorybased ion transport following Eq. (7) is deterministic, the internal system state q(t) also tends to become deterministic. However, q(t) changes randomly as the mobile ions move randomly across the solid-state electrolyte for a physical device. So to account for stochasticity, we artificially perturb the internal state q(t) of the device [36]. Notably, the internal state is changed indirectly by randomly perturbing the position of the mobile ions at each iteration of ion motion. So the position of the i th ion, x i(s) is given bȳ where r is a uniform (pseudo) random number between 0 and 1, and δ denotes the maximum random displacement. The latter is chosen so that the underlying physical processes do not change and remain stable (i.e., 1% − 5% ofx i(d) ). By modifying the position of the mobile ions, the average distanced also changes, which in turn makes q(t) stochastic. Thus, the C2C variability feature of the memristive devices is mimicked.

C. Implementation of current mechanisms
The calculation of current through the different regions of the DBMD is explained in the following sections:

Tunnelling current
The electronic current flow through a metal-insulatormetal (MIM) system can be given by the electric tunnel effect. The tunnelling phenomenon of electrons through the Al 2 O 3 tunnel barrier can be expressed using a set of equations derived by Simmons [37]. The general Simmons current equation is given as where Φ t eff is the effective tunnel barrier height, d eff is the effective tunnel barrier width, A d is the device area, and β is a correction factor. e, m, and h are the elementary charge, the free electron mass, and the Planck constant, respectively. The absolute value of V TB is considered here, assuming that Eq. (12) is symmetric between positive and negative voltage biases. Since the Simmons formula is only considered for elastic tunnelling, the resulting electronic current mainly depends on the local ion concentration at the Al 2 O 3 /Nb x O y contact. So, the effective tunnel barrier width and the effective tunnel barrier height are used here, which are given by and Φ t eff = Φ t (1 + λ t q(t)). (14) d 0 is the actual width of the tunnel barrier, Φ t is the actual tunnel barrier height. λ d and λ t are the fitting parameters. The fitting parameters are chosen in a way that the simulation outcome leads to specific experimental behaviour, and their values can be between 0 and 1.

Ohmic current
The current across the Nb x O y solid-state electrolyte (I SE ) can be calculated directly using Ohm's law as, where σ is the conductivity of Nb x O y and l SE is the length of Nb x O y . Also, from Eq. (2), I SE = I TB , so

Schottky current
In the presence of a high electric field, the oxygen ions move close to the Nb x O y /Au Schottky interface and back to their inertial position. The concentration of these ions at the Schottky contact considerably affects the Schottky barrier height and the ideality factor. The current conduction mechanism across this Schottky contact is described by thermionic emission theory [38]. According to this theory, the Schottky diode current is given by where n eff is the effective ideality factor which describes the deviation from an ideal current, k B is the Boltzmann constant, and T is the temperature. The reverse current, I R for different voltage polarities, is given by forward bias: reverse bias: eff is the effective Schottky barrier height, and A * = 1.20173 × 10 6 A/(m 2 K 2 ) is the effective Richardson constant. α r is a device-dependent parameter that describes the voltage dependency of the reverse current as measured experimentally. The reverse current during the set process is dominated by lowering the Schottky barrier height, decreasing gradually with the applied positive bias. While the Schottky barrier height increases during the reset process. By considering the spatial rearrangement of mobile ions due to the high electric field, the effective Schottky barrier height and effective ideality factor are defined as n eff = n 0 (1 + λ n q(t)).
Φ b0 is the initial Schottky barrier height, and n 0 is the initial ideality factor.λ b and λ n are the fitting parameters.

IV. RESULTS AND DISCUSSION
The current-voltage characteristics (I -V curve) of the double barrier memristive device calculated using the proposed stochastic model is compared to the I -V curve measured experimentally [20] and the I -V curve calculated using the kMC simulations [27] in Figs. 4 for a single cycle. For the simulation, a linear voltage sweep was applied from 0V to 3V and back to 0V, in order to set the device from HRS to LRS. To reset the device to its initial HRS, the voltage was ramped linearly from 3V to -1.5V and then back to 0V. Then the resultant I -V curve, shown in Fig. 4(a), resembles a typical memristive hysteresis loop. As observed, the I -V curve extracted from the 1D stochastic model (shown in Fig. 4(b)) is qualitatively in good agreement with the measured I -V curve and also the kMC simulated I -V curve, but quantitatively the three I -V curves show slight variation. How-ever, this minute variation is acceptable considering that the device and the simulation models are stochastic.
To investigate the movement of mobile ions in the solid electrolyte under the influence of an applied voltage and their influence on the device resistance, the Spatiotemporal plot of ion transport for a single applied voltage cycle is shown in Fig. 4(c). The coloured lines in the plot represent 20 different ion movements that are distributed randomly across the solid-state electrolyte. The absolute average position of the ions,d displayed as a black dashed line in Fig. 4(c), indicates the internal state of the device. Additionally, the four significant positions highlighted in red correspond to analogous points on the I -V curve in Fig. 4(b). Initially, at point (i), all the ions are randomly scattered across the simulation domain (corresponding tō d = 1.1 nm). When a positive voltage is applied, the ions start to drift towards the Schottky interface (indicated as point (ii)), thus changing the device state from HRS to LRS. For a negative voltage ramp, the ions drift in the opposite direction, and the value ofd is then close to the initial value at -1.5V (point (iii)). The voltage is then increased from -1.5V to 0V to reset the device completely (point (iv)).
Since the process is stochastic, the value ofd is variable but should be within a range of approximately ±δ of d (deterministic). This can be evaluated using the phasespace plot ofd and the derivative ofd with respect to the time, as shown in Fig. 5. The plot is obtained for the movement of ions recorded for ten consecutive applied voltage cycles (each colour corresponds to a single cycle of the applied voltage). It is a kind of Poincaré plot that describes the periodic behaviour of the model. Although the trajectories ofd in the phase-space plot are distinct and random for all applied voltage cycles, their ensemble is still within an arbitrary range of ±5% ofd (black ellipse shown in Fig. 5), hence proving the periodicity and stochasticity of the proposed model. The range determined by the black ellipse depends on the value of δ, which is taken as 5% in the manuscript.
The stochasticity of the proposed model can be further determined by observing the I -V curves in Fig. 6 and 7. The plots show the simulated and measured C2C variability and D2D variability of DBMD obtained under comparable conditions. Firstly, to check for the C2C variability, the simulation was performed for four consecutive cycles of the applied voltage. As a result, four different I -V curves are observed in Fig. 6(b), which follow the same trend yet have slightly varying R on /R off ratios. The deviation from one cycle to another can be reduced or increased by considering a smaller or larger value for δ, respectively. Similarly, to verify the D2D variability, the simulation is performed four times with different initial ion arrangements. The difference in the ion arrangements here corresponds to having four different real-world DB-MDs. The resultant I -V curves in Fig. 7(b) are distinct from each other with different R on /R off ratios. Thus proving the stochasticity of the proposed model, which very well reproduces the inherent stochastic behaviour of real-world DBMDs shown in Figs. 6(a) and 7(a).
In comparison to the state of the art simulation models, the proposed model shows a significant C2C and D2D variation, similar to experimental results. This makes the proposed model well suited for performing circuit simulations in artificial intelligence (AI) computing, reconfigurable logic computing, or hardware security primitives, which mainly rely on the stochastic behaviour of memris- tive devices. Moreover, the set current is nearly the same for experimental and calculated I -V curves in both C2C and D2D plots (see Figs. 6 and 7). This is because, during the set process (LRS), most of the ions are present at the Au interface, regardless of their initial position and applied voltage cycle. However, due to their random movement, the ions do not migrate back to their original position during the reset process. In other words, the value ofd for LRS is relatively similar for all cycles, butd is always different for HRS. Therefore, the reset current, which depends ond, is always different. This is clearly noticeable in both experimental and calculated I -V curves.
The applied voltage used to obtain the plots in Fig. 4 follows an unsymmetrical triangular waveform, particularly a sawtooth waveform. The sawtooth waveform was chosen solely to study the set and reset process in a device. Using the proposed 1D stochastic model, the response of the device to other input signals with different set and reset voltages can also be studied. Fig. 8 shows the I -V curve of a DBMD obtained using the proposed model for a sinusoidal applied voltage. The I -V curves are obtained for four different maximum applied voltages (2V, 2.5V, 3V and 5V). As expected, the I -V curves show smooth transitions at the maximum positive and negative bias. However, for the positive branch of the hysteresis, asymptotic behaviour is observed for voltages more than 3V. Since the set voltage for DBMD is 3V (as mentioned in [20]), the device does not show any significant change for voltages more than 3V. Moreover, it is also verified experimentally that the width of the hysteresis increases as the maximum applied voltage increases, whereas hysteresis effects collapse at lower voltages [27,39]. The physical devices suffer a dielectric breakdown at high voltages around ± 5V and even for lower absolute values of negative voltages.

V. CONCLUSION
A circuit simulator compatible 1D stochastic model for simulating interface type memristive devices (such as DBMD, BFO) is proposed in this paper. The model integrates the stochastic behaviour of a real-world DBMD with a realistic ion transport by consistently coupling Newton's laws with a Poisson solver, inspired by the Cloud-In-a-Cell scheme. The DBMD current-voltage characteristics obtained using the proposed model are in good agreement with experimental results. Furthermore, different hysteresis relations are obtained while performing several simulations, describing the stochastic nature of the model. Despite being a straightforward 1D model with minimum mathematical formulations, the proposed model successfully reproduced all the essential features of a real-world DBMD. These features include the nonlinear analogue behaviour, intrinsic stochastic behaviour, an asymmetry between positive and negative bias, and high voltage current saturation. Therefore, the proposed stochastic model is perceived as a powerful tool to simu-late large-scale memristive circuits in areas such as bioinspired neural networks, hardware security or electrical networks. Additionally, the model could be used to study the dynamics of other non-filamentary memristive devices with minimal computational resources.