A Compact Model of Interface-Type Memristors Linking Physical and Device Properties

Memristors are an electronic device whose resistance depends on the voltage history that has been applied to its two terminals. Despite its clear advantage as a computational element, a suitable transport model is lacking for the special class of interface-based memristors. Here, we adapt the widely-used Yakopcic compact model by including transport equations relevant to interface-type memristors. This model is able to reproduce the qualitative behaviour measured upon Nb-doped SrTiO$_3$ memristive devices. Our analysis demonstrates a direct correlation between the devices' characteristic parameters and those of our model. The model can clearly identify the charge transport mechanism in different resistive states thus facilitating evaluation of the relevant parameters pertaining to resistive switching in interface-based memristors. One clear application of our study is its ability to inform the design and fabrication of related memristive devices.


Introduction
One of the greatest challenges faced in the future of computation is in how to extract information from and act upon the ever-increasing quantity of data generated. On one end of the spectrum, the amount of big data being produced is constantly increasing [1] while, on the other, low-power autonomous decision-making is becoming increasingly important in edge computing applications. Creating truly intelligent systems with brain-like functionality has so far been beyond our present-day capabilities; doing so in an energy-efficient manner will be another monumental challenge.
The limitations we face in matching a brain's performance across a variety of tasks stem not only from an algorithmic shortcoming, but also from the comparatively low energy efficiency of our systems. Even our best supercomputers cannot solve many tasks that are natural for us [2]. Nearly all computers are built within the von Neumann framework, which enforces a strict physical and functional division between memory and computation; in most cases this leads to the connection becoming a bottleneck for overall performance, because of the data transfer back and forth during a computation. Given that data creation is accelerating faster than computational capabilities, the von Neumann bottleneck is set to become an ever-greater limiting factor. Relinquishing the von Neumann framework in favour of newer paradigms as neuromorphic computing [3] could pave a way to address both the algorithmic limitations that make our computers unable to be truly intelligent [2] and at the same time enable us to start improving their energy efficiency.
One of the most promising emerging technologies that can be used in starting to emulate the brain is a class of electronic device know as memristor. Memristive devices can switch their electrical resistance between two or more levels through the application of external stimuli. They are typically composed of a metal/insulator/metal (MIM) stack which in which microscopic changes occur when an electrical bias is applied across it, resulting in a macroscopic change in resistance. The existence of multilevel conductance states, together with low switching speeds and energy requirements, enables memristors to be used as a remarkably efficient computational substrate which could quite naturally support brain-inspired algorithms and approaches. Memristors can enable the co-location of memory and computation when used to implement synaptic weights in a non-Von Neumann architecture [4].
Different memristive mechanisms are possible depending on the material composition of the MIM cell with phase-change (PCM), filament formation and rupture (ReRAM), magnetic (MRAM), and ferroelectric processes being the most technologically advanced. Interface-type memristors, where electric field controlled resistive switching results from changes occurring at interfaces, are a lesser studied class of memristors but offers distinct advantage for integration in an architecture. Filamentary memristors typically require forming processes before they can be used, which is typically unfavourable for device performance and integration [5,6]; an attractive feature of interface memristors is that switching behaviour is present in the as-fabricated device. Additionally, they also show reproducible and gradual multi-level resistive switching at room temperature.
It has previously been shown that metal/Nb-doped SrTiO 3 (Nb:STO) devices exhibit resistance dynamics in response to repeated voltage pulses describable by a simple power-law of the form ( , ) = 0 + 1 + [7]. Given that this resistance follows an exponential trajectory, the rate of change slows as its bounds 0/1 are approached; such a soft-bounded device can be productively used to model synapses in an ANN while solving life-long learning and catastrophic forgetting issues [8]. While linear synaptic weight changes are often considered ideal, they lead to a response that is either hard-bounded or unbounded; in the continual-learning setting that all biological cognitive systems operate in, this would actually be a disadvantage as it would lead to hindered performance and memory capacity limitations [9]. The resistive switching in Nb:STO Schottky junctions has been documented in literature [10][11][12][13][14] but so far there have been no models that capture the hysteretic current-voltage response.
Large-scale application of memristors is still in its infancy [15] but, like any other component used in computer-aided integrated circuit design, it is paramount to have reliable models of the devices' electrical behaviour, which enable simulation and prediction of the behaviour of both individual memristors as well of the whole integrated circuit. Consequently, a great deal of effort has gone into modelling of the hysteresis curves of these types of memristors: PCM [16][17][18], ReRAM [19][20][21][22][23], MRAM [24][25][26][27], and ferroelectric [28][29][30][31]. The Yakopcic model is a generalised, compact representation of the I-V characteristics of a memristor with its core assumption being that a memristor can be represented as two resistors in series, with an internal state variable mixing between the two.
Here we extend the widely-used Yakopcic generalised memristor model [32] to incorporate the often neglected charge transport through metal/insulator interfaces, which is especially important to model systems where the resistive switching occurs at the interface. While the Yakopcic model has been successfully applied to reproducing the electrical behaviour of many memristive devices [33][34][35][36][37][38][39][40], it has not been used to replicate interface-type memristors as the Nb-doped SrTiO 3 device we select in this work.
In order to derive a functional description of the resistive switching in interface memristors we modify the electron transmission equations in the Yakopcic model to more closely align with the physical mechanisms relevant for charge transport through Schottky junctions in our Nb:STO device. By fitting this physicallyinformed compact model to data obtained from real devices fabricated at different sizes we show clear trends in the model parameters, which match the variations expected by our understanding of the underlying physical mechanisms. Given the generality of the model, it could be applied to other interface-type memristors in order to provide insight into their charge transport parameters and help uncover unaccounted physical mechanisms.

