1 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports Physics inspired compact modelling of BiFeO 3 based memristors Sahitya Yarragolla 1*, Nan Du 2,3*, Torben Hemke 1, Xianyue Zhao 2,3, Ziang Chen 2,3, Ilia Polian 4 & Thomas Mussenbrock 1 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, BiFeO 3 (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. 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 compromised1,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 solutions3. Researchers and engineers have shifted their focus towards finding hardware-level solutions to address secu- rity-related challenges in recent times. The hardware-level solutions include the new nano-electronic devices, such as memristive devices or memristors, spintronics, or carbon nanotubes4. 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 reliability5–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 functions11. 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 field12. 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 BiFeO3 (BFO)-based memristors14–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. OPEN 1Chair of Applied Electrodynamics and Plasma Technology, Ruhr University, Bochum, Germany. 2Institute for Solid State Physics, Friedrich Schiller University, Jena (FSUJ), Jena, Germany. 3Department of Quantum Detection, Leibniz Institute of Photonic Technology (IPHT), Jena, Germany. 4Institute of Computer Science and Computer Engineering, University of Stuttgart, Stuttgart, Germany. *email: sahitya.yarragolla@rub.de; nan.du@ uni-jena.de https://orcid.org/0000-0002-2973-4943 https://orcid.org/0000-0002-7775-7795 https://orcid.org/0000-0003-2436-5840 https://orcid.org/0000-0002-6563-2725 https://orcid.org/0000-0001-6445-4990 http://crossmark.crossref.org/dialog/?doi=10.1038/s41598-022-24439-4&domain=pdf 2 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ 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-dimen- sional 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 behaviour18–20. However, they are computationally very expensive and, therefore, cannot be used for perform- ing 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 devices21–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. 126. 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 (mem-PUFs) or memristive-TRNGs (mem-TRNGs) with a SPICE-like circuit simulator. It is important to mention that this paper serves as a basis for conducting circuit simulations of mem-PUFs and mem-TRNGs and is mainly intended to demonstrate the capabilities of the proposed model. Therefore, the scope of the paper is initially limited to the compact modelling and parametric investigation steps shown in Fig. 1. These steps provide circuit designers and device engineers with an insight into device physics. A simulation campaign and experiments of mem-PUFs and mem-TRNGs are planned as the next step. 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. Figure 1. A flowchart illustrating the goal of the proposed work. The non-grey areas indicate the main steps addressed in the current manuscript. 3 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ Simulation approach A polycrystalline perovskite structured Au/BiFeO3 (BFO)/Pt/Ti memristive device and its equivalent circuit are shown in Fig. 2a15. A BFO memristive device is regarded as an interface-type memristive device with BiFeO3 as the primary layer, a rectifying Au/BFO top Schottky contact with nearly unflexible barrier height ( Dt ), and a rectifying and/or Ohmic BFO/Pt/Ti bottom Schottky contact with flexible barrier height ( Db ). Previous studies have shown that the resistive switching behaviour in BFO can be described using the ferroelectric polarisation mechanism27–29 or the oxygen vacancy migration mechanism15,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 bias 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 Db 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 Dt is almost unaffected by the V+ O vacancies migration. The other fixed ions in the BFO are the Fe3+ ions and Ti4+ ions. These two fixed ions do not participate directly in the switching process but mainly determine the activation energy ( UA ) in the BFO. The Ti4+ ions replace the Fe3+ 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 Ti4+ ions. It is important to note here that with the diffusion of Ti4+ ions into the BiFeO3 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 memristor31. According to You et al.30 and Du et al.15, the V+ O migration can be explained as: (a) the back and forth move- ment 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 Ti4+ 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 BiFeO3 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 distribu- tion. Then an equal number of positively charged and negatively charged defects are randomly distributed across the computational domain. The initial activation energy (UA) 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 ( VDevice ) is specified, and the effective activation energy, which changes as a function of VDevice , is calculated using the following line equation: Au Pt/Ti BFO Ti4+ Fe3+ V+ o V− o I I V Device VSC,t(I) VBFO(I, RBFO) Dt −V SC,b(I) Db (a) 0.55 eV U A,eff U A,eff SET RESET BFO (600 nm) Ti+4 region (50 nm) Au (180 nm) Pt/Ti (190 nm) ≈ 0.75 eV 0.6 eV 0.55 eV (b) Figure 2. Simulation scenario of BFO memristor. (a) A Au/BiFeO3 (BFO)/Pt/Ti memristive device and it’s equivalent circuit with back-back schottky diodes and a variable resistor. (b) Ion transport and the potential structure across a 1D BFO memristor for set and reset process. 4 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ where x is a position in the BFO layer, lBFO is the length of BFO, and �U is the fitting parameter that determines the rate of change of UA,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. 2a25. where VSC,t/b,VBFO are the voltage drops across the top and bottom Schottky barrier and BFO, respectively. Similarly, ISC,t/b, IBFO 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 = VDevice − VSC,t and φBFO/Pt/Ti = −VSC,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 Dt and reverse-biased Db . 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 Db , thus making it conducting. Db is forward biased during a RESET process and changes its mode from conducting to rectifying, while Dt 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 memristive device with strong rectifying 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 tun- nelling does not lead to rectifying properties, it can be excluded as one of the dominant mechanisms. Therefore, the current across Dt and Db is calculated using the Richardson–Schottky equation32, which is based on the TE mechanism and further modified by considering an effective ideality factor (neff ) 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, IR for forward and reverse bias conditions are respectively given by (1)UA,eff = UA + �UVDevice ( 1− x lBFO ) , (2) d dx ( ε dφ dx ) =− ρ, (3)E =− dφ dx . (4)VDevice =VSC,t + VBFO − VSC,b, (5)ISC,t = IBFO = −ISC,b = I , (6)ISC = IR ( exp { eVSC neff kBT } − 1 ) . Table 1. Details of the simulation parameters.15 Physical quantity Symbol Value Temperature (K) T 298 Phonon frequency (Hz) ν0 1× 10 12 Lattice constant (m) d 0.56× 10 −9 Device area ( mm 2) Ad 0.04 Relative permittivity of BiFeO3 εr 52.0 Conductivity of BiFeO3 ( �m) σ 7.0× 10 −4 Length of BFO (m) lBFO 600× 10 −9 Defect density ( cm−3) ρ 2× 10 16 Top Schottky barrier height (eV) �t 0 0.8 Top Schottky barrier ideality factor n t 0 4.2 Bottom Schottky barrier height (eV) �b 0 0.85 Bottom Schottky barrier ideality factor n b 0 4.5 5 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ Here kB is the Boltzmann constant, T is the temperature, and A∗= 1.20173× 106 A/(m2 K2) is the effective Richardson constant. The effective ideality factor, neff = n0(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 n0 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 lBFO is the length of BFO. Using Eq. (5), the above Eq.  (9) can be modi- fied as Using the electric field from Eq. (3) and the effective activation energy, the drift velocity of V+ O is given as follows34,35: (7)IR,VSC>0 =AdA ∗T2exp { −�eff kBT } , and (8)IR,VSC<0 =AdA ∗T2exp { −�eff kBT } exp { αr √ |VSC| kBT } . (9)IBFO = σAd VBFO lBFO , (10)VBFO = ISC,t lBFO σAd . (11)vD = ν0dexp ( − UA,eff kBT )( exp { |z|edE 2kBT } − exp { − |z|edE 2kBT }) , 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). 6 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ where ν0 is the phonon frequency, UA,eff is the effective activation energy, d is the jump distance (lattice constant), z is the charge number, kB 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 ith 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, d̄r , 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 cur- rent–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 through 8.5 V to 0 V was applied to set the device to an LRS, and then from 0 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 ( Vdevice ) with time and the resultant I–V curves are shown in Fig. 3a. 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 mem- ristors 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 pres- ence of simultaneous resistive and capacitive switching in BFO, with HRS corresponding to the low capacitance state and vice versa. Figure 3b shows the change in activation energy ( UA ) across the BFO memristor as a function of time and Vdevice shown in Fig. 3a. As mentioned earlier, the activation energy across the BFO is not homogeneous but gradually increases from 0.55 to 0.75 eV between 550 and 600 nm.This inhomogeneous UA is due to the presence of the Ti region between 550 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. 3c 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+ O vacancies between 2 and 3.2 s sets the device to an LRS with a nearly constant average relative distance ( ̄d(t) ≈ 562 nm). (12) x̄i ︸︷︷︸ stochastic = x̄i ︸︷︷︸ deterministic + (r − 0.5)δx̄i ︸ ︷︷ ︸ stochastic . (13)q(t) = d̄(t)− d̄r d̄r . 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. 7 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ 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 ( neff ) 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. Table 2. The parameters of the four devices used to determine the I–V curves in Fig. 4. Device d̄r (nm) lBFO (nm) ρ ( cm−3) 1 295.4 600 2× 10 16 2 306.89 588.2 2.6× 10 16 3 299.001 601.5 2.1× 10 16 4 302.23 586.9 1.5× 10 16 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. 3b), 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. 4a show the temporal variability of the BFO. The curves were obtained for four consecutive voltage sweeps shown in the inset of Fig. 4a, and the initial condi- tions 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. 4b, Table 3. The parameters used to simulate the I–V curves in Fig. 6. T (K) εr σ ( �m) �T 298 52 8× 10 −4 0.0 313 60 9× 10 −4 0.005 323 72 1× 10 −3 0.01 333 100 3.5× 10 −3 0.04 343 132 5× 10 −3 0.06 348 145 7× 10 −3 0.062 8 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ 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. 4c 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 compared to the spatial variability for interface-type memristive devices such as BFO and DBMD25. This difference in the temporal and spatial variability of the BFO is mainly essential for implementing mem- PUFs because multiple runs of the same PUF should give almost the same response, but it should be different for different copies of the same circuit. It means that (a) each memristive device in the PUF crossbar should be different from each other (i.e., more spatial variability), and (b) each memristive device should produce nearly the same output every time for the same supply voltage (i.e., comparatively less temporal variability)6. Moreover, the difference in the I–V curves showing spatial variability could be mainly due to the different initial condi- tions 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. 4d. 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. 5a to investigate the rectifica- tion process and the physical mechanisms behind the strong voltage dependence. The different regions include the two Schottky contacts at the Au/BFO ( Dt ) and BFO/Pt ( Db ) interfaces, and the BiFeO3 layer. For a positive bias, a back-to-back rectification is observed with a forward-biased Dt and reverse-biased Db . As can be seen in Fig. 5a, although there is some voltage drop across the Dt and BiFeO3 layer, most of the voltage is blocked by the reverse-biased Db . Furthermore, for a negative bias, almost all the applied voltage is blocked by the reverse-biased Dt , with a small voltage drop across the Db 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 fac- tor (n) during the sweep time, are also shown in Fig. 5b. 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, Db becomes non-rectifying, with �b eff decreasing from 0.85 to 0.62 eV and nbeff 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 nbeff increase as the vacancies move away from the BFO/Pt interface and Db becomes rectifying. On the other hand, �t eff and nteff increase with a positive bias and decrease with a negative bias. However, their change is almost negligible; therefore, Dt can be considered non-flexible and constantly rectifying. Moreover, the initial Schottky barrier heights are considered as one of the entropy sources in hardware secu- rity 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. 5c 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, respectively39. Increasing �t/b 0 from 0.75 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 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 (b). (d) The calculated average drift velocity ( νD ) with and without �T at different temperatures. 9 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ Figure 7. (a) Input voltage source ( VDevice ) with different amplitudes. (b) Experimental I–V curves, and (c) simulated I–V curves for different maximum VDevice . (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. 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 Db , and the reset current by the barrier height and the voltage drop across Dt. 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 neu- romorphic 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 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. 6b. However, as observed experimentally in Fig. 6a, 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 frictional force could be due to the following reasons: (a) With the increasing temperature, more Ti4+ ions may diffuse into the BiFeO3 layer due to thermal diffusion. The more Ti diffuses into the BFO layer, the higher the activation energy, resulting in a decrease of νD , (b) The collision rate between particles increases with increasing temperature, which affects the motion of particles in a device, and (c) the presence of temperature-dependent ferroelectric polarisation switching. To account for the frictional force in BFO, we modify νD by using a fitting parameter to match the experimental results. So, νD = νD (1 − �T ) where �T could be any number between 0 and 1. In this way, we can reduce νD of all vacancies for different temperatures. This can be seen in Fig. 6d, which shows the average νD with and without including the fitting parameter for different temperatures. The parameters used to simulate I–V curves in Fig. 6c are given in Table 3. Therefore, by considering the fitting term, we are able to mimic the temperature-dependent resistive switching in BFO, as shown in Fig. 6c. 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. 7a. The observations from Fig. 7b 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 Vdevice . The simulated I–V curves in Fig. 7c also show a broadening of I–V curves with increasing maximum Vdevice and quantitatively match quite well with the experimental results. However, there is a discrepancy in the simulated I–V curves between voltages − 1.5 and 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. 7d, which shows the change in average distance for different maximum applied volt- ages. 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 Vdevice from 3 to 5 V, d̄(t) is less than 520 nm, which means that certain vacancies do not reach 10 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ the BFO/Pt interface and are not trapped. This can be better understood from Fig. 7e 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 particu- lar 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 Ti4+ ions that increases the activation energy at the BFO/Pt interface. This high activation energy limits the movement of vacancies i.e., the vacancies get trapped. However, BFO shows a degradation in the LRS current until approximately the 700th cycle, i.e., 2 h before stabilizing. Through simula- tions, 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 to 542 nm, which increased �b eff from 0.62 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 fluence31. The second possibility would be to improve the BFO surface using low-energy Ar+ ion irradiation, as suggested by Shuai et al.37. Finally, a cycle number-dependent plasticity measurement, according to Du et al.40. is carried out to assess the model’s performance further. The nominal thickness of the BFO used here is 500 nm, with the top and νD = νD, updated for 0 > β < 0.33, νD = νD, present for 0.33 > β < 1. Figure 9. Cycle-number dependent plasticity in BFO. The calculated and experimental synaptic change of (a) potentiation and (b) depression dynamics. G1 is the conductance measured during the first spike. 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. 11 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ bottom contacts thickness of 180 nm. For the potentiation spike sequence, an initialization pulse of − 6 V is used for the RESET process. Similarly, for the depression spike sequence, a pulse of 6 V is used for the SET process. A normalized synaptic weight change is calculated in Fig. 9 for 250 similar spikes. The amplitude of the pulse used for potentiation dynamics is 5 V, and for depression dynamics, it is − 5 V. A pulse width of 100 ms and the interval between two pulses of 20 ms is applied. As observed in Fig. 9 the experimental and calculated normalized synaptic weight change follow a similar trend. As mentioned by Du et al.40, it was also observed via simulations that with an increase in the potentiation spikes, the oxygen vacancies tend to move towards the Pt/ Ti electrode and thus increasing the conductance. For the depression spikes, the change in conductance is not so noticeable, and the oxygen vacancies tend to move back to their inertial position, causing the simultaneous switching of the device into HRS. 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 BiFeO3-based memristive devices are deemed to be suitable for such applications. In the proposed work, a physics-inspired compact 1D model of an Au/BiFeO3 (BFO)/Pt/Ti memristor is developed for circuit- level simulations in the field of hardware security applications and neuromorphic circuits. The model success- fully 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–volt- age 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 dep- osition on Pt/Ti/SiO2/Si substrates. Circular Au top contacts with thickness of 180 nm have been magnetron sputtered on the BFO thin films using a shadow mask41,42. 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. The codes developed and implemented are currently research codes that are not yet ready to be used by non-experts. Data availibility The datasets used and/or analysed during the current study available from the corresponding author on reason- able request. Received: 6 October 2022; Accepted: 15 November 2022 References 1. Chen, B. & Willems, F. M. J. Secret key generation over biased physical unclonable functions with polar codes. IEEE Internet Things J. 6, 435–445. https:// doi. org/ 10. 1109/ JIOT. 2018. 28645 94 (2019). 2. Weber, R. H. Internet of Things Vol. 12 (Springer, 2010). 3. Rajendran, S. & Rehman, M. Security of Internet of Things Nodes: Challenges, Attacks, and Countermeasures (Chapman and Hall, 2021). 4. Rose, G. S. et al. Nanoelectronics and Hardware Security 105–123 (Springer, 2014). 5. Pang, Y., Gao, B., Lin, B., Qian, H. & Wu, H. Memristors for hardware security applications. Adv. Electron. Mater. 5, 1800872. https:// doi. org/ 10. 1002/ aelm. 20180 0872 (2019). 6. Du, N., Schmidt, H. & Polian, I. Low-power emerging memristive designs towards secure hardware systems for applications in internet of things. Nano energy materials and devices for miniaturized electronics and smart systems.. Nano Mater. Sci. 3, 186–204. https:// doi. org/ 10. 1016/j. nanoms. 2021. 01. 001 (2021). 7. Rajendran, G., Banerjee, W., Chattopadhyay, A. & Aly, M. M. S. Application of resistive random access memory in hardware security: A review. Adv. Electron. Mater. 7, 2100536. https:// doi. org/ 10. 1002/ aelm. 20210 0536 (2021). 8. Uddin, M., Shanta, A. S., Badruddoja Majumder, M., Hasan, M. S. & Rose, G. S. Memristor crossbar puf based lightweight hardware security for iot. In 2019 IEEE International Conference on Consumer Electronics (ICCE), 1–4. https:// doi. org/ 10. 1109/ ICCE. 2019. 86619 12 (2019). 9. Ibrahim, H. M., Abunahla, H., Mohammad, B. & AlKhzaimi, H. Memristor-based puf for lightweight cryptographic randomness. Sci. Rep. 12, 8633. https:// doi. org/ 10. 1038/ s41598- 022- 11240-6 (2022). 10. Jiang, H. E. A. A novel true random number generator based on a stochastic diffusive memristor. Nat. Commun. 8, 882. https:// doi. org/ 10. 1038/ s41467- 017- 00869-x (2017). https://doi.org/10.1109/JIOT.2018.2864594 https://doi.org/10.1002/aelm.201800872 https://doi.org/10.1016/j.nanoms.2021.01.001 https://doi.org/10.1002/aelm.202100536 https://doi.org/10.1109/ICCE.2019.8661912 https://doi.org/10.1109/ICCE.2019.8661912 https://doi.org/10.1038/s41598-022-11240-6 https://doi.org/10.1038/s41467-017-00869-x https://doi.org/10.1038/s41467-017-00869-x 12 Vol:.(1234567890) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ 11. Azriel, L. & Kvatinsky, S. Towards a memristive hardware secure hash function (memhash). In 2017 IEEE International Symposium on Hardware Oriented Security and Trust (HOST), 51–55. https:// doi. org/ 10. 1109/ HST. 2017. 79517 97 (2017). 12. Waser, R. & Aono, M. Nanoionics-based resistive switching memories. Nat. Mater. 6, 833–840. https:// doi. org/ 10. 1038/ nmat2 023 (2007). 13. Du, N. et al. Electroforming-free memristors for hardware security primitives. In 2019 IEEE 4th International Verification and Security Workshop (IVSW), 67–70. https:// doi. org/ 10. 1109/ IVSW. 2019. 88543 94 (2019). 14. Shuai, Y. et al. Key concepts behind forming-free resistive switching incorporated with rectifying transport properties. Sci. Rep. 3, 2208. https:// doi. org/ 10. 1038/ srep0 2208 (2013). 15. Du, N. et al. Field-driven hopping transport of oxygen vacancies in memristive oxide switches with interface-mediated resistive switching. Phys. Rev. Appl. 10, 054025. https:// doi. org/ 10. 1103/ PhysR evApp lied. 10. 054025 (2018). 16. You, T. et al. Exploiting memristive bifeo3 bilayer structures for compact sequential logics. Adv. Funct. Mater. 24, 3357–3365. https:// doi. org/ 10. 1002/ adfm. 20130 3365 (2014). 17. Hansen, M. et al. A double barrier memristive device. Sci. Rep. 5, 13753. https:// doi. org/ 10. 1038/ srep1 3753 (2015). 18. Dirkmann, S., Hansen, M., Ziegler, M., Kohlstedt, H. & Mussenbrock, T. The role of ion transport phenomena in memristive double barrier devices. Sci. Rep. 6, 35686. https:// doi. org/ 10. 1038/ srep3 5686 (2016). 19. Dirkmann, S., Kaiser, J., Wenger, C. & Mussenbrock, T. Filament growth and resistive switching in hafnium oxide memristive devices. ACS Appl. Mater. Interfaces 10, 14857–14868. https:// doi. org/ 10. 1021/ acsami. 7b198 36 (2018). 20. Abbaspour, E., Menzel, S. & Jungemann, C. Kmc simulation of the electroforming, set and reset processes in redox-based resistive switching devices. In 2016 International Conference on Simulation of Semiconductor Processes and Devices (SISPAD), 141–144. https:// doi. org/ 10. 1109/ SISPAD. 2016. 76051 67 (2016). 21. Solan, E. et al. An enhanced lumped element electrical model of a double barrier memristive device. J. Phys. D Appl. Phys. 50, 195102. https:// doi. org/ 10. 1088/ 1361- 6463/ aa69ae (2017). 22. Ambrogio, S., Balatti, S., Gilmer, D. C. & Ielmini, D. Analytical modeling of oxide-based bipolar resistive memories and comple- mentary resistive switches. IEEE Trans. Electron Devices 61, 2378–2386. https:// doi. org/ 10. 1109/ TED. 2014. 23255 31 (2014). 23. Bengel, C. et al. Variability-aware modeling of filamentary oxide-based bipolar resistive switching cells using spice level compact models. IEEE Trans. Circ. Syst. I Regul. Pap. 67, 4618–4630. https:// doi. org/ 10. 1109/ TCSI. 2020. 30185 02 (2020). 24. Laux, S. E. On particle-mesh coupling in Monte Carlo semiconductor device simulation. IEEE Trans. Comput. Aided Des. Integr. Circ. Syst. 15, 1266–1277. https:// doi. org/ 10. 1109/ 43. 541446 (1996). 25. Yarragolla, S. et al. Stochastic behavior of an interface-based memristive device. J. Appl. Phys. 131, 134304. https:// doi. org/ 10. 1063/5. 00840 85 (2022). 26. Saha, S. Compact Models for Integrated Circuit Design: Conventional Transistors and Beyond (Taylor and Francis Group, 2015). 27. Choi, T., Lee, S., Choi, Y. J., Kiryukhin, V. & Cheong, S.-W. Switchable ferroelectric diode and photovoltaic effect in BiFeO3. Science 324, 63–66. https:// doi. org/ 10. 1126/ scien ce. 11686 36 (2009). 28. Wang, C. et al. Switchable diode effect and ferroelectric resistive switching in epitaxial bifeo3 thin films. Appl. Phys. Lett. 98, 192901. https:// doi. org/ 10. 1063/1. 35898 14 (2011). 29. Lei, Y. et al. Ferroelectric and flexible barrier resistive switching of epitaxial bifeo3 films studied by temperature-dependent current and capacitance spectroscopy. J. Mater. Sci. Mater. Electron. 27, 7927–7932. https:// doi. org/ 10. 1007/ s10854- 016- 4784-y (2016). 30. You, T. et al. Bipolar electric-field enhanced trapping and detrapping of mobile donors in bifeo3 memristors. ACS Appl. Mater. Interfaces 6, 19758–19765. https:// doi. org/ 10. 1021/ am504 871g (2014). 31. You, T. et al. Engineering interface-type resistive switching in bifeo3 thin film switches by ti implantation of bottom electrodes. Sci. Rep. 5, 18623. https:// doi. org/ 10. 1038/ srep1 8623 (2015). 32. Sze, S. M. & Ng, K. K. Physics of Semiconductor Devices (Wiley, 2007). 33. Tyagi, M. S. Physics of Schottky Barrier Junctions 1–60 (Springer, 1984). 34. Bruce, P. G. Solid State Electrochemistry. Chemistry of Solid State Materials (Cambridge University Press, 1994). 35. Meyer, R. et al. Oxide dual-layer memory element for scalable non-volatile cross-point memory technology. In Proceedings—9th Annual Non-Volatile Memory Technology Symposium, NVMTS. https:// doi. org/ 10. 1109/ NVMT. 2008. 47311 94 (2008). 36. Sun, B. et al. Non-zero-crossing current-voltage hysteresis behavior in memristive system. Mater. Today Adv. 6, 100056. https:// doi. org/ 10. 1016/j. mtadv. 2020. 100056 (2020). 37. Shuai, Y. et al. Improved retention of nonvolatile bipolar bifeo3 resistive memories validated by memristance measurements. Phys. Status Solidi C 10, 636–639. https:// doi. org/ 10. 1002/ pssc. 20120 0881 (2013). 38. Shuai, Y. et al. Coexistence of memristive and memcapacitive effects in oxide thin films. Jpn. J. Appl. Phys. 57, 121502. https:// doi. org/ 10. 7567/ jjap. 57. 121502 (2018). 39. Hudait, M. & Krupanidhi, S. Doping dependence of the barrier height and ideality factor of au/n-gaas schottky diodes at low temperatures. Phys. B 307, 125–137. https:// doi. org/ 10. 1016/ S0921- 4526(01) 00631-7 (2001). 40. Du, N. et al. Synaptic plasticity in memristive artificial synapses and their robustness against noisy inputs. Front. Neurosci.https:// doi. org/ 10. 3389/ fnins. 2021. 660894 (2021). 41. Schmidt, H. et al. Integrated non-volatile memory elements, design and use. US Patent 9520445 (2016). 42. Jin, L. et al. Transport properties of ar+ irradiated resistive switching bifeo3 thin films. Appl. Surf. Sci. 336, 354–358. https:// doi. org/ 10. 1016/j. apsusc. 2014. 12. 136 (2015) (E-MRS 2014 Spring Meeting. Symposium J. Laser Interaction with Advanced Materials: Fundamentals and Applications). Acknowledgements This work was supported by the DFG (German Research Foundation) Priority Program Nano Security, Projects MU 2332/10-1, DU 1896/2–1, PO 1220/15–1. N.D. acknowledges the Young Scientist Project in the Priority Program Nano Security. Y.S., T.H. and T.M. acknowledge the DFG in the frame of Research Grant SFB 1461 (Project-ID 434434223). We acknowledge Dr. Danilo Bürger and Prof. Dr. Heidemarie Schmidt for preparing and providing the physical BFO memristors. We acknowledge support by the Open Access Publication Funds of the Ruhr-Universität Bochum. Author contributions S.Y. prepared the data sets, implemented the methodology, conducted the simulations, analysed the results, and prepared the first draft of the manuscript. T.M., N.D. and I.P. conceived and directed the conceptual ideas of the work. S. Y., T.H. and T.M. developed the methodological concept for memristor modelling. N.D. designed the experimental work. X.Z. and Z.C. carried out the electrical experiments. T.M. acquired the funding, and administered the project. All authors contributed with interpreting and discussing the results as well as manu- script writing and revision. https://doi.org/10.1109/HST.2017.7951797 https://doi.org/10.1038/nmat2023 https://doi.org/10.1109/IVSW.2019.8854394 https://doi.org/10.1038/srep02208 https://doi.org/10.1103/PhysRevApplied.10.054025 https://doi.org/10.1002/adfm.201303365 https://doi.org/10.1038/srep13753 https://doi.org/10.1038/srep35686 https://doi.org/10.1021/acsami.7b19836 https://doi.org/10.1109/SISPAD.2016.7605167 https://doi.org/10.1088/1361-6463/aa69ae https://doi.org/10.1109/TED.2014.2325531 https://doi.org/10.1109/TCSI.2020.3018502 https://doi.org/10.1109/43.541446 https://doi.org/10.1063/5.0084085 https://doi.org/10.1063/5.0084085 https://doi.org/10.1126/science.1168636 https://doi.org/10.1063/1.3589814 https://doi.org/10.1007/s10854-016-4784-y https://doi.org/10.1021/am504871g https://doi.org/10.1038/srep18623 https://doi.org/10.1109/NVMT.2008.4731194 https://doi.org/10.1016/j.mtadv.2020.100056 https://doi.org/10.1016/j.mtadv.2020.100056 https://doi.org/10.1002/pssc.201200881 https://doi.org/10.7567/jjap.57.121502 https://doi.org/10.7567/jjap.57.121502 https://doi.org/10.1016/S0921-4526(01)00631-7 https://doi.org/10.3389/fnins.2021.660894 https://doi.org/10.3389/fnins.2021.660894 https://doi.org/10.1016/j.apsusc.2014.12.136 https://doi.org/10.1016/j.apsusc.2014.12.136 13 Vol.:(0123456789) Scientific Reports | (2022) 12:20490 | https://doi.org/10.1038/s41598-022-24439-4 www.nature.com/scientificreports/ Funding Open Access funding enabled and organized by Projekt DEAL. Competing interests The authors declare no competing interests. Additional information Correspondence and requests for materials should be addressed to S.Y. or N.D. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. © The Author(s) 2022 www.nature.com/reprints http://creativecommons.org/licenses/by/4.0/ Physics inspired compact modelling of  based memristors Simulation approach Results and discussion Conclusion Methods Experiments. Simulations. References Acknowledgements