Multiple slopes in the negative differential resistance region of NbOx-based threshold switches

Niobium oxide devices exhibit threshold switching behavior which enables their use as selectors in memory arrays or as locally active devices for neuromorphic computing. Among the basic dynamical phenomena appearing in non-linear circuits, the oscillations generated in a relaxation oscillator, which is making use of the negative differential resistance (NDR) effect of a threshold switching device, are of special significance for the design of neuromorphic electronic systems. Here, the necessary requirements for the emergence of oscillations of this kind in a simple relaxation oscillator circuit and their influence on the shape of the measured quasi-static I– characteristic of the threshold switch are examined. In the corresponding experiments multiple NDR regions were found to appear in the quasi-static I– characteristic of the threshold switch concurrently with the occurrence of oscillations. The observed ‘multiple NDR phenomenon’ is therefore merely a measurement artefact due to the averaging effect associated to the operating principles of the source measure unit (SMU) utilized to measure the device current and voltage. In this work, we analyzed how the emergence of oscillatory behavior in the relaxation oscillator depends upon the device layer stack composition. The probability of the appearance of oscillations within a large current range can be increased by decreasing the oxygen content in the sub-stoichiometric bottom layer of a niobium oxide bi-layer stack. It is shown that this trend is caused by the resulting decrease in the value of the product between thermal capacitance and thermal resistance of the threshold switching device. Furthermore, the changed stack composition reduces the variability and changes the forming voltage, which goes hand in hand with a change of the threshold voltage.


Introduction
Threshold switching in niobium oxide-based devices is a well-known phenomenon that was already discovered in the 1960s [1]. The threshold switching effect typically reveals itself through the appearance of a hysteresis loop in the quasistatic I-V m characteristic of the device [2]. Recently it was shown that the threshold effect originates from the existence of an inherently unstable NDR region in the aforementioned I-V m characteristic. The NDR region can be stabilized by connecting a suitable resistor in series to the device [3]. However, in the first publications, the physical mechanism behind the threshold switching effect in niobium oxide devices has remained unclear and still is the subject of continuous discussions in the scientific community. Nowadays there are only two different widely accepted theories. One theory, based on a metal-insulator-transition (MIT), assumes a transition from the insulating to the conductive phase due to a strong heating inside a previously formed filament [4,5]. Such a MIT would need a very high local temperature of about 1070 K. However, in the analysis of temperature-dependent measurements, a maximum temperature inside the filament of about 550 K was estimated, contradicting the MIT assumption [3]. Therefore, it is more reasonable that the underlying physical mechanism is based upon a temperature-activated Frenkel-Poole conduction phenomenon [3,[6][7][8].
The occurrence of multiple NDR regions in the quasistatic I-V m characteristic of a threshold switch was observed in several publications [9,10]. The temperature-activated Frenkel-Poole conduction theory, at the basis of our model describing the mechanisms behind the emergence of threshold switching in niobium oxide nano-devices, struggles to explain this phenom enon, which the authors from [9] attributed to a progressive MIT with growing metallic nano-domains.
As a thorough understanding of the physical mechanisms underlying the threshold switch locally active behaviour [2] is necessary to enable proper modeling and prediction of the electrical behavior of the device for its adoption in complex circuits, it is necessary to clarify the origin of those multiple NDR regions.
The feasibility of an oscillator circuit based on niobium oxide threshold switching devices was already proven in one of the very first papers dealing with threshold switching in niobium oxide [10]. In 2012 oscillator circuits of this kind received a lot of attention when their oscillations were connected to the dynamics in neuromorphic circuits [4]. An indepth understanding of the necessary requirements to induce the emergence of the desired oscillatory behavior in circuits based upon NbO x threshold switches is crucial for the optimization of the sample stack and paves the way to the future development of large mem-computing or neuromorphic hardware platforms.

Methods
To build oscillator circuits, niobium oxide-based threshold switching devices were fabricated, mounted into packages and connected electrically by wire bonding. For this purpose, threshold switching devices with an area below 1 µm 2 were structured applying electron beam lithography and connected to larger electrodes structured by photolithography. In a first study, aimed at gaining a good understanding of the correlation between dynamic behavior and shape of the quasi-static I-V m characteristics (see section 3.1) each threshold switching device consists of a two layer stack with a niobium bottom layer covered by an oxygen-rich niobium oxide layer deposited by reactive sputtering. Those niobium and niobium oxide layers were sandwiched between two sputtered electrodes consisting of a Ti (5 nm)/Pt (35 nm) stack each.
The oscillator circuit depicted in figure 1 was constructed using the packaged threshold switching device, a resistor R ser with 5540 Ω resistance, and different parallel capacitors. The contact resistance R con of the memristor was assumed to be 30 Ω. The oscillatory behavior of this circuit was then investigated through electrical measurements. In order to generate oscillations in the circuit shown in figure 1, a constant voltage of 18 V was applied to node 1 and the oscillating voltage at node 2 was recorded using an oscilloscope. With regards to the experimental determination of the quasi-static I-V m characteristics, instead of a constant voltage, a voltage sweep was applied to node 1 in the circuit of figure 1, and the current I was measured using the same source measure unit (SMU) adopted to generate the voltage sweep. The voltage drop V m over the memristor was measured at node 2 in the oscillator circuit by performing a four-terminal measurement in order to cancel out the effects of the electrode and contact resistances. Throughout the manuscript, for simplicity, we shall refer to I and V m as memristor current and memristor voltage, respectively.
To investigate the influence of the stack composition on the appearance of multiple NDR regions in the quasi-static I-V m characteristics, various niobium oxide bi-layer stacks were sandwiched between two platinum electrodes. In this second study (see section 3.2) each bottom electrode consists of an unstructured evaporated 5 nm thick titanium adhesion layer underneath a 20 nm thick evaporated platinum layer. As top electrodes, evaporated 40 nm thick platinum dots with a diameter of 50 µm were used. Further, here, the niobium oxide stack consists of two layers, which were both reactively sputtered from a niobium target. Particularly, the bottom layer was sputtered with oxygen mass flows of 0 sccm (acronym for standard cubic centimeters per minute), 0.5 sccm, and 1 sccm corresponding to 0%, 3.3%, and 6.7% of oxygen inside the argon atmosphere. The thicknesses of those layers are approximately 4.5 nm. On top of those bottom layers 10 nm thick oxygen-rich niobium oxide layers were deposited using 4 sccm oxygen mass flow, corresponding to 26.7% of oxygen inside the argon atmosphere. In order to create the conducting filament an electro-forming step was preliminarily performed by sweeping the applied voltage V s from 0 V to 10 V while using a current compliance of 1 mA. Afterwards the I-V m characteristics were measured by applying a voltage sweep between 0 and 9 V, and setting the current compliance to 0.7 mA. Throughout all measurements in this second study, including the forming step, a series resistance of 9950 Ω was used inside a experimental circuit setup comaparable to the one depicted in figure 1. The only parallel capacitances during the measurements were the parasitic cable capacitance of the probe station, amounting to around 300 pF, and a parasitic capacitance of the metal-insulator-metal structure itself, which was found to vary negligibly from sample to sample as compared to the overall parasitic capacitance, having an average value of 67 pF. During the measurements, the SMU used for the voltage sweep was also used to determine the current I (see figure 1) while the voltage drop V m over the memristor was measured separately by a second SMU.

Correlations between dynamic behavior of the threshold switching device and shape of its quasi-static current-volt age characteristics
In this section the dynamics and the quasi-static behavior of NbO x -based threshold switching devices are studied and the necessary requirements for the emergence of oscillations in the circuit depicted in figure 1 are determined. The quasistatic I-V m characteristics of the considered memristor were modelled by assuming a temperature-activated Frenkel-Poole conduction in the filament and an additional parasitic conduction path in parallel to the filament. This parasitic conduction path features a non-linear dependence on the applied electric field, but, in contrast to the conduction mechanism inside the filament, it does not show any activation due to elevated temper atures. Such a parasitic conduction mechanism could be described either by hopping conduction or by Frenkel-Poole conduction [11]. Since the conduction mechansim inside the parasitic current path has no direct influence on the threshold switching behavior, but only leads to an increased leakage current below the threshold voltage, we did not investigate it in more detail, and assumed a Frenkel-Poole conduction mechanism for the fitting procedure. The parasitic conduction path might be caused either by some niobium oxide film inhomogeneity at the edges of the bottom electrode, that could result from the topography due to the lift-off process used to structure the bottom electrode, or by some partially formed filaments created during the electroforming process. The parasitic conduction path has a less confined geometry, explaining why its temperature can be assumed to stay constant at a value equal to the ambient temperature T amb . We thus model the parasitic conduction path similarly as in [3] assuming its variable resistance R p changes with the memristor voltage according to where R 0 is a constant resistance, E a describes the activation energy, q is the elementary charge, V m is the voltage across the memristor, t ox is the effective oxide thickness, ε 0 and ε r are the vacuum and relative permittivities, and k is the Boltzmann constant. A slightly modified formula, as compared to the mathematical expression published in [11], was used to describe the variable resistance in the filament: In this formula the barrier lowering effect caused by a large donor state concentration N D as well as the fitting parameter λ, which describes the de-centering of the maximum of the potential well between two donor states due to the impact of the electric field, were included. The resistance R 0 depends on the effective oxide thickness, the density of states N C in the conduction band, the electron mobility µ and the diameter of the filament d fil according to ( The rate of change Ṫ of the temperature inside the memristor filament is governed by Joule heating effects and can be expressed bẏ The internal heating, resulting from the power dissipation inside the filament, depends on the thermal capacitance C th . Additionally the thermal conductance Γ th plays an important role as soon as a significant temperature difference ∆T between the filament temperature T and the ambient temperature T amb arises. The fitting parameters for the parasitic conduction path as well as for the 'virgin' memristor conduction path are given in table 1. The simulation results in comparison to the measured quasi-static I-V m characteristics are depicted in figure 2(a).
With regard to the memristor filament, a very large donor state concentration as well as a small relative permittivity were obtained from the fitting procedure. This suggests that, after the electroforming step is completed, the filament contains a high concentration of oxygen vacancies, which act as donor states, leading to a decrease in the relative permittivity [12]. Although the fit shows that the donor state concentration in the filament is high whereas its dielectric constant is small as compared to the respective parameters of stochiometric Nb 2 O 5 , those values do not represent absolute results, since a change in ε r could for example be compensated by a corresponding change in N D , while none of these parameters could be determined by independent experiments. The Memristor-based oscillator circuit using a series resistance R ser , a memristor M with a contact resistance R con and a parallel capacitance C par .
fitting results for the parasitic conduction path show a larger value for d fil . This is expected since the parasitc path does not correspond to one clearly defined filament. In contrast to the parasitic conduction path, where the temperature remains constant at room temperature, the temperature inside the filament increases significantly during the voltage sweep (see figure 2(b)). Nevertheless, according to the results obtained by simulations, a temperature of only 420 K is reached at the threshold voltage, confirming that the filament temperature is still well below the temperature necessary for a MIT when the threshold switching sets in. Due to the constant temperature T, the influence of the parasitic conduction path on the current regime of the NDR region becomes negligible.
The thermal capacitance cannot be determined by fitting the quasi-static I-V m characteristics, because it only influences the dynamic properties of the device. These dynamic properties will manifest themselves when applying the device in the relaxation oscillator circuit as shown in figure 1.
Overall the oscillator circuit consists of a constant voltage source, coupled in series with a resistance R ser , and the memristor M with its parallel parasitic capacitance C par . The voltage at node 2 will start to oscillate if the memristor is biased in the locally-active NDR region of its DC I-V m loci [2]. Under this condition the capacitor will be progressively charged, and, as soon as the voltage V m exceeds the threshold voltage of the memristor, the resistance of the latter will decrease initiating the discharging process of the capacitor. This discharging continues with the power dissipation in the filament getting smaller at low voltages until the internal temper ature decreases so much that the memristor switches back to its high resistance state. Then the next oscillation cycle starts by charging the capacitor again. To enable an experimental investigation of the oscillatory behavior without the influence of parasitic cable capacitances from the measurement setup, the NbO x -based threshold switching devices were mounted into DIL-packages using wire bonding, and the oscillator circuit was realized on a breadboard. In that way the parasitic parallel capacitance could be reduced to about 10 pF. To bias the memristor in its locally active NDR region and thus to enable oscillations, a series resistor R ser of 5.5 kΩ was chosen, and an external DC voltage V s of 18 V was applied to node 1 (see figure 1). The series resistance was chosen so large as to ensure the stabilization of the entire NDR region according to the small value of the slope of the load line of the resistor [2], allowing a precise determination of its shape. When no parallel capacitor was used, only negligible instable oscillations could be experimentally observed (see figure 3(a)). Furthermore, a minimum parallel electrical capacitance is necessary to ensure that the product between electrical resist ance and electrical capacitance in the oscillator circuit (R ser C par ) is significantly larger than the product between thermal capacitance and thermal resistance of the threshold switching device (R th C th with R th = 1/Γ th ) [13,14]. If the latter condition is not fulfilled the internal heating and cooling rate of the filament at a given Joule heating is not Table 1. Values for the coefficients used to fit the I-V m characteristics of the threshold switching device resulting from a superimposition of the characteristics of the memristor and the characteristics of the parasitic current path. The symbol for the fitting parameter d fil (see equation (3)) is denoted as d filp and d film for the parasitic and 'virgin' memristor current paths, respectively.

Coefficient
Value large enough to follow the oscillation frequency that is given by the inverse of the electrical R ser C par product. Thus, the dynamic range of the temperature T inside the filament of the thermally-activated memristor device is not large enough for the device to leave the equilibrium state and therefore no oscillation emerges [14]. To achieve stable oscillations the parallel electrical capacitance was varied from 220 to 820 pF and the voltage oscillations at node 2 in figure 1 were recorded (see figures 3(b)-(e)), whereby the oscillation frequency decreases with an increase in the parallel capacitance C par . The value for C th in the fit of the measured device characteristics was adjusted to match the dynamic behavior that was exper imentally observed showing first oscillations with a parallel capacitance of 220 pF. Using the relaxation oscillator circuit the quasistatic I-V m characteristics of the threshold switch were measured with various parallel capacitances, utilizing only one SMU to apply a voltage sweep and to measure the current at the same time, whereby the minimum integration time of the SMU is 20 ms. As a result, the measured value of the cur rent plotted in figure 4 equals the current through the series resistor (at node 1 in figure 1). The resulting I-V m characteristics depicted in figure 4 show a gradual transition from a single NDR region towards multiple NDR regions with different slopes. In case of the single NDR region, which is present for small parallel capacitances, the voltage over the threshold switch decreases gradually with increasing current within the NDR region (see figure 4(a)). For the case without an external parallel capacitor we assume the parasitic capacitances of the setup and the device to be negligible small, and, therefore, in the following we will refer to this case as the scenario without a parallel capacitor. When the parallel capacitance is increased to a value that causes the emergence of oscillations (see figure 4(b)) four different regions with negative differential resistance (NDR) arise. The first region (visualized by red background coloring) has a very small slope and the voltage over the threshold switch decreases while the current stays almost constant. In the second region (green) the slope is even larger than the slope observed in case of a single NDR region. The third region (orange) is similar to the first (red) one where the voltage decreases significantly while the current stays almost constant. After the third NDR region the curve shows the same behavior as the characteristic measured without a parallel capacitance. Region 4, visualized by blue background coloring, ends as soon as the slope of the quasi-static I-V m characteristic turns back to a positive value (V h ). For an increasing parallel capacitance up to 820 pF the first (red) and the third (orange) NDR regions contain less measurement points, while the size of the second NDR region increases. Further, the size of the fourth NDR region decreases (see figures 4(c)-(e)). If one would increase the parallel electrical capacitance even further beyond 820 pF, the second NDR region would expand until the fourth NDR region completely disappears, while the slope of the third NDR region directly changes into a positive slope. According to these observations it seems logical to assume that the oscillations arising during the measurements with large parallel capacitances correlate with the changes in the shape of the NDR region. To provide evidence for this assumption, the time-dependent voltage course during the linear voltage sweep was simulated for the different parallel capacitances used in the measurements (see figure 5) applying LTspice transient simulations [14,15]. Analyzing the resulting time-dependent voltage course reveals that the time span over which oscillations occur increases with increasing parallel electrical capacitance. The oscillation frequencies observed in the simulated curves range from 2 MHz to 6 MHz which is in a similar range compared to the measurements depicted in figure 3. Since a linear voltage sweep was applied during the measurements, the increased time window for oscillatory behavior corresponds to a larger range of V s in which oscillations occur. From the simulated oscillations, Figure 3. Oscillations of the voltage at node 2 in the circuit shown in Figure 1 under an applied bias voltage V s of 18 V, R ser = 5540 Ω and a variety of capacitors connected in parallel with the memristor: C par = 0 pF, C par = 220 pF, C par = 330 pF, C par = 470 pF, and C par = 820 pF.  figure 4 can therefore be compared with the simulation results as shown in figure 6. This comparison shows a good agreement concerning the development of a sharp transition into the second NDR region as the size of the parallel capacitance gets larger. The measurement without a parallel capacitor, which has been considered as a reference to the one taken for the model parameter determination, is included as a blue curve in figures 6(b)-(e).   The increase in the current with increasing V s deviates only slightly from a linear behavior even for a large parallel capacitor (see figure 7). This leads to an almost direct transfer of the curve shape observed in figure 6 into the I-V m characteristics, shown in figure 8. The comparison between the simulation and the measurement results of the I-V m characteristics of the threshold switch, depicted in figures 8(a)-(e), shows a good qualitative agreement in terms of the basic shapes of the four possible NDR regions for each of the parallel capacitances under study. The slight shift of the NDR regions towards lower currents in the measurement results, observed for large parallel capacitances, is most likely due to an undesired nonvolatile switching effect occuring during the course of the measurements [16]. Overall, on the basis of the simulated time-dependent course of the average voltage drop over the memristor, it was possible to show that the occurence of multiple NDR regions in measurements which were performed using a SMU is caused by oscillations during the measurements. Therefore, our findings show that the occurrence of multiple NDR regions in the quasi-static I-V m characteristic of the device can be in accordance with a temperature-activated Frenkel-Poole conduction as the physical mechanism responsible for the threshold switching effect. Up to now only the influence of the parallel capacitance on the shape of the NDR region in the measured quasi-static I-V m characteristic was discussed, while the memristor parameters were kept fixed to the values given in table 1. As discussed before, oscillations occur provided the condition R ser C par R th C th is fulfilled. This means that the shape of the NDR region can be influenced also by keeping C par and R ser constant while changing R th or C th . A variation of R th influences the shape of the whole I-V m characteristic making the study of the shape of the NDR region harder. Therefore, only the influence of the thermal capacitance will be investigated in the following. The influence of the thermal capacitance on the shape of the NDR region of the quasi-static I-V m characteristics is displayed in figure 9 by using the previously applied memristor model and the oscillator circuit depicted in figure 1, while varying C th at a fixed value for C par of 470 pF. When analyzing the simulation results that are depicted in figure 9 it becomes obvious that the influence of a decreasing C th on the NDR shape is comparable to the influence of an increasing C par (refer to figure 4). For a large value of C th , equal to or larger than 45 fJ K −1 , no oscillations occur and the resulting I-V m characteristics exhibit only a single smooth NDR region. When decreasing C th to 35 fJ K −1 the V s -range in which oscillations occur gets sufficiently large to lead to the emergence of multiple NDR regions with different slopes. These multiple NDR regions get more pronounced when the value of C th is further decreased. This is caused by oscillations occuring in an increasing V s -range. At some point the V s -range in which oscillations occur cannot increase any further, even if the value of C th is further decreased, because it is limited by the hold voltage V h . This effect becomes obvious when comparing the resulting I-V m characteristics for values of C th of 0.001 fJ K −1 and 0.000 01 fJ K −1 , looking similar one to the other despite the fact that one of the C th values is two orders of magnitude smaller than the other one. According to the conclusions drawn from figure 9, it can be assumed that the analysis of differences in the shape of the NDR region in the quasi-static I-V m characteristics of the devices measured with the same experimental setup (with fixed C par and R ser ) can give insight into the underlying sample-to-sample differences in C th and/or R th . This finding will be adopted in the next section to analyze the thermal properties of threshold switching devices having different stack compositions.

Impact of the sample composition on the appearance of multiple NDR regions
The influence of a variation in the composition of a bi-layer stack, consisting of a substoichiometric niobium oxide layer underneath an oxygen-rich niobium oxide layer on the shape of the NDR region in the quasi-static I-V m characteristic was analyzed systematically. As discussed in detail in the previous  Here, I denotes the current through the series resistor. The capacitance of the capacitor in parallel to the memristor was kept at a fixed value of C par = 470 pF and the thermal capacitance of the memristor was varied among the following values: C th = 45 fJ K −1 , C th = 35 fJ K −1 , C th = 30 fJ K −1 , C th = 25 fJ K −1 , C th = 20 fJ K −1 , C th = 15 fJ K −1 , C th = 10 fJ K −1 , C th = 5 fJ K −1 , C th = 0.1 fJ K −1 , C th = 0.001 fJ K −1 and C th = 0.000 01 fJ K −1 . chapter, variations in the shape of the NDR region are caused by a distinct oscillatory behavior, which is most likely caused by a difference in the thermal capacitance and/or in the thermal conductance of the memristor. The device stack composition was modified by changing the oxygen flow (the percent age of oxygen in the mixture of oxygen and argon) during the reactive sputtering of the substoichiometric bottom layer between 0 sccm, 0.5 sccm (3.3%) and 1.5 sccm (6.7%). For each of the three different oxygen contents in the bottom layer the behavior of 20 devices has been determined in measurements. The measurements were performed through a probe station with a cable capacitance of around 300 pF acting as the parasitic parallel capacitance and a resistance of 10 kΩ connected in series with the threshold switch in order to stabilize the NDR region. As it was concluded from figures 4 and 9, the occurrence of multiple NDR regions in the quasi-static I-V m characteristic of the device may be attributed to oscillations in the NDR region during the measurement for the given series resistance and parallel capacitance, and depends most likely on changes in the product between thermal capacitance and thermal resistance. Based on this knowledge, the measured quasi-static I-V m characteristics of the 60 devices with the three distinct stack compositions considered in this phase of the study, were divided into three categories, according to the shape of their NDR region, as shown in figure 10. One group consists of characteristics that clearly exhibit multiple NDR regions, like the red curve in figure 10, corresponding to a device with oscillations in a large V s -range. Whereas characteristics with only slight deviations from the ideal shape with a single NDR region due to oscillations in a small V s -range, like the blue curve in figure 10, were classified in the group with intermediate NDR region form. The third group consists of curves that display only one smooth NDR region as one would expect from our model for the quasi-static I-V m characteristic (see for example the black curve in figure 10).
To determine the influence of the stack composition on the oscillatory behavior, we analyzed how many of the I-V m characteristics pertaining to the 60 devices under focus fall into each of the three different categories. The resulting histogram, shown in figure 11, reveals that a lower oxygen content in the bottom layer leads to an increased number of curves displaying multiple NDR regions. This effect can be explained by the reduction in thermal capacitance and the increase in thermal conductance resulting from the use of a more metallic instead of a substoichiometric oxide bottom layer. The lower thermal resistance-capacitance product R th C th leads to oscillations already for a smaller electrical resistance-capacitance product R ser C par , and, therefore, also enables oscillations with larger frequencies. For an increasing oxygen flow during sputtering of the bottom layer the amount of devices exhibiting intermediate NDR region form or a single NDR region increases, pointing towards an increase in R th C th . The most important parameters for the operation of the threshold switching devices are the threshold voltage, and the threshold current. The hold voltage is an important device parameter as well because it provides the upper limit of the possible oscillation range. Since in a quasi-static measurement the large currents required to determine the hold voltage would lead to severe device degradation, this parameter will not be evaluated in the following. The characteristic parameters are marked in figure 12, which depicts an example measurement of the quasi-static I-V m characteristics of a device during and after the electroforming step. In figure 13 the statistical distributions of forming voltage, threshold voltage and threshold current for the three different stack compositions under investigation   are summarized. Only a small variability was observed in the forming voltage of the devices with a pure niobium bottom layer, while for the sample where the bottom layer was sputtered with 1 sccm oxygen flow the variability was significantly larger. The sample sputtered with only 0.5 sccm oxygen flow exhibited the largest variability in the forming voltage. The smallest variability observed in the forming voltage of the 20 devices with a niobium bottom layer could be partially attributed to a reduction in the median value of the forming voltage, caused by a chemical reaction between the niobium bottom layer and the niobium oxide top layer, leading to a smaller effective niobium oxide thickness. However, this effect cannot explain the huge variation in the forming voltage for the devices whose bottom layer was sputtered with 0.5 sccm oxygen flow. One possible explanation for this other effect is that for the devices, sputtered with low oxygen flow, the distribution of oxygen in the bottom layer could be inhomogeneous, leading to local variations in the forming voltage as well as in the thermal properties. The inhomogeneity in the oxygen distribution in the bottom layer might become more pronounced for very low oxygen flows (e.g. 0.5 sccm) during sputtering. The same trend in variability is also apparent when comparing the distribution of the threshold voltage for the different samples. The smallest variability is again observed for the sample with a niobium bottom layer. The considerably smaller variability in the threshold voltage of these devices can be attributed to their reduced forming voltage. In fact it can be assumed that a more uniform and gentle forming leads naturally to a closer resemblance between the shapes of the filaments formed within the distinct devices, with a consequent reduction in the variability of their threshold switching properties. However, with regards to the variability in the threshold current, no significant difference was observed between the three different stack compositions.
All in all, the comparison between the three stack compositions with different oxygen contents in the NbO x -based bottom layer indicates that a reduction in the forming voltage can be attained by using a purely metallic bottom layer. This further results in better control and narrower distribution of the forming and threshold voltages. However, despite the fact that there are clear indications that both the values and the variability in forming and threshold voltages are affected by the thermal capacitance and by the thermal conductance of the device, another factor influencing them is the effective oxide thickness. From the data displayed in figure 13, it is not possible to distinguish between these two effects. In order to separate the influence of the two different factors, the measured quasi-static I-V m characteristics of the 60 measured devices were classified in the three categories independently from the device stack composition: curves featuring multiple clearly visible NDR regions (first category), curves with an intermediate NDR shape, differing only slightly from a single NDR region (second category), and curves exhibiting a single NDR region (third category). Since the numbers of devices with I-V m characteristics falling in the distinct categories are different, a device variability study would not be possible. Therefore, the  distribution of the aforementioned characteristic parameters was analyzed for 15 devices per category (see figure 14). The variability in the forming voltage increases from devices out of the first category, featuring I-V m characteristics with multiple NDR regions, through devices out of the second category towards devices out of the third category, whose I-V m curves feature a single NDR region. Additionally, also the median value of the forming voltage increases from 1.25 V for devices in the first category through 1.64 V for devices in the second category to 2.15 V for devices in the third category. As for the threshold voltage, a similar trend can be observed. The variability in this parameter increases from devices in the first category through devices in the second category to devices in the third category. Similarly, the median value of the threshold voltage increases from 1.14 V for devices in the first category through 1.21 V for devices in the second category to 1.36 V for devices in the third category. In case of the threshold current, the variability does not significantly change for the three different categories but the median value is increasing from 315 µA for devices in the first category through 328 µA for devices in the second category to 414 µA for devices in the third category. Overall, the classification of the measurement data on the basis of the observed NDR region shape makes trends for all the characteristic device parameters more clearly recognizable. Based on this classification, and through experimental and numerical studies presented in this paper, we found out that there exists a strong correlation between the variability in the characteristic parameters of the device, specifically forming and threshold voltages, and the oscillatory behavior of the device, which is most likely controlled by the thermal device parameters, namely thermal resistance and thermal capacitance.

Conclusions
The investigation of a simple oscillator circuit with a niobium oxide-based threshold switching device confirmed that a minimum value for the product between the parallel electrical capacitance and the series resistance is necessary to achieve stable oscillations. Particularly, the electrical time constant R ser C par has to be significantly larger than the thermal time constant R th C th .
Furthermore, in numerical simulations of our device model, based upon the temperature-activated Frenkel-Poole conduction theory, in those cases in which the values of electrical resistance and electrical capacitance enabled the emergence of oscillatory dynamics in the aforementioned circuit, oscillations were found to appear in the NDR region of the quasi-static I-V m characteristic of the device. However, in the corresponding experiments, due to the averaging effect of the SMU, those oscillations were not observed in the experimental quasi-static I-V m characteristic of the memristor, leading instead to the appearance of multiple NDR regions with their shape depending upon the value of the parallel parasitic capacitance.
This was further confirmed in LTSpice simulations, showing that the average voltage drop over the memristor remains nearly constant over the whole current range in which stable oscillations occur. Importantly, our research study reveals that the 'multiple NDR phenomenon', recently discussed in [9], is a mere measurement artefact.
By varying the oxygen content in the bottom layer of a niobium oxide bi-layer stack, it was possible to influence the oscillatory behavior of the threshold switching device. This observation is most likely caused by a change in the product between thermal resistance and thermal capacitance, resulting from a variation of the oxygen content in the device stack bottom layer. A transition from a sub-stoichiometric niobium oxide bottom layer towards a more metallic niobium oxide bottom layer results in an increasing percentage of I-V m characteristics exhibiting multiple NDR regions caused by a larger range with oscillatory behavior due to a reduction in the thermal time constant R th C th . Furthermore, a transition towards a more metallic bottom layer reduces both the median values and the variability in forming and threshold voltages. This seems to point to a correlation between the reduction in the median value and in the variability in V th and V form and a reduction in the product between thermal capacitance and thermal resistance.
In summary our results show clearly that an optimization of the device stack structure can lead to an improvement of forming voltage, threshold voltage and threshold current. Furthermore it was shown that changes in the aforementioned characteristic device parameters go hand in hand with a variation in the oscillatory behavior over the course of the measurement of the quasi-static I-V m characteristics of the device. Finally, with reference to the dynamics of the circuit depicted in figure 1, optimizations of the device stack can lead to larger oscillation frequencies and may even obviate the need for an additional capacitor in parallel to the memristor.