System Modelled
To model a memristive system, we first have to take account of the relevant underlying transport and switching mechanisms. An important consideration that calls for attention and is often neglected is the role of Figure 1: (a) Schematic representation of devices used for obtaining experimental data. (b) Device structure assumed for modelling, highlighting the interface between the cobalt electrode and the insulating Nb:STO layer. (c) Voltage sweeping sequence used to obtain I-V data. The sweep can be split into four parts by considering the polarity (forward bias, ≥ 0, or reverse bias, < 0) and the resistance state (LRS ℎ 1 , or HRS ℎ 2 ). the metal/insulator interfaces. It can be expected that Schottky barriers form at these interfaces due to the significantly different electronic properties of the contacting materials, which can have a considerable influence on the overall device characteristics. A complete physical device model should include these interfaces to describe the full system. The contact resistance between a metal electrode and oxide is controlled by a space-charge region which forms due to the differences in the Fermi levels of the contacting materials. Typically a Schottky barrier forms at a metal/insulator or metal/semiconductor interface, giving rise to a high contact resistance due to the formation of a space-charge (or depletion) layer in which the mobile carriers are depleted.
The effect of such an interface becomes especially important when considering material systems where interface effects play the dominant role in mediating the resistive switching mechanism. This class of memristive devices consist of Schottky interfaces that form between doped wide-band gap oxide, such as SrTiO 3 , and high work-function metals. These devices show hysteretic current-voltage characteristics with bipolar switching and large resistance windows together with continuously alterable resistance states. They do not require an initial forming step, as is typically the case for filamentary memristors, which is a significant advantage for practical applications. Both the SET and RESET transitions are gradual and highly modifiable enabling analogue switching between different states within a window. This distinguishes them from the more advanced resistive switching devices such as phase-change materials, in which the RESET operation is abrupt and filamentary memristors which tend to have abrupt SET transitions.
Charge transport is controlled by the Schottky interface and resistive switching is accompanied by changes in the effective Schottky barrier height and width. This is mediated by processes that occur when an electric field is applied: (i) the movement of ionic species, such as oxygen vacancies, and/or (ii) the trapping of electronic charges. In Nb:STO, for example, oxygen vacancy migration under applied electric fields has been observed to be an important factor in mediating the resistive switching [41].
There are various ways in which electrons can be transported across a Schottky interface under bias. Here we consider two of the most important mechanisms which are thermionic emission of electrons over the top of the barrier and quantum-mechanical tunnelling through the barrier. In ideal diodes thermionic emission is the dominant mechanism and this tends to be the most important process in forward bias; in reverse bias, however, tunnelling is expected to strongly contribute to the flow of charge.
There are a vast number of mechanisms through which electrons can traverse the barrier by tunnelling. Some of the widely considered ones include direct tunnelling -which is possible if the barrier is sufficiently thin -or -when the barrier is thicker -through thinner regions of the barrier, for example at higher energies (Fowler-Nordheim) or through defect states (Poole-Frenkel). The relation between the current and voltage depends strongly on system-specific mechanism(s). In this work we consider the theory proposed by Simmons describing the current flow through a generalised barrier [42]. This model approximates tunnelling current through an arbitrary barrier shape by way of a hyperbolic sinusoid function and is hence more widely applicable as it does not assume one specific mechanism.
Here we focus on the interface between Co (work function 5 eV) and Nb:STO (0.1 wt%) (see figure 1(b)). STO is an insulator that becomes an n-type semiconductor when doped with Nb. When Nb:STO is contacted to a metal with a high work function, a Schottky barrier forms at the interface. It has been found that while this doping concentration has no significant effect on the bulk structure and does not induce oxygen vacancies throughout the single-crystalline substrate, it is sufficiently large to make the Nb:STO degenerate. This ensures that the resistance of the semiconductor is relatively low and the Schottky interface contributes the most to the measured resistance.

Measurements
We obtained experimental results from circular cobalt (Co) electrodes of different areas, contacted to an Nb:STO single crystalline substrate with a doping density of 0.1 wt%. The Co contacts were fabricated using a two-step electron beam lithography process using aluminium oxide as an insulation layer to define the contact areas and to prevent electronic cross talk. The bottom of the substrate serves as a back contact for the devices. A device schematic highlighting the different contact areas is shown in figure 1(a).
were done from the virgin state; the voltage sweep is shown in figure 1(c). For this study, we measured devices of three different areas and radial dimensions of: 10 µm, 32 µm, and 100 µm. To incorporate device-to-device variation, we obtained results from several devices of each area; the results from individual devices are shown in supplementary figure S1.

Yakopcic Generalised Memristor Model
One model which has seen widespread use is the Yakopcic memristor model [32], which is a generalised, compact mathematical representation of the relationship that exists between input voltage (V) and measured output current (I) in a memristor. The model follows in the steps of [19], in that its fundamental representation of the memristor is as two resistors in series: one with low resistance 1 and the other of high resistance 2 . As in [43], this model has been generalised [44] so that the resistors are actually electron transmission equations ℎ 1/2 ( ), making the current a non-linear function of the voltage. Given that a memristor's resistance is analogue and programmable, the two electron transmission equations are weighted and mixed by a dynamic internal state variable ∈ [0, 1], as proposed by [45]. The internal state variable encapsulates a set of physical variables that influence the material's resistance. Given the known importance of oxygen vacancies, we consider the movement of these ions to be the most relevant internal state. Hence, we model the change in as non-linear ion motion. Consequently, the I-V curve describing the characteristic relationship between output and input of the device is given by: with ℎ 1 modelling the behaviour in the low-resistance state (LRS) and ℎ 2 in the high-resistance state (HRS).
The dynamics of the internal state variable are based on two functions ( ( )) (Eq. 2) and ( ( )) (Eq. 3); the former implements the threshold voltage effect by virtue of which does not change when the voltage drop is sufficiently small, while the latter models the non-linearity of ion motion whereby the change is slowed when approaching its boundaries: Finally, the state variable dynamics are modelled by integrating the ordinary differential equation (ODE):

Application to Interface Memristors
We expand on the approach presented in [46] by accounting for the fact that in our Nb:STO device different transport mechanisms are dominant in forward and reverse bias [14]. In order to do so, we define the electron transmission equations ℎ 1 and ℎ 2 (from Eq. 1) as: The parameters min, , min, , max, and max, relate to Schottky transmission and are compact representations of the quantities present in the standard Schottky equation [47]: where is the junction area, * the Richardson constant, the electronic charge, the Boltzmann constant, the temperature and the Schottky barrier height. The ideality factor, quantifies the deviation from ideal thermionic emission. min, , min, , max, and max, , on the other hand, are related to tunnelling as described by the generalised Simmons model [42]: where and are fitting parameters and V is the voltage.
Our choice of ℎ 1 and ℎ 2 reproduces the model proposed in [45] in that it models the interaction between tunnelling and thermionic emission processes taking place at the interface between metal and insulator (see figure 1(b)). While the ordering of the equations in each bias is not essential, it must be congruent with parameters 0 and in order to correctly capture the dominance of thermionic emission and tunnelling in forward and reverse bias, respectively. Strictly speaking our model does involve the physical parameters of the material, but captures the relative contribution of thermionic and tunnelling in the transport process within the Nb:STO devices. An overview of all equations and parameters used by our model is shown table S1.

Parameter Fitting
In order to fit the parameters of our model to experimental data , we simulate its response to the same experimental procedure we used to obtain the I-V characterisation for Nb:STO memristors (section 2.2). We do so by integrating Eq. 5 in order to simulate the memristor models' state evolution in response to the voltage sweeps; at the same time, the current through the device is calculated using Eq. 1. With a slight abuse of notation we can say that each parameter set yields a model which will be in state at each simulation timestep ; we can calculate the current ( , ) using Eq. 1 and compare it to the experimentally determined current at time . The objective function used to optimise the fitting is the mean absolute error (MAE): From the MAE we also derive a mean percent error (MPE) that is useful as a normalised measure of fitting quality between different models: This value is normalised with respect to the number of datapoints and consequently any influence from the larger number of reverse bias samples is removed.
We use a two-step automated regression to meaningfully optimise the model so that it can reproduce the currents measured on the Nb:STO devices. The dataset used as basis to derive the models consists of the currents resulting from repeated measurements (meaning a +1/-2 V sweep as detailed in section 2.2 and shown in figure 1(c)) on the memristive devices of different areas (10 µm, 32 µm, and 100 µm area as shown in figure 1(a)). We thus have multiple Time (s) Voltage (V) Current (A) triples for each of the three devices.
From an initial parameter estimate 0 we regress the model parameters upon each available measurement, which are obtained as described in section 2.2. , are omitted from the optimisation and set to zero because interface-type memristors do not present a voltage threshold, given that they exhibit analogue modulation of their resistance [14]. Three averaged models are obtained by taking the mean of each parameter value for each device area. The parameters , , show little to no variation between the three models and are thus excluded from the next regression step by freezing them to the averaged value. The second regression starts from the same initial parameter estimate as the first but ignores , , -in addition to , -for the reasons just mentioned. By grouping by device area the resulting fitted parameters and fitting a Gaussian distribution to each parameter we obtain three models, each one describing the expected I-V relationship for a device fabricated at each area. We can simulate the mismatch in electrical behaviour of different devices by independently sampling parameter sets from the distributions, thus reproducing the device-to-device variation present among the real memristors.

Sensitivity Analysis
Memristors are prone to device mismatch and cycle-to-cycle variation [48]. It is important to understand and characterise the impact that the variations on the underlying physical characteristics have on the electrical behaviour of the devices as circuits and algorithms built upon memristors need to be aware, or at least robust, to these discrepancies in behaviour.
In order to quantify the impact that each parameter has on the overall behaviour of the models we run a sensitivity analysis. In the three models representing a device area, each parameter is increased and decreased until the MPE (Eq. 11) between the current obtained with the model with the varied parameter and the current generated by the original model is 10%. Intuitively: the percentage change in the parameter necessary to see a 10% change in the measured current is in inverse association with the sensitivity of the model to it; the stronger the influence of a parameter, the less it will need to be changed to see a 10% variation in the overall function behaviour of the model. We exclude and from the analysis as we have elected to fix these to zero when fitting the model, as specified in section 2.5. We also arbitrarily choose to stop the procedure when the change needed in a parameter exceeds 100%, as we consider their impact too small to be relevant after that point.

Parameter Fitting and I-V simulations
Each parameter in the three models is characterised by a Gaussian distribution whose mean and spread (SD) are reported in table 1. The parameter sets are obtained using the two-step optimisation procedure described in section 2.5. Figure 2 shows how the simulated I-V curves for the models fitted to the 10 µm, 32 µm, and 100 µm area devices compare to the experimentally measured I-V curves for those memristors. Each model is simulated 100 times with its parameters independently sampled from the Gaussian distribution reported in table 1. The currents collected by measuring the Nb:STO devices are averaged (red curves); the ones resulting from simulating the models display the range of variation upon repeated trials (green curves). The mean absolute error (Eq. 10) and mean percent error (Eq. 11 between the average experimental and simulated I-V curves are presented in table 2.

Relating Model and Device Parameters
In order to start relating model and device parameters we look for trends in the parameters of the average models for the three different areas. While most existing Schottky barrier models do not take device area into consideration, it is known that device edges in finite device areas will affect transport properties. We can gain insight into this rarely considered area variation through our model.
We expect thermionic emission to be mainly present in forward bias while tunnelling is expected to mostly define reverse bias. The resistive switching process occurs under the application of an external electric field due to oxygen vacancy migration and electronic (de)trapping. If the current density is constant, we would expect the current to increase with device area. It has however been observed that in smaller devices the average electric field controlling the charge transport is larger and decreases with increasing device size. This was shown to significantly impact the current-voltage characteristics; notably, smaller devices were shown to have larger resistance windows [49]. We can interpret systematic variations in different model parameters in light of this knowledge and gain more insight into the role of the electric field.
We distinguish two classes of results we can expect from this analysis: (i) trends in model parameters that validate physical mechanisms we are aware of; (ii) correlations in model parameters that let us infer the presence of overlooked physical mechanisms and that will guide future empirical inquiry.  Given that the parameters used in these equations are strongly tied to the physical mechanisms that govern that transport mechanisms as well as the switching mechanism, we can relate the parameters to device parameters as shown in table 3. This relationship will be further discussed in section 4. Inversely proportional to ideality factor (forward bias).
min, , min, Related to tunnelling magnitude (reverse bias). Table 4 reports the results of the sensitivity analysis obtained following the methodology described in section 2.6. Each of the parameters within the average models ("Mean" parameters in table 1) for the three device areas is increased and decreased until a 10 % change in the MPE (Eq. 11) between the original and varied model is observed. The table is organised from most to least impactful parameter in the model, as defined by its average % of variation necessary to see a 10 % change in the MPE across all tested cases. Where a parameter needs to be varied by more than 100 % to see the aforementioned change in MPE we report it as ">100.0", but consider it as 100 % when calculating the averaging as that is the cutoff point for the search.

Sensitivity Analysis
A preliminary observation of the sensitivity analysis in light of table 3 indicates that the model depends most strongly on the reverse bias parameters and -only in second order -on the forward bias.

Discussion
In this work we set out to find a physical correlations and interpretation of the parameters in the widelyused Yakopcic memristor model. Starting from empirical data measured on memristive devices fabricated at different scales we defined models able to capture their I-V hysteretic behaviour. We were able to find clear correlations between the models' parameters fitted to the different devices, enabling us to reason back and identify the connection to the underlying physical mechanisms.

Connection Between Model and Physical Device
The parameters , , , , , , , and relate to the switching properties of a memristive device. and represent the positive and negative voltages ( ) that need to be exceeded in order to induce a change in the internal state variable and -in turn -the resistance, while factors and control how fast this change is. and are the points where the state variable's motion starts being non-linearly dampened, with and controlling the amount of this damping.
There have been previous attempts to connect the variation in these four pairs of parameters to changes on the filamentary physical device being modelled [44]; it was suggested that , , , most probably vary as a function of the number of oxygen deficiencies, while , , , depend strongly on the degree of crystallinity of the dielectric film and of the electrodes. [44] also looked at the effect of variation of the parameters relating to the I-V equation, namely 1, 2, and as defined in the MIM electron transfer equation: 1 and 2 were found to relate to the thickness of the dielectric layer while to the proportion of Ohmic-totunnelling conduction mechanism. Looking at figure 3(a) we can appreciate that parameters and are inversely correlated to device area. This results in a more gradual change in the models' state in reverse bias for larger devices, which therefore presents a less abrupt transition from thermionic emission towards tunnelling. Given that tunnelling probability is inversely proportional to the potential barrier width [42] this likely indicates that the smaller the device, the thinner its potential barrier.

Parameter
is also negatively correlated to devices' area, with this suggesting that electrons are able to tunnel at lower reverse bias voltages in smaller devices. A higher value for leads to a model where the internal state changes faster. These results are consistent with the expected inverse relation between area and average electric field. The dielectric constant of SrTiO 3 , , has an inverse relation with the electric field strength [50,51]. Given that the barrier width ∝ √ , a narrower barrier is expected for smaller devices.
In reverse bias, the interfacial field is further increased; this is accompanied by a further reduction in the barrier width, facilitating tunnelling at lower voltages.
Turning to look at the reverse bias, we can note that the dual parameters to the ones we have just considered i.e., , , and , are constant with respect to device area. This means that the change in internal state in response to SET pulses is independent of device size. This can be expected for thermionic emission as temperature is an important factor governing at which voltage electrons are able to traverse the barrier, and hence no significant differences are to be expected when comparing the transition to thermionic emission for different devices.
In analysing the parameters that control the magnitudes of the current ( max, / , min, / ), we have to take care to consider that differences in device area are expected to give rise to larger currents in larger devices.
As can be gleamed from figure 3(b), min, and min, correlate to device area in forward bias. These two parameters define thermionic emission behaviour, which is the dominant transport mechanism in this regime. If we compare the relevant functional form of thermionic emission as we implement it (Eq. 7, ( ) ≥ 0) to the Schottky emission equation (Eq. 8), we immediately note that min, is proportional to the height of the Schottky barrier ( min, ∝ ) while min, is inversely proportional to the ideality factor ( min, ∝ 1 ). This suggests that smaller devices may have a lower barrier height and a higher ideality factor. Both these factors account for the increased current flow observed; a lower barrier would allow more transport through thermionic emission while the larger ideality factor indicates the presence of current flow through other transport mechanisms.
On the converse, min, and min, control electron tunnelling in reverse bias, which is the dominant transport mechanism in this other regime. As can be discerned from figure 3(b), min, is proportional to device area, which would lead us to expect more current to flow through larger devices. This is exactly what we see in figure 2, with a higher current at −2 V in the larger devices. However, if we normalise the current flow to account for the larger device area, the current density is larger for smaller devices; this is captured by min, (figure 3(a)), which is inversely proportional to the area.

Sensitivity Analysis
From Table 4, it can be seen that the model is most sensitive to changes in min, followed by min, : the parameters related to reverse bias tunnelling. The third most impactful parameter is max, describes tunnelling in forward bias. From the discussions above and experimental results in [41], tunnelling is largely responsible for the additional current in the LRS compared to the HRS, especially in reverse bias where it is the dominant mechanism. This explains why the model is most sensitive to the parameters controlling the tunnelling behaviour. Similarly, the model is also highly sensitive to changes in , a parameter describing the rate of change of transport from thermionic emission to tunnelling in reverse bias, particularly for the smaller devices. This again alludes to the importance of electronic tunnelling to the switching process and indicates that this mechanism may be more strongly present in smaller devices; this is the probable origin of the larger resistance ratios, which can be appreciated by comparing figure 2(a) to (b) and (c).
The MPE (equation 11) -which we use to evaluate the change in the model's behaviour -is a normalised measure w.r.t. the number of datapoints generated; this means that even though more datapoints are generated in RB due to the deeper voltage sweep, there is not a bias towards the parameters impacting the RB behaviour being reported as more impactful.

Potential Applications and Significance
The model we present in this paper is a first step towards a physical model, as we have started to find correlations between model parameters and corresponding device properties. Defining -and then using -a true physical model would allow us to find an even clearer bijective mapping between the abstract and actual parameter spaces. Moreover, our current model cannot account for many important factors that impact real electrical circuit behaviour as -for example -temperature.
Nonetheless, we can extrapolate some guidelines regarding the expected impact that device design choices may have on the electrical behaviour If the parameters are determined for a particular metal-semiconductor combination, the hysteresis resulting from small variations in experimentally controllable parameters can be simulated. For example, if the simulations indicate the memristive properties can be improved by a change in the Schottky barrier height, experimental strategies can be adopted to achieve this. Specifically, the barrier height could be increased by selecting a metal electrode material with a higher work function, while it could be decreased using a lower work function metal or by introducing a thin oxide barrier in between the metal and semiconductor. Therefore, the model could aid in determining an optimal electrode material or interfacial oxide thickness.
The specific quantitative results of our analyses (as those presented in figure 3 and table 4) are relegated to the properties of Nb:STO and to other interface-type memristors with similar charge transport mechanisms.The model equations presented and our methodology, however, can readily be extended and applied to other types of memristors to include relevant effects that may occur at their metal/insulator interfaces and fully describe the system. Upon defining a physically-informed model of a device it is possible to relate its parameters to the those characterising its physical attributes and/or charge transport mechanisms, for example by running a differential analysis as that we present in section 3.2. After observing how the model parameters change as a function of some device property (in our case, device area) and combining this with a theory of the underlying physical mechanisms, one can define a mapping between the model and device parameter spaces, similarly to what we show in table 3. As mentioned in section 3.2, once this mapping is defined it can productively be used to both verify that predictions made from physical theory are correct and also help in coming up with new theories for the empirical behaviour of the device. This analysis could also prove valuable in aiding device fabrication itself. Designing and testing electronic components is expensive and time consuming process that may call for a large experimental effort in repeatedly fabricating and testing devices to explore a large design space. Physical theory sets the groundwork for model definition, but its predictions can't be verified without observing the real device's (virtuous loop A in figure 4) or the model's behaviours (virtuous loop B in figure 4). Any of the discrepancies between the model and the device behaviours, and the theory's predictions can be separately used to positively improve the theory itself; integrating the contribution from both model and device can lead to even greater refinements (virtuous loop A+B in figure 4). On the other hand, both the model's and fabricated device's behaviours could be integrated to not only inform theory, but each other directly. A sufficiently detailed physical model (or one with a strong mapping to the device parameters) could be used to extrapolate physical parameter sets that promise to yield a device with the desired response; one could then only fabricate devices to those specifications. Model-aided fabrication would therefore be able to inform the creation of the devices themselves, in order to aid in smartly fabricating computational components. It would then be ideal to imagine enacting a virtuous closed-loop procedure (virtuous loop C in figure 4) whereby a physical model would guide fabrication of a device to a certain specification. Model predictions could also be extended to improving physical theory through the results measured on the devices fabricated based on the model's predictions (virtuous loop A+C in figure 4). The measurements on the as-realised device could be compared to the predictions originally made by the model and any error between device response and model predictions could guide optimisation of the latter's parameters. All these relationships could and should be integrated in order to iteratively and synergistically improve theory, model, and physical devices (virtuous loop A+B+C in figure 4). An overview of the relations between physical theory, model, and experimental work are shown in figure 4. Figure 4: Depiction of the relations defining virtuous loops between physical theory, model, and device fabrication and experiments: A: Physical theory guides fabrication and experiments and is verified by the outcome of these. B: Physical theory guides modelling and is -to a certain degree -verified by its predictions. C: Models can also guide experimentation and device fabrication when a strong enough link exists between the model and device parameter sets. These virtuous closed loops can be combined (A+B, B+C, A+B+C) in a synergistic manner.

Conclusions
In this paper, we have extended the widely-used Yakopcic memristor model to describe Schottky interfaces, such as those that form at metal-semiconductor/insulator interfaces. Specifically, we consider thermionic emission and tunnelling to be the possible transport mechanisms and the switching is mediated by nonlinear ion motion. We have tested the model on empirical data collected on cobalt/Nb-doped SrTiO 3 memristors of different areas and found that the I-V behaviour and hysteresis loop are well-captured in all cases. The model accurately describes the significant differences in transport in forward and reverse bias and is able to capture correlations between the underlying device physics and the model parameters.
These trends are further validation of the physical mechanisms we believe define the electrical characteristics of the devices and also highlight the presence of overlooked physical mechanisms that will guide future empirical inquiry. The generality of the model makes it well-suited for different kinds of interface-type memristors whereby the variables provide insight into device parameters.