Physics inspired compact modelling of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {BiFeO}_3$$\end{document}BiFeO3 based memristors

With the advent of the Internet of Things, nanoelectronic devices or memristors have been the subject of significant interest for use as new hardware security primitives. Among the several available memristors, BiFe\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{O}_{3}$$\end{document}O3 (BFO)-based electroforming-free memristors have attracted considerable attention due to their excellent properties, such as long retention time, self-rectification, intrinsic stochasticity, and fast switching. They have been actively investigated for use in physical unclonable function (PUF) key storage modules, artificial synapses in neural networks, nonvolatile resistive switches, and reconfigurable logic applications. In this work, we present a physics-inspired 1D compact model of a BFO memristor to understand its implementation for such applications (mainly PUFs) and perform circuit simulations. The resistive switching based on electric field-driven vacancy migration and intrinsic stochastic behaviour of the BFO memristor are modelled using the cloud-in-a-cell scheme. The experimental current–voltage characteristics of the BFO memristor are successfully reproduced. The response of the BFO memristor to changes in electrical properties, environmental properties (such as temperature) and stress are analyzed and consistant with experimental results.


Introduction
It is highly appreciated how the Internet of Things (IoT) has inevitably integrated into our lives, making it more convenient and efficient.However, with the expansion and vast diffusion of connected devices in the IoT, cybersecurity concerns have also increased.The privacy of individuals, companies and institutions have been highly compromised 1,2 .Unfortunately, the classical security solutions (software-level mathematical or algorithmic solutions) are no longer sufficient to secure modern-day applications.The increasing physical and side-channel attacks necessitate the need for alternative solutions 3 .
Researchers and engineers have shifted their focus towards finding hardware-level solutions to address security-related challenges in recent times.The hardware-level solutions include the new nano-electronic devices, such as memristive devices or memristors, spintronics, or carbon nanotubes 4 .Explicitly, memristive devices are foreseen as promising candidates for future hardware security applications mainly because of their special properties, such as low power consumption, scalability to the nano grade, fast switching, large off/on ratio, good endurance and reliability [5][6][7] .Also known as the resistive switching random access memory (ReRAMs), these memristive devices are two-terminal devices whose resistance can be changed by applying a suitable electrical input.Apart from the features mentioned above, the switching mechanisms in these devices are intrinsically stochastic, which make these devices highly suitable for hardware security applications like physical unclonable functions (PUFs) 8,9 , true random number generators (TRNGs) 10 , and hash functions 11 .
So far, many devices have been reported that exhibit resistive switching behaviour; however, in the present work, we focus on the devices where the resistive switching is triggered by ionic motion driven by an electric field 12 .These devices can be either filamentary-type devices involving filaments' formation or interface-type (also called non-filamentary) devices involving the movement of charged defects.As mentioned by Du et al. 6,13 , the high currents induced in filamentary devices during the electroforming process can damage or destroy the device via thermodynamic dielectric breakdown, reducing the reliability of the device.In order to avoid the electroforming process, interface-type devices such as BiFeO 3 (BFO)-based memristors [14][15][16] , double-barrier memristive devices (DBMD) 17 are preferred.The BFO memristive devices have been intensively studied in the memristive community because they exhibit excellent characteristics such as electroforming free switching, long retention time, good endurance, and also offer multistage switching.These features make the BFO device highly recommended for implementing future hardware security applications such as PUFs and TRNGs.
Furthermore, the development of existing and new memristive devices for hardware security applications requires a precise understanding of their physical behaviour.It is often challenging to determine the exact switching mechanism using experimental or diagnostic methods.Therefore, simulation models are developed that can contribute significantly to understanding the behaviour of such devices.On the one hand, multi-dimensional computational models (such as 3D kinetic Monte-Carlo) are exploited for an in-depth understanding of resistive switching and moderately include real-world devices' physical and chemical processes and stochastic behaviour [18][19][20] .However, they are computationally very expensive and, therefore, cannot be used for performing circuit simulations.On the other hand, there are compact or concentrated models based on mathematical formulae.They are fast but do not include the physical and chemical processes in the device, and often, they do not include the intrinsic stochasticity found in these devices [21][22][23] .
Unlike the state-of-the-art models, we propose a circuit simulator-compatible distributed model for BFO memristor that considers the advantages of both the models mentioned above.It is a one-dimensional (1D) compact model including more or less realistic physics and the experimentally observed stochastic behaviour i.e., the cycle-to-cycle (C2C) variability and device-to-device (D2D) variability observed in BFO memristors.A kinetic Cloud-in-a-cell (CIC) 24,25 scheme is used to simulate the resistive switching mechanism based on the ion/vacancy transport.Although it is a distributed model, because we resolve it in a 1D space, it is computationally less demanding and fast.The model is primarily considered to provide an interface between circuit designers and device developers, as shown in Fig. 1 26 .It is used to explore the electrical properties of BFO as entropy sources and the effects of physical variables such as temperature and voltage on the entropy sources.The model can be extended further to investigate the performance of BFO memristive devices for hardware security applications by performing circuit simulations of memristive-PUFs or memristive-TRNGs with a SPICE-like circuit simulator.
The manuscript is divided into three sections.First, the simulation approach is discussed in detail, explaining the BFO memristive device and its current mechanisms.Then, the simulation results based on the experimentally determined electrical parameters of BFO are discussed and compared with the experimental findings.Finally, an overall summary and significant findings from the current work are provided in the conclusion section.

Simulation approach
A polycrystalline perovskite structured Au/BiFeO 3 (BFO)/Pt/Ti memristive device and its equivalent circuit are shown in Fig. 2a 15 .A BFO memristive device is regarded as an interface-type memristive device with BiFeO 3 as the primary layer, a rectifying Au/BFO top Schottky contact with nearly unflexible barrier height (D t ), and a rectifying and/or Ohmic BFO/Pt/Ti bottom Schottky contact with flexible barrier height (D b ).Previous studies have shown that the resistive switching behaviour in BFO can be described using the ferroelectric polarisation mechanism [27][28][29] or the oxygen vacancy migration mechanism 15,30 .However, as reported by You et al. 30 and Du et al. 15 , the resistive switching behaviour in BFO is primarily due to the migration of positively charged defects, and the ferroelectric switching can be excluded based on polarisation-electric-field measurements.
The positively charged defects are identified as the oxygen vacancies (V + O ) present at the interstitial sites.When a positive write biased is applied to the device, the V + O vacancies migrate to the bottom BFO/Pt interface under the influence of the electric field, thereby lowering the barrier height of D b and setting the device to low resistance state (LRS).These V + O vacancies migrate back to their equilibrium positions for a negative write bias, re-setting the device to an high resistance state (HRS).The barrier height of D t is almost unaffected by the V + O vacancies migration.The other fixed ions in the BFO are the Fe 3+ ions and Ti 4+ ions.These two fixed ions do not participate directly in the switching process but mainly determine the activation energy (U A ) in the BFO.The Ti 4+ ions replace the Fe 3+ ions close to the BFO/Pt interface.Due to this, the activation energy is not uniform but increases gradually from 0.55 eV at Au/ BFO interface to 0.75 eV at BFO/PT interface, depending on the occupation of Ti 4+ ions.It is important to note here that with the diffusion of Ti 4+ ions into the BiFeO 3 layer, the activation energy near the bottom electrode is increased to the point where vacancies can be easily trapped, thus improving the retention and endurance of the BFO memristor 31 .
According to You et al. 30 and Du et al. 15 , the V + O migration can be explained as: (a) the back and forth movement of V + O vacancies in the presence of an electric field with a certain drift velocity and (b) the trapping or de-trapping of V + O vacancies in traps or potential wells formed by the Ti 4+ fixed ions.A comparable V + O migration is implemented in this paper using the CIC scheme-based simulation model, which is explained in detail below.The parameters and their values used in the simulation of the BFO are given in Table .1.
The simulation approach is adapted from Yarragolla et al. 25 .The approach combines Newton's laws of motion with the kinetic cloud-in-a-cell (CIC) scheme to couple the vacancy transport with the electric field in the BiFeO 3 layer.The simulation scenario is depicted in Fig. 2b.The figure shows a BFO modelled on a 1D computational grid with equally and randomly arranged mobile positively charged and fixed negatively charged defects.The negatively charged defects are the fixed oxygen ions (V − O ), assumed only for having a neutral charge in BFO.Initially, all parameters listed in Table .1 are mostly constants and are defined manually, except for parameters such as defect density and length of the device, which vary slightly depending on the device fabrication process.These variables are selected from a collection of randomly generated values with a truncated Gaussian distribution.Then an equal number of positively charged and negatively charged defects are randomly distributed across the computational domain.The initial activation energy (U A ) throughout the computational domain is also set manually so that it gradually changes from 0.55 eV at the Au interface to 0.75 eV at the Pt/Ti interface.After the initialisation process is complete, the input voltage (V Device ) is specified, and the effective activation energy, which changes as a function of V Device , is calculated using the following line equation: where x is a position in the BFO layer, l BFO is the length of BFO, and λ U is the fitting parameter that determines the rate of change of U A,eff .λ U can be any number between 0 and 1.In the next step, the electric potential (φ ) and electric fields (E) within BFO are then calculated using the following equations, Here, ρ is the charge density, and ε is the permittivity of BFO.For this, Dirichlet boundary conditions are used, which are calculated using the voltage drops at the top and bottom interfaces by applying Kirchhoff's voltage law and Kirchhoff's current law to the equivalent circuit shown in Fig. 2a 25 .
where V SC,t/b ,V BFO are the voltage drops across the top and bottom Schottky barrier and BFO, respectively.Similarly, I SC,t/b , I BFO are the currents across the top and bottom Schottky barrier and BFO, respectively.An iterative scheme mentioned by Yarragolla et al. is used to calculate the voltage drops and currents.Followed by this, the electric potentials at Au/BFO and BFO/Pt interfaces are given by φ Au/BFO = V Device −V SC,t and φ BFO/Pt/Ti = −V SC,b , respectively.The currents mentioned in Eq.( 5) can be calculated as follows.The resistive switching in BFO is mainly attributed to the change in Schottky barrier properties at the metal/oxide interfaces.As seen in Fig. 2a, both diodes are initially in rectifying mode with forward-biased D t and reverse-biased D b .When a positive biased is applied, the V + O vacancies drift in the presence of an electric field, changing the V + O vacancies concentration at the BFO/Pt interface.This change in V + O vacancy concentration significantly decreases the barrier height of D b , thus making it conducting.D b is forward biased during a RESET process and changes its mode from conducting to rectifying, while D t is reversed biased and still in the rectifying mode.
In general, the current through the Schottky barrier is determined by thermionic emission (TE), thermionic field emission (TFE), direct tunnelling or a combination of these mechanisms.However, for a BFO  properties due to the two metal-semiconductor contacts at the top and bottom, it is sufficient to consider the mechanisms of thermionic emission and thermionic field emission.Since direct tunnelling does not lead to rectifying properties, it can be excluded as one of the dominant mechanisms.Therefore, the current across D t and D b is calculated using the Richardson-Schottky equation 32 , which is based on the TE mechanism and further modified by considering an effective ideality factor (n eff ) greater than one to account for the additional current tunneling through the barrier (i.e., described by TFE) 32,33 .The final current equation for the Schottky barriers is given by The reverse current, I R for forward and reverse bias conditions are respectively given by Here k B is the Boltzmann constant, T is the temperature, and A * = 1.20173 × 10 6 A/(m 2 K 2 ) is the effective Richardson constant.The effective ideality factor, n eff = n 0 (1 + λ n q(t)) and the effective Schottky barrier height, Φ eff = Φ 0 (1 + λ b q(t)), depend on the internal state of the device, q(t).Φ 0 is the initial Schottky barrier height, and n 0 is the initial ideality factor.λ b and λ n are the fitting parameters that define the rate of change of effective barrier height and effective ideality factor, respectively.The current flowing through the central BFO region can be calculated by applying Ohm's law.
where σ is the conductivity of BFO and l BFO is the length of BFO.Using Eq.( 5), the above Eq.( 9) can be modified as Using the electric field from Eq.( 3) and the effective activation energy, the drift velocity of V + O is given as follows 34,35 : where ν 0 is the phonon frequency, U A,eff is the effective activation energy, d is the jump distance (lattice constant), z is the charge number, k B is the Boltzmann constant, T is the temperature, and e is the elementary charge.Solving this velocity equation, we obtain a deterministic position of the V + O vacancies.But, in order to account for the intrinsic stochastic behaviour of the BFO device, the position of every i th V + O vacancy is artificially perturbed and given by Here r is the random number between 0 and 1, and δ is the maximum random displacement.After every iteration of the ion movement, the average relative distance between the current position of the vacancies, d(t) and their initial position, dr , determines the internal state of the device,

Results and discussion
First, the 1D compact model of BFO device described above was validated by comparing the calculated current-voltage characteristics (I-V curve) with the experimentally obtained I-V curve.For this purpose, the simulation model was initialized with the parameters listed in Table 1.In this state, the device is in its HRS.A voltage sweep with a ramp from 0 V through 8.5 V to 0 V was applied to set the device to an LRS, and then from 0 V through -8.5 V to 0 V to reset it back to HRS.For both set and reset process, a sweeping velocity of 0.36 V per 100 ms was applied.The change in applied voltage (V device ) with time and the resultant I-V curves are shown in Fig. 3(a).From the I-V curves, it can be seen that the model reproduces the analogue behaviour of BFO very well.The calculated I-V curve for the SET process (0 V to 8.5 V to 0 V) agrees with the experimental results both qualitatively and quantitatively.However, although the I-V curve for RESET agrees quite well quantitatively with the experimental I-V curve, the model does not entirely capture the effect of non-zero crossing (at -3 V) observed in the experimental I-V curve.Sun et al. 36 attributed this type of non-zero crossing behaviour of memristors to three mechanisms.They are the capacitive effect, the ferroelectric polarisation effect, and the presence of an internal electromotive force.Since the ferroelectric polarisation effects are already ruled out for a BFO, the possible reason for non-zero crossing in a BFO can most likely be the capacitive effects.The capacitance-voltage measurements of Shuai et al. 37,38 showed the presence of such capacitive effects in BFO.They indicated the presence of simultaneous resistive and capacitive switching in BFO, with HRS corresponding to the low capacitance state and vice versa.Fig. 3(b) shows the change in activation energy (U A ) across the BFO memristor as a function of time and V device shown in Fig. 3(a).As mentioned earlier, the activation energy across the BFO is not homogeneous but gradually increases from 0.55 eV to 0.75 eV between 550 nm and 600 nm.This inhomogeneous U A is due to the presence of the Ti region between 550 nm and 600 nm and plays an essential role in the vacancy transport responsible for the resistive switching behaviour of BFO.Vacancy transport is illustrated in Fig. 3(c .2. (d) The change in internal state of the device (q(t)) for D2D.The maximum applied voltage in both plots is ±8.5 V vacancies between 2 s and 3.2 s sets the device to an LRS with a nearly constant average relative distance ( d(t) ≈ 562 nm).Moreover, the change in the position of the V + O vacancies is so small that they can be assumed to be almost in a trapped state.When a negative write bias is applied after 4 s during the reset process, the activation energy at the BFO/Pt interface drops to about 0.6 eV (as shown in Fig. 3(b)), which is sufficient to return the V + O vacancies to their equilibrium position and reset the device to HRS.At the end of the reset process (at 8 s), the average relative distance also drops to an initial value of about 268 nm.
The intrinsic stochastic behaviour of BFO is measured in terms of spatial (device to device) and temporal (cycle to cycle) variability.First, the four I -V curves in Fig. 4(a) show the temporal variability of the BFO.The curves were obtained for four consecutive voltage sweeps shown in the inset of Fig. 4(a), and the initial conditions used to simulate these curves are given in the first row of Table .2. Although, the I -V curves calculated for each applied voltage cycle follow a similar trend, they slightly vary from each other.As observed from Fig. 4(b), the slight variation in the I -V curves is likely due to the change in internal state of the device (q(t)) from cycle to cycle.This change in q(t) is actually due to the random movement of the vacancies as calculated using Eq.( 12).Second, the spatial variability in BFO is illustrated using the four I -V curves in Fig. 4(c) for a single voltage sweep shown in the inset.The initial conditions used to simulate the four curves are given in Table .2.As can be observed, the I -V curves showing spatial variability are more clearly separated from each other than the I -V curves showing temporal variability.In general, based on experimental observations, the temporal variability, is much lower for interface-type memristive devices such as BFO and DBMD 25 .Moreover, the difference in the I -V curves showing spatial variability could be mainly due to the different initial conditions used to simulate each curve.The initial conditions more or less directly effect q(t), so for this reason we observe a change in q(t) as seen in Fig. 4(d).However, apart from the above-mentioned reasons, it is predicted that, several other unknown physical or chemical phenomena may also contribute to this spatial and temporal variability in a true BFO memristor.
The voltage drops across different BFO memristor regions are plotted in Fig. 5(a eff and n t eff increase with a positive bias and decrease with a negative bias.However, their change is almost negligible; therefore, D t can be considered non-flexible and constantly rectifying. Moreover, the initial Shottky barrier heights are considered as one of the entropy sources in hardware security applications.So, it is also important to check the behaviour of a BFO memristor to change in initial Schottky barrier heights.The graph in Fig. 5(c) shows the I -V curves with different combinations of Φ t 0 and Φ b 0 .Ideally, the barrier height can be increased or decreased by reducing or increasing the doping concentration, respectively 39 .Increasing Φ t/b 0 from 0.75 eV to 0.9 eV increases the energy required for the electrons to cross the barrier, which reduces the current through the BFO.In contrast, if Φ t/b 0 is decreased, the electrons can move more easily, which increases the current through the BFO.Also, a change in Φ b 0 mainly affects the right loop of the I -V curves, while a similar change in Φ t 0 , affects the left loop of the I -V curve.Therefore, as mentioned by Du et al. 15 and observed in Fig. 5, we can conclude that the set current is mainly determined by the barrier height and the voltage drop across D b , and the reset current by the barrier height and the voltage drop across D t .
The response of BFO memristor to changing environmental conditions (e.g., temperature) and stress (e.g., excessive voltage) is illustrated in Figs. 6 and 7.These factors are important when considering BFO for neuromorphic circuits, hardware security and non-volatile memory applications.Fig. 6 shows the simulated and experimental temperature-dependent I-V curves.The results are obtained for a writing bias of ±11 V and different temperatures increasing from 298 K to 348 K.According to Eq.( 11), the velocity of the vacancies (ν D ) increases with increasing temperature.As ν D increases, the maximum voltage required to switch the device to LRS reduces as shown in Fig. 6(b).However, as observed experimentally in Fig. 6(a), the maximum voltage required for switching is the same for all temperatures.No change in the switching voltage suggests that there might be some frictional force acting on the vacancies that decrease their velocity (ν D ) with increasing temperature.This Device dr (nm) l BFO (nm) ρ (cm     The measured and simulated I -V curves at room temperature for different maximum applied voltages are shown in Fig. 7.The voltage profiles used to obtain the plots are shown in Fig. 7(a).The observations from Fig. 7(b) indicate that the shape of the I -V curve in the set and reset direction strongly depend on the maximum applied voltage.At higher applied maximum voltages, the shape of the hysteresis is relatively wider compared to the hysteresis obtained at lower maximum V device .The simulated I -V curves in Fig. 7(c) also show a broadening of I -V curves with increasing maximum V device and quantitatively match quite well with the experimental results.However, there is a discrepancy in the simulated I -V curves between voltages -1.5 V to 1.5 V.As mentioned earlier, this discrepancy could be due to capacitive effects or the presence of an internal electromotive force.Furthermore, when the applied voltage is very low, the vacancies do not receive enough energy to drift to the BFO/Pt interface, resulting in switching failure, i.e. the device does not switch to LRS and remains in a HRS.This can be seen in Fig. 7(d), which shows the change in average distance for different maximum applied voltages.The average distance is only tracked for the positive V −device cycle since we are interested in observing the vacancies drift towards the BFO/Pt interface, which mainly contributes toward switching the device to LRS.For low maximum V device from 3 V to 5 V, d(t) is less than 520 nm, which means that certain vacancies do not reach the BFO/Pt interface and are not trapped.This can be better understood from Fig. 7(e) that shows the vacancy transport in terms of charge density for a positive applied voltage cycle.The retention time of the device is also severely affected due to the untrapped vacancies in the 3 V and 5 V plots.
The retention characteristics of BFO memristor are investigated using the proposed stochastic model.The retention tests were performed by switching the device to LRS or HRS, which are the initial states for this particular study.Then the externally applied voltage was switched off, and the diffusion of the vacancies was recorded.Since the experimental results are obtained for a real BFO device (i.e., 3D), we had to interpolate the 1D model to 3D.In a 3D model, the vacancies can generally move in six directions, while in a 1D space, the vacancies can move either back or forth.So, to fit the simulated results with the experimental results, we use a fitting parameter that restricts the movement of vacancies in BFO by reducing their jump attempts.For this, we randomly pick a number, β between 0 and 1, and use the following relation for moving the vacancies: where ν D is the drift velocity of the vacancies given by Eq. (11).The development of the device current was recorded every 10 s with a read voltage of 2 V at room temperature.The simulated and experimental results for a total of 3000 cycles/pulses are shown in Fig. 8.As observed, BFO shows good retention characteristics.The HRS is stable, and no significant change was observed during the 3000 cycles.For the LRS, the good retention is primarily due to the diffusion of Ti 4+ ions that increases the activation energy at the BFO/Pt interface.This high activation energy limits the movement of vacancies .ie., the vacancies get trapped.However, BFO shows a degradation in the LRS current until approximately the 700 th cycle, i.e., two hours before stabilizing.Through simulations, this degradation was found to be due to the diffusion of some vacancies away from the BFO/Pt interface that were not trapped.The relative average distance, d(t), decreased from 564 nm to 542 nm, which increased Φ b eff from 0.62 eV to 0.68 eV, thereby decreasing the current.One possible way to improve retention could be to ensure that all the vacancies are properly trapped by increasing the Ti fluence 31 .The second possibility would be to improve the BFO surface using low-energy Ar + ion irradiation, as suggested by Shuai et al. 37 .

Conclusion
With billions of electronic devices connected to the Internet worldwide, low-power, nanoscale memristive devices are considered favourable devices for a more secure Internet of Things.Due to their stochastic behaviour, these memristive devices are ideally suited for hardware security applications such as PUFs, TRNGs, and cryptographic algorithms.The BiFeO 3 -based memristive devices are deemed to be suitable for such applications.In the proposed work, a physics-inspired compact 1D model of an Au/BiFeO 3 (BFO)/Pt/Ti memristor is developed for circuit-level simulations in the field of hardware security applications and neuromorphic circuits.The model successfully simulates resistive switching based on electric field-driven migration of oxygen vacancies and accounts for the intrinsic stochastic nature of the BFO memristor.A cloud-in-a-cell scheme is used in which Newton's laws are consistently coupled with the Poisson solver.The simulated current-voltage characteristics of the BFO memristor obtained with this scheme agree well with the experimental results.It was found that the set current is mainly determined by the Schottky barrier height and the voltage drop across the BFO/Pt interface, while the reset current is determined by the Schottky barrier height and the voltage drop across the Au/BFO interface.In addition, based on the observations of the simulated and experimental temperature-dependent current-voltage characteristics, we anticipate the presence of a frictional force acting on the oxygen vacancies that increases with temperature.The simulated and experimental results illustrating the effects of temperature, stress, and the retention characteristics of BFO show reasonable agreement.The proposed model is highly efficient and reliable as it consists of various parameters that can be easily tuned to match the experimental results, and the degree of stochasticity can also be adjusted.To further comprehend the switching process in the BFO memristor, the 1D model could be extended to a 2D or 3D model that better represents the real-world BFO memristor.

Methods
Experiments.Polycrystalline, 600 nm thick BFO functional thin film have been grown by pulsed laser deposition on Pt/Ti/SiO 2 /Si substrates.Circular Au top contacts with thickness of 180 nm have been magnetron sputtered on the BFO thin films using a shadow mask 40,41 .All the electrical measurements, illustrated in this work, were recorded using a Keithley source meter 2400, which were connected to the PC via GPIB cables and can be controlled through LabVIEW program.
Simulations.An in-house model developed using C programming language.

Figure 1 .
Figure 1.A flowchart illustrating the goal of the proposed work.

Figure 3 .
Figure 3.The resistive switching process in BFO memristor.(a) calculated and experimentally obtained I-V characteristics of the BFO memristor, (b)the change in activation energy with time across the BFO device during the set and reset process, and (c) the ion transport in BFO memristor shown as a change in charge density (ρ) over time.The black dashed line indicates the absolute average distance d(t)

O 6 / 13
) by considering how the charge density in BFO changes during the simulation time of 8 s.Initially, the vacancies are randomly placed across the BFO computational domain, and when a voltage is applied, the vacancies move towards the BFO/Pt interface.Due to the very high activation energy (about 0.75 eV) at the BFO/Pt interface, the change in position of the V + O vacancies near the Au/BFO interface is insignificant.This insignificant change in the position of the V +

Figure 4 .
Figure 4.The simulated IV curves showing (a) cycle-to-cycle (C2C) variability obtained for four consecutive voltage sweeps with initial conditions given in the first row of Table.2.(b) The change in internal state of the device (q(t)) for C2C.(c) device-to-device (D2D) variability obtained for a single voltage sweep but with different initial conditions given in Table.2.(d) The change in internal state of the device (q(t)) for D2D.The maximum applied voltage in both plots is ±8.5 V ) to investigate the rectification process and the physical mechanisms behind the strong voltage dependence.The different regions include the two Schottky contacts at the Au/BFO (D t ) and BFO/Pt (D b ) interfaces, and the BiFeO 3 layer.For a positive bias, a back-to-back rectification is observed with a forward-biased D t and reverse-biased D b .As can be seen in Fig.5(a), although there is some voltage drop across the D t and BiFeO 3 layer, most of the voltage is blocked by the reverse-biased D b .Furthermore, for a negative bias, almost all the applied voltage is blocked by the reverse-biased D t , with a small voltage drop across the D b and a negligible one across the BFO layer.The changes in other properties of the Schottky contacts, such as the barrier height (Φ) and the ideality factor (n) during the sweep time, are also shown in Fig.5(b).The change in Φ and n are strongly influenced by the vacancy transport in the BFO.As vacancies move toward the BFO/Pt interface during positive bias, D b becomes non-rectifying, with Φ b eff decreasing from 0.85 eV to 0.62 eV and n b eff from 4.5 to 3.3.This allows electrons to flow easily across the barrier, increasing the current through the BFO memristor.With a negative bias, Φ b eff and n b eff increase as the vacancies move away from the BFO/Pt interface and D b becomes rectifying.On the other hand, Φ t

Figure 5 .
Figure 5. (a) Voltage drop across different regions of the BFO memristor.(b) The change in top and bottom effective Schottky barrier heights (Φ eff ) and effective ideality (n eff ) factor as a function of time.(c) The simulated I-V curves showing the effect of top and bottom Schottky barrier height on the set and reset process.Curves with Φ 0,t =0.8 eV overlap in the negative bias region.

Figure 6 .
Figure 6.(a) The experimentally obtained temperature-dependent current-voltage characteristics (I-V-T curves).(b) Simulated I-V-T curves obtained using the fitting parameter (λ T ).The parameters used in the simulation are given in Table.3.(c) I-V-T curves obtained for positive applied voltage and without λ T .The dashed lines indicate the voltage required to switch the device from HRS to LRS.The legend is same as mentioned in Fig.6(b).(d) The calculated average drift velocity (ν D ) with and without λ T at different temperatures.
Figure 6.(a) The experimentally obtained temperature-dependent current-voltage characteristics (I-V-T curves).(b) Simulated I-V-T curves obtained using the fitting parameter (λ T ).The parameters used in the simulation are given in Table.3.(c) I-V-T curves obtained for positive applied voltage and without λ T .The dashed lines indicate the voltage required to switch the device from HRS to LRS.The legend is same as mentioned in Fig.6(b).(d) The calculated average drift velocity (ν D ) with and without λ T at different temperatures.

3 .
Therefore, by considering the fitting term, we are able to mimic the temperature-dependent resistive switching in BFO, as shown in Fig.6(c).

Figure 7 .
Figure 7. (a) Input voltage source (V Device ) with different amplitudes (b) Experimental I -V curves, and (c) simulated I -V curves for different maximum V Device .(d) The change in average distance d(t) with time for a positive applied voltage cycle with different amplitudes.(e)The change in charge density over time for a positive applied voltage cycle with maximum voltage of 3 V, 5 V and 7 V.The legend for all plots is shown on the top right corner.

Figure 8 .
Figure 8.The comparison between simulated and experimental retention characteristics of a BFO memristor.The current evolution is recorded every 10 s for 3000 cycles for both LRS and HRS.

Table 1 .
15mristive device with strong rectifying Details of the simulation parameters15

Table 2 .
The parameters of the four devices used to determine the I-V curves in Fig.4

Table 3 .
The parameters used to simulate the I-V curves in Fig.6