Insight into physics-based RRAM models – review

: This article presents a review of physical, analytical, and compact models for oxide-based RRAM devices. An analysis of how the electrical, physical, and thermal parameters affect resistive switching and the different current conduction mechanisms that exist in the models is performed. Two different physical mechanisms that drive resistive switching; drift diffusion and redox which are widely adopted in models are studied. As for the current conduction mechanisms adopted in the models, Schottky and generalised hopping mechanisms are investigated. It is shown that resistive switching is strongly influenced by the electric field and temperature, while the current conduction is weakly dependent on the temperature. The resistive switching and current conduction mechanisms in RRAMs are highly dependent on the geometry of the conductive filament (CF). 2D and 3D models which incorporate the rupture/formation of the CF together with the variation of the filament radius present accurate resistive switching behaviour.


Introduction
Scaling limitations of conventional non-volatile memory lead to the invention of emerging technologies like phase change memory (PCM) [1], resistive random-access memory (RRAM), and spin transfer torque random-access memory (STTRAM). Amongst these devices, RRAMs have demonstrated potential for adoption as next generation solid-state non-volatile memory devices because of its simple structure [2], compatibility with CMOS [3], high endurance, retention property and nanoscale dimension [4].
RRAMs are two-terminal devices with a top and bottom electrode sandwiching a transition metal oxide (TMO) layer in a metal-insulator-metal (MIM) structure. The TMO layer is capable of switching between stable oxidation states with different resistances, allowing data storage in the form of a high '1' and low '0'. Typically, RRAMs require a forming step to form a conduction path in the TMO layer which shunts the two electrodes. The existence of this conductive filament (CF) lowers the resistance of the RRAM and can be ruptured to raise the resistance. The rupture and formation of the CF is repeatable and the RRAM is now able to switch between a low resistance state (LRS) and high resistance state (HRS).
Further studies on engineering the material and structure of RRAMs to fine-tune its performance have been performed [5,6] and RRAMs can now be further categorised by its structure and the resistive switching characteristics. RRAMs can be divided into two categories, oxide-based RRAMs (OxRAM) [7] and conductive bridge RRAMs (CBRAM) [8,9] depending on the mechanism behind their RS behaviour. RS in OxRAMs stems from the driftdiffusion of oxygen vacancies (V O ), while the migration of metal ions contributes to RS in CBRAMs [10,11]. As the general RS mechanisms and the models for both are different, this paper will focus on the modelling of OxRAMs. Any mention of RRAMs will refer to OxRAMs here from now onwards.
To facilitate studies on the RRAM behaviour in electronic circuits, different types of models have been developed namely electrical, analytical, mathematical, physics-based models and more general compact models [12]. Among these models, physicsbased models explain device behaviour with physical concepts, which is more acceptable to define the working of the device [13][14][15][16].
This paper is organised as follows, the second section discusses the general behaviour of the RRAM, the third section describes the RS process and the modelling parameters used in different approaches, the fourth section discusses the current conduction mechanisms in recent models and the parameters utilised and the paper is concluded with the analysis and findings from the discussions.

SET and RESET switching behaviour
The switching behaviour of the devices during state transitions from LRS to HRS (RESET) and from HRS to LRS (SET) can be either abrupt or gradual depending on the electric field strength, internal temperature, thermal conductivity, hopping distance, and various other properties of the RRAM material. A high electric field across the filament gap length (g) defines the abrupt SET and a low field across g causes gradual RESET in [5,13], whereas temperature influences the rate of oxygen vacancy formation/ recombination to control the changes in CF radii and gap length in [17]. Similarly, thermal feedback mechanisms affect SET and RESET switching in [18,19]. Positive feedback loops are generated from the self-acceleration of ions to cause a sharp SET process by enhancing the electric field and temperature of the filament tip, while negative feedback is generated from the accelerated ionic migration and subsequent creation of the filamentary gap to form the basis of gradual RESET. The physical mechanism reported in [20] is the drift-diffusion of V O s in the filament with changes in field and temperature. In this case, gradual SET switching is a result of the electric field and temperature-driven drift-diffusion processes acting in same direction, whereas the opposing drift and diffusion fluxes produce a cancelling effect to account for gradual RESET. The model presented in [15] reports the abrupt SET and RESET switching governed by the change in the volume of oxygen vacancy doped region and its corresponding variation with Schottky barrier.
In addition to the electric field and temperature, the material properties such as hopping distance and thermal conductance are vital in predicting the abrupt and gradual changes in the switching behaviour. A smaller hopping distance yields analogue switching characteristics, while a larger hopping distance causes abrupt transition during switching. RRAMs with high thermal conductivity electrodes tend to exhibit gradual switching behaviour by drawing more heat and thereby decreasing the temperature to impact the drift and diffusion.

Onset of origin of Gap opening in conductive filament
The modelling of RRAM devices by adopting a physical approach provides an insight on the rupture of CF during RS. The prominent factors which affect the site of CF rupture are; the location of oxygen/metallic ion reservoir [21], temperature profile along the CF [18] and thermal properties of oxide and electrode.
The region near the oxygen/metallic ion reservoir layer provides favourable conditions for the growth and destruction of filaments as there are better chances of oxygen vacancy recombination and release or migration of ions. For instance, filament rupture starts near top electrode in [22] and at interface between the switching and bulk layer in bi-layered structures [21,23].
In temperature-driven models, the location of the depletion gap starts at the point of highest temperature in the CF. The CF initially ruptures at middle of the CF is reported in [18] and at top electrode interface is given in [13,20]. The onset of gap can be determined by selecting the suitable materials for electrode and oxides. Since the metal electrodes are great heat conductors, they act as heat sinks. The thermal conductivity of both top and bottom electrodes plays a role in the CF rupture site formation in single-layered devices, whereas the thermal properties of the electrode next to switching layer determines the location of gap formation in bilayered devices. Therefore, to determine the change in temperature during RS, the heat flow in CF is calculated by taking the ratio of thermal conductivity of electrode and oxide [20,24]. Moreover, in a few models, the Schottky contact formed at metal-semiconductor junction facilitates as an interface for the onset of CF rupture [5,25,26].
From the analyses above, it can be understood that CF formation or rupture is modelled based on two major physical mechanisms namely: drift-diffusion and reduction-oxidation (redox). Both mechanisms are driven by electric field or temperature or a combination of both. The following sections analyse the physics behind the two approaches in detail as well as the physical concepts adopted to model the current conduction.
3 Resistive switching mechanisms adopted in analytical models 3.1 Approaches adopted to model resistive switching 3.1.1 Drift-diffusion model: In this approach, RS is modelled through defect migration induced by an electric field and an increase in the localised temperature during SET and RESET. The application of an electric field drives the migration of ions with migration direction dependent on the ionic charge and field strength direction. The ionic migration is a result of the interaction between two fluxes: Fick's diffusion and ionic drift. Fick's diffusion is dependent on the concentration gradient while electric field-induced barrier lowering influences ionic drift. These fluxes contribute to the change in oxygen vacancy (V O ) concentration in the CF which is given by [19]: where D∇n D and μFn D are the Fick's diffusion and ionic drift terms, respectively. n D is the doping concentration, D is ion diffusivity, is the ion mobility, and F = −∇ψ where ψ is the local electrostatic potential demonstrating the relationship between drift and electric field. D and are based on Mott-Gurney ion hopping law. The equation above can be solved with the carrier continuity (2) for electronic conduction to incorporate the temperature effect [19]: where σ is the electrical conductivity, and with the steady-state Fourier equation [19]: where k th is the thermal conductivity and T is the local temperature.
In [20],a fitting parameter γ is added to the rightmost term to become γ ⋅ σ ∇ψ 2 where γ = 1 and γ = 2 are used for DC and AC programming conditions, respectively, for the simulation. Solving these expressions provides the change of n D , ψ, and T during RS. Additional physical terms included in (1) are incorporated in [20,27] to fine-tune the model. In [27], the migration or generation of V τ through thermally activated hopping during SET is considered termed G: As a result, the growth of the CF is considered proportional to the rate of ion migration and the generation term, G is given by [27]: where A is the pre-exponential constant, E b is the energy barrier for ion hopping and K is the Boltzmann constant. The qβE term is the barrier lowering effect of the applied electric field, E where β is the simulation's mesh size. The addition of the G term is shown in (E2) of Table 1.
The G term in (5) describes V O generation during SET and is not applied during the RESET phase of the model.
A temperature influenced Soret's effect or thermophoresis term is added to (1) to account for the contribution from Joule heating in [20]. Thermophoresis which describes migration of vacancies along the temperature gradient is given as [20]: where S is the Soret coefficient. This term describes the tendency for V O to move towards a hotter region [20]. The revised form of (1) is given as (E3) in Table 1.
The change in n D during RS affects the electrical and thermal conductivity, and k th as well as the activation energy for conduction, E AC through the thermally activated electrical conductance Arrhenius equation: where 0 is the pre-exponential factor, E AC is the activation energy for conduction, and k is the thermal conductivity. Both σ 0 and E AC are dependent on n D as the increase in V O increases the device conductance and reduces the activation energy. The change in σ 0 and E AC during RS are plotted in Fig. 1. The PDEs from (1)-(3) should be solved for all layers including the electrode layers to account for the contribution of the entire device to RS. In [19], the electrode layers act as ideal heat sinks with boundary condition T = T 0 = 300 K. This is a reasonable assumption as the electrode area is generally large compared to the CF. Since [19,20,27]modelled a RRAM with inert electrodes, they do not contribute to ion migration and, therefore, have constant k th and no V O flux at the top electrode (TE)/oxide and oxide/ bottom electrode (BE) interfaces. In [19,20], bi-layered RRAMs are modelled and the PDE is likewise solved for the bulk layer. As the bulk layer acts as an oxygen vacancy reservoir, only Fick's diffusion was considered and the drift term, μFn D in (1) is zero.
The I-V curves generated from this approach in [19,20,27] are shown in Fig. 2.

3.1.2
Reduction-oxidation (redox)-based modelling approach: In this approach, redox mechanisms are applied to model defect evolution during CF formation and rupture. The redox mechanism is a result of the generation and recombination (G-R) of V O with oxygen ions, oxide phase transition (P-T), as well as the hopping and migration of generated defects. The SET process comprises the P-T from semiconductor to metal phase and the generation of oxygen vacancies [21,22], while the RESET process involves the reverse P-T from metal to semiconductor phase and the release and hopping of oxygen ions (12) into the switching layer to recombine with V O [24,29]. The variation of the filament gap is given as: where g is the gap length, α is the field enhancement factor, and f is the attempt-to-escape frequency. When a gap in the filament is present, the activation energy for conduction, E A is lowered due to the Poole-Frenkel effect by the electric field in the gap given as: where α 0 is the activation energy barrier lowering factor, E A0 is the barrier at zero field, and V gap is the voltage across the gap. In addition to the field-induced evolution of the filamentary gap, the redox models in [18,22,24] consider the change in filamentary radius during RS. The I-V curves and process schematic from [22,24] are shown in Fig. 3a and 4a, and Fig. 3b, and Fig. 4b, respectively. The radii evolution is calculated as: where r is the CF radius and Z is the charge of ion/vacancy. The I-V curve generated is shown in Fig. 3a.
Another physical mechanism governing RS is the Butler-Volmer redox reaction which considers the role of oxidation and reduction on material property. This reaction is given as [28]:  [19] showing different SET switching voltages as a result of different RESET stop voltages, (b) [20] showing the effect of different EA values, and (c) [27]

Fig. 3 I-V curve and filament growth model of Hfτ x -based RRAM
(a) Simulated I-V curve, (b) Schematic of CF growth from [22] 4646 J. Eng where r and r max are the changing radius and maximum radius of the CF, respectively, redox is the nominal redox rate, and V cell is the voltage across the device. The RS process can be modelled with the same defect kinetics for both SET and RESET by considering the rate of migration of ionised defects from one part of filament to another with the Arrhenius rate expression (8) and (12) [18,29]. Gap length variations due to migration can be modelled with (8) and the radial rate of growth can be implemented with (10) or (11). The movement of oxygen ions and G-R of vacancies can be included to form a single expression: to model SET and RESET [29]. Similarly, the thin film growth rate is described with the concept of oxidation in the Mott-Gurney model in [5] which is utilised to model the CF gap length modulations during RS process. From the expression, γ can be considered as a constant or variable [5,29] depending on the models. Since CF geometry determines the state during RS, the filament growth can be modelled three dimensionally by considering the volume of the CF [17] which takes the form (13) for SET process while RESET (14) is given as: where α a and α g are the electric field enhancement factors.

Activation energy:
The list of activation energies used in different models that differ according to the switching material are given in Table 1. It is found that the activation energy required for redox mechanism models are much higher than drift-diffusion models.

Contribution of thermal conductance of electrode and oxide layers in RRAM:
Thermal conductivity is a material property which captures the change of defect concentration in the gap length [18]. The change in vacancy concentration in CF affects its thermal conductivity whereby the thermal conductivity at minimum vacancy concentration takes the value of the insulating oxide (e.g. HfO 2 , Ta 2 O 5 ) which linearly increases to the conductivity of the metallic CF (Hf, Ta) at maximum vacancy concentration. The enhancement of the thermal conductivity from the increase in vacancy concentration is a result of the free-carrier contribution to heat conduction. To solve (3), the change in thermal conductivity is modelled using Wiedemann-Franz law as [19]: where k HfO/Ta2O5 is the thermal conductivity of the oxide at T = T 0 = 300K and is the linear thermal coefficient. The values for k th differ according to the metal. The linear evolution of k th is shown in Fig. 2.
In addition to thermal conductivity, resistivity of oxide material and CF forming material has an influence on the device characteristics. Resistivity of the oxide layer in the gap region directly affects Joule heating and influences V SET and thus R on /R off ratio [30]. With an increase in resistivity of oxide layer, the device can achieve increase in V SET and decrease in the leakage current. These two electrical parameters significantly enhance device functionality and efficiency. Likewise, resistivity of base layer and electrodes also contribute to RS. Resistivity of base layer has direct dependency on providing self-compliance in bi-layer devices [5,21] during SET process by facilitating a large resistance series to low resistance of CF. In [17], the resistance of electrodes was considered in addition to the resistance of switching and bulk layers (Fig. 5).

Barrier height reduction factor:
The lowering of potential barrier due to the application of electric field is shown to be an important parameter in RS approaches. The term α a ZeE denotes general form of barrier reduction term. [18,22,24,29] where α a(/h) , is the field enhancement factor which can be a constant or variable quantity based on the models [18,22]. A few models formulated variable enhancement factor for barrier lowering [17,29]. α a is function of g in [17], whereas in [29], enhancement factor (γ) which depends on polarisability of the material [29]. Z is the charge number of oxygen ions/vacancies [22].

CF geometry-dependent physical parameters:
State of RRAM devices is derived from RS process which is demarcated with variations in CF geometry as well as shape. Important physical parameters which facilitate RS are, namely, the gap length (g) between top electrode (TE) and CF tip, the width (w) or radius (r) of CF and the rate of change of g or w, (dg/dt or dw/dt). The temporal evolution of g [27] is the common state parameter in 1D models while the change in filamentary r and/or w is the state parameter adopted in 2D and 3D models [22,24,28]. The volume of the CF is also considered in the 3D model in [17]. The SET  process in many models depends only on the gap length [29], while in a few models, both radius (10) and conduction filament length are taken into account [22,24]. The RESET process in most models are strongly dependent on g [18,22,24,29].
The shape of conduction filament is another factor to be considered while RS modelling. RS in the filamentary switching models in [18] is obtained by adopting a cylindrical-shaped CF, while a conical-shaped CF with variable upper and lower radius is considered in [17].

Temperature:
RS in many models apply Fourier heat flow equations and its developed forms [17,21,22] to calculate the effect of temperature during resistive switching. Major parameters required to implement these processes are thermal conductivity and specific heat per unit volume of the oxide material. In [21], the local temperature for the Pt/Ta 2 O 5 /TaO x /Ta RRAM is calculated from where C represents the specific heat per unit volume of the Ta 2 O 5 layer, k is the thermal conductivity of Ta 2 O 5 (9.6 W K −1 m −1 ), and Q is the Joule heat power density. Fig. 6 shows the simulated HRS and LRS transitions where A-H corresponds to the points in I-V curve. The structure layout is from left to right: top electrode, Ta 2 O 5 , TaO x (not shown), bottom electrode (not shown). The TaO 2 layer exists as an intermediary layer between the Ta 2 O 5 and TaO x interfaces due to the migration of oxygen ions (Fig. 7). The RRAM is initially in a LRS state at A. Reset begins at B when the supplied voltage is increased and the temperature can be seen to raise along the entire CF. This induces V τ migration towards the bottom electrode and the filament ruptures at C and the gap region increasing in D. The high temperature is now focused on the gap region as the electric field is concentrated there. The CF gap prevents high temperatures during SET (F) but enough thermal effect is present to induce V O migration and the reformation of the CF. The temperature profile in G is similar to B from the current flowing through the now connected CF.
In [24], the localised temperature effect induces both vertical and lateral migration of V O s where the Joule heating effect on the local CF temperature is calculated as where T 0 is the ambient temperature and R th is the thermal resistance of the CF. During SET switching, the CF in this model grows vertically to connect the TE with the BE and laterally once the connection is formed (Figs. 3b and 4b).
The effect of temperature on RS is studied in [17] where a model with fixed temperature (T = 300 K) is plotted against a model dependent on the heat equation [31] (Fig. 8). No RS occurs with fixed CF temperatures.
A simplified heat (21) is used to model the filament temperature variation in [5,24,29], while a 1D steady-state Fourier equation given as where z is the space coordinate along the CF (z = 0 at the injecting electrode and t ox = 20 nm is the oxide thickness and total CF length), k th is the thermal conductivity, ρ is the resistivity, and J is the current density is utilised to obtain the temperature variation in [18].
To obtain the temperature along CF, solution of 1D steady-state Fourier equation is utilised in [18].

Schottky-based current conduction:
Several models followed Schottky current conduction for modelling the interface between metal-semiconductor layers [5,25,26]. Schottky current expression takes general form as [14]:  In the given expression, A is area of the electrode, η is the ideality factor, k is the Boltzmann's constant, T is the temperature, V sch is the voltage across the Schottky junction, A * is the Richardson constant. exp − ξδ is the tunnelling probability parameter which includes the effect of tunnelling [32], ξ is the effective barrier in eV and δ is the interface layer thickness, ϕ Bi is the effective barrier height which depends on barrier height modulation due to oxygen vacancy doping effect ϕ Bn0 and barrier lowering due to image force and the electric field ϕ del given as: Different approaches adopted in modelling the current conduction are Richardson-Schottky and Simmon's methods [32][33][34].
Richardson-Schottky conduction is state dependent and is well suited to models with the mean free path of electron comparable to the thickness of the oxide film [33].

Generalised hopping conduction:
To correlate the gap length with applied voltage, a generalised hopping current conduction is used in conduction models [17,22,24,29]: where g T and V T are the characteristic length and voltage, g is the hopping/tunnelling length, V is the applied voltage across the cell.

Parameters required to model current conduction mechanisms 4.2.1 Barrier height:
The barrier height is the most critical electrical parameter in the thermionic emission process [34] which is determined by the physical properties of the oxide layer, structure of RRAM devices, image force lowering effect and electric field. Effective barrier height can be calculated using (19) by taking the effects of image force lowering and doping (Table 2).

Image force lowering effect with uniform electric field within the depletion layer:
In the presence of an electric field, the opposite charges accumulated at the interface of the metaldielectric have an effect of reducing the energy required for the carriers to overcome the barrier and is termed as image-forcelowering. Different methods are adopted to implement this concept as shown in Table 2, under the category of image force lowering [13,15,25,26]. These models derive parameters based on different concepts. For example, the permittivity of the semiconductor layer ε S is obtained from the product of optical dielectric constant ε ∞ and insulator dielectric ε ins in [15].

Temperature:
As a result of the strong current flow and electric field present in RRAMs, Joule heating is responsible for the temperature variations in the CF. For example, the models in [15,24] use a simple Joule heating temperature model controlled by electronic current in the filament region given in as [22]: However, many Schottky-based models approximate room temperature for current conduction (300 K) [5,25,26]. The model outputs for [25,26] are shown in Figs. 7 and Fig. 8, respectively. Models based on generalised hopping do not consider temperature variations for current conduction. Ii can thus be concluded that CF temperature variation does not play a significant role in RRAM current conduction.

Physical parameters of CF and structure of device:
Physical parameters like channel length, filament thickness, location of Schottky junction and oxygen vacancy doping concentration have significant impact on the variation in Schottky barrier height whereby influence electronic current. The length of the doped or un-doped section of the filament is spotted as a variable to control the Schottky barrier height in Schottky conduction models. The application of input voltage enables the oxygen vacancies to move towards the TE to form a CF and this channel doping effect ϕ Bn0 is the major factor for reducing the barrier height. To incorporate ϕ Bn0 , several approaches are adopted in several models such as simply a constant value has been assigned on the basis of material property in [15] whereas ϕ Bn0 has been proportionally related to insulating volume, l t in [26], to doped region length in [25] and to undoped tunnelling length (g) in [13] which is given in the table, where, L 0 is the initial filament length and ϕ Bi is the initial barrier height, w is the length of doped region, D is the switching layer thickness. Ideality factor η is used in RRAM current conduction models to assimilate the deviation of junction characteristics from its ideal behaviour. Different approaches are used to assign the ideality factor in models. A constant value is assigned in [26] while two distinct values are assigned to this parameter to show the trap effect in ON state and OFF state of device operation in [25] (Table 2), where D its and D itm are interfacial trap densities in the semiconductor and insulator [25]. Value '1' for η indicates an ideal Schottky junction and '2' represents non-ideal behaviour due to interface traps. η is used in [5] to represent the variation in trap charging and discharging throughout the cycles. In the same model, the tunnelling probability factor (TPF) is introduced as a critical parameter which incorporates the effect of tunnelling in the model [5]. TPF depends on the thickness of the oxide layer while, other models with Schottky approach ignored this parameter.
Alternatively, current through the device is implemented with hopping conduction in [24,29] models, with gap length g and applied input voltage V as the control variables. However, in filamentary models, the total current in the device is associated with the shape of CF. In [17], a conical shape is approximated to the CF where the current is not only contributed by the top face of CF, but the lateral sides as well.

Table 2 Parameters used in Schottky modelling Parameter
Model from [15] Model from [26] Model from [14] Model from [13] image force lowering effect ϕ del q 3 Nψ 8π 2 ε Sϕb The position of Schottky contact in the device structure also has contribution in modelling approaches. Schottky junction concept can either be applied to the interface between the top electrode and insulator layer [15,25,26,28,35]. However, Schottky barrier junction is assigned at the interface between the switching and bulk layers in Pt/Ta 2 O 5−x /TaO 2−x /Ta device at HRS [23].
The position of Schottky junction in the device structure also has a contribution in modelling approaches. Schottky junction concept can either be applied to the interface between the top electrode and insulator layer [15,25,26].

Comparison of models with experimental results
This section compares the models discussed in the previous sections with experimental data to understand the efficiency of the previously discussed physical mechanisms to obtain the model.
The experimental results reported for Pt/Ta 2 O 5 /TaO x /Pt RRAM with a 4 nm switching layer is given in Fig. 9a and the corresponding model curves are shown in Figs. 9b-d. In Figs. 9bc, the current conduction mechanisms adopted for the models are the Schottky current conduction at HRS and ohmic conduction at LRS. The model in Fig. 9d further considers a state-independent trap-assisted tunnelling through the region around the CF along with Schottky at HRS. The Schottky barrier at the Pt and Ta 2 O 5 interface produces the asymmetric current during forward and reverse biases at HRS. The HRS current is a result of both the doping effect of V τ s and tunnelling in [13], while it is primarily dependent on the Schottky conduction in [25].
The multilevel capability of the Ta 2 O 5 /TaO x has been proven experimentally as in Fig. 10a and theoretically modelled shown in Fig. 10b as given in [26]. The multiple resistance states are the reflections of variations in the insulating volume of conduction filament with Schottky barrier height. Hence, Schottky conduction is derived as the major current mechanism in the model, while the resistive switching is resulting from the formation and rupture of CF originated from the physical mechanism modelled with a differential equation derived from the equation of motion of ions.
In [15], HRS and LRS are modelled by a Schottky contact with tunnelling current and Poole Frenkel conduction areal current and the generated I-V are symmetrical as a result (Fig. 11). The simulation result is shown to match well with the experimental measurement as shown in Fig. 11. Meanwhile, the model presented in [15] defined the total current as the sum of state dependent Schottky current contributed by ionic and electronic components in addition to the state-independent areal current that activated only at HRS. The state variable is considered as a function of average oxygen vacancy concentration which determines the resistive switching process in the model.
Alternatively, the experimental results and models based on HfO 2 [29] are plotted in Fig. 12, while instead of adopting Schottky mechanism, different physical approaches are implemented such as hopping current conduction and redox-based bipolar resistive switching, where RS is demonstrated with the gap dynamics due to the generation recombination and migration-based redox modelling approach. However, the gap evolution exponentially relates to hopping conduction along with the hyperbolic sinusoidal function of applied voltage. Similarly, the characteristics of AlO x /HfO 2 are also modelled with the same physical concepts of generalised hopping conduction [29].
To validate the model presented in [18], the plot of measured and calculated data given in Fig. 13 is analysed and found that the model is fitting well with experimental results where the current profile is derived by the temperature and material-dependent current conduction model where the variable thermal conductivity of the gap region is achieved with exponential law and variable resistivity is derived from the Poole conduction law, while constant values are assigned for resistivity and thermal conductivity of the metal. Further, the RS is modelled with thermally activated ion kinetics falling under the redox reaction mechanisms.
Similarly, the model shown in Fig. 14 is based on the driftdiffusion RS mechanism and electronic current is represented with carrier continuity expression and Fourier equation is used for Joule Fig. 9 Comparison between experimental and simulated I-V curves for Ta 2 τ 5 /Taτ x RRAM with Schottky conduction (a) Experimental result from [14], (b) Simulation result from [14], (c) Simulation result from [13], (d) Simulation result from [25]   Experiment and simulated I-V curves for TiN/HfO x /Hf:AIO x /TiN RRAM from [29] heating. However, in contrast to above-mentioned model, instead of resistivity, variable thermal conductivity as well as thermally activated electrical conductance has been utilised to model the electronic current where the magnitude of electronic current changes with transition of oxide material from insulator to metallic phase [30].
The simulation results of model based on Pd/Ta 2 O 5 /TaO x /Pd structure are given in Fig. 15, where the current conduction is modelled with the basic current continuity approximation where the variation of current is linear with the applied voltage and the RS is obtained by the drift diffusion-based approach. Even though the model describes the fine tuning of the hopping distance parameter (a) to select analogue or digital switching behaviour, the measured data shown in Fig. 15 depict analogue switching characteristics, while the model displays a digital switching behaviour [20].
To validate the physics-based compact model of TiN/HfO2/TiN stack with real device characteristics, an analysis has been done and found that the model is in well agreement with the experimental data showing the suitability of applied physical mechanisms such as the redox-based resistive switching and the modelling of co-existing current mechanisms such as tunnelling through the pristine oxide structure and the ohmic conduction through the conduction filament as well as sub-oxide phase. Additionally, the model reports the importance of temperature on RS process and modelled with heat equation. In addition to the modelling of SET/RESET, the model given in [28] included the characterisation of electroforming stage in which highly resistive pristine oxide is converted to conductive channel on the application of high bias and the concept has been modelled with the rate of growth of sub-oxide defined with the redox mechanism (Fig. 16).
The models presented in this review have been compared with their respective experimental results. Due to the different materials and physical mechanisms adopted, a direct comparison between models is unfeasible. An error analysis is however carried out to determine each model's accuracy compared to the experiments. The results are listed in Table 3.

Conclusion
A comprehensive analysis of physics-based models has been carried out by identifying various approaches and parameters adopted in resistive switching and current conduction mechanisms. The investigated approaches to modelling RS are drift-diffusion and redox-based mechanisms. As for current conduction, Schottky conduction and a general hopping mechanism are analysed along with their required parameters. This review has found that the RS process is dependent on the electric field and temperature regardless of the mechanism. Meanwhile for current conduction mechanisms, the variation in temperature does not contribute significantly in the models. The CF geometry is found to have an important role on RS behaviour as well as current conduction, with 2D and 3D models which calculate both radius and gap length variations generally exhibiting more accurate results.