Device fabrication
Memristive devices were fabricated by sandwiching a SiO2 insulator in between a platinum bottom electrode and a silver top electrode. The choice of silver as the top electrode is based on its electrochemical activity which allows dissolution of silver atoms and migration of Ag+ ions at lower voltages compared to other noble metals, while platinum was chosen as the counter-electrode because it is electrochemical inert (details on the selection of materials and device configuration can be found in Supplementary Section 5). The pad structure devices were fabricated on a thermally oxidized silicon substrate starting with the d.c. magnetron sputtering deposition (power, 200 W) of a TiO2 (10 nm) adhesion layer and a platinum (100 nm) bottom electrode. The homogeneous SiO2 film (20 nm) with a purity of 8N was deposited by radiofrequency (RF) magnetron sputtering with a sputtering power of 150 W in a processing gas mixture of 9 sccm argon and 1 sccm oxygen at 150 °C. Note that the choice of the 8 N SiO2 matrix is related to the very low level of impurities, a potential disturbing factor for achieving controlled conductance states. Also, the resulting SiO2 is rather stoichiometric and chemical and physical interactions with silver are not thermodynamically favourable. Following switching layer deposition, feature sizes of 50 × 50 μm2 were patterned by negative photolithography. Then, the Ag (20 nm) active top electrode was deposited by e-beam evaporation with a deposition rate of 0.01 nm s−1, followed by a d.c.-sputtered platinum (50 nm) capping layer. The role of the capping layer is to prevent degradation of the silver active electrode over time as required for long-term use of the device. A standard lift-off process was utilized for the final cleaning of devices, obtaining an Ag/SiO2/Pt cell with a top electrode size of 50 × 50 μm2.
Device modelling
Departing from the experimental observation of well-defined conductance jumps and states, we model the RESET transition (the SET transition is also considered for generality) as a random generation of events related to the destruction of single quantum mode channels with conductance ∼G0. This is a stochastic version of a continuous behavioural compact model47 which has been successfully applied to memristors with different material systems, different switching modes (bipolar, unipolar, complementary and threshold switching) and for the SPICE simulation of neuromorphic circuits. The stochastic version of the model presented here was recently applied to valence change memory devices which show variability, but not quantum conductance jumps48.
The stochastic resistive switching model follows Chua’s approach49 to memristors and is based on two equations, one for the current and one for the internal memory variable. In our case, the memory state variable is the number of conducting channels, nch, each of these channels contributing ∼G0 to the filament conductance. In a naive interpretation, each of these channels can be considered either as ‘atomic chains’ or as ‘quantized quantum transport modes’ in the filament constriction. This is a simple implementation of the Landauer theory for ballistic transport through an atomic-size constriction50. We consider that the SET/RESET transitions occur by successive discrete conductance jumps (events) corresponding to the creation/destruction of single conduction channels. For simplicity, we assume that each switching event increases or decreases the conductance by the same amount. However, this might not be completely realistic because several channels can be created/destroyed at the same time. During the RESET transition, we will consider that each jump is |ΔG| = G0. Given the experimental results, we impose that the first SET event is abrupt so that the device reaches the compliance limit in a single conductance jump. The creation/destruction of single channels will occur at random times during the application of the external electrical signal (voltage/current). For the sake of generality, we limit the number of channels to nmax. This parameter is related to the maximum area of the filament created during electroforming. Under these conditions, the proposed memory equation is:
$$\frac{\text{d}{n}_{\mathrm{ch}}}{\text{d}t}=\frac{{n}_{\max }-{n}_{\mathrm{ch}}}{{\tau }_{{\rm{S}}}}-\frac{{n}_{\mathrm{ch}}}{{\tau }_{{\rm{R}}}}$$
(1)
where the two terms of the right-hand side (RHS) represent the SET and RESET transitions, and τS and τR are the SET and RESET characteristic times, respectively. Because the SET transition resembles the dielectric breakdown process and is strongly accelerated by the electric field, an exponential voltage dependence for τS is assumed:
$${\tau }_{{\rm{S}}}(V\,)={\tau }_{{\rm{S}}0}\exp \left[-{\gamma }_{{\rm{S}}}(V-I{R}_{{\rm{S}}})\right]$$
(2)
where γS is the acceleration factor, τS0 is the time scale prefactor, I is current, V is voltage and RS is the series resistance. On the other hand, consistently with the electropolishing interpretation, the RESET transition is assumed to be controlled by the oxidation/reduction dynamics and/or by the out-diffusion of species to the filament surroundings. Because both processes are strongly accelerated by temperature, we neglect voltage acceleration (as discussed within the electropolishing interpretation) and we only consider the local temperature rise related to the power dissipated in the filament, \(P=I\left(V-I{R}_{{\rm{S}}}\right)\). Assuming an Arrhenius temperature dependence as a first-order approximation, the characteristic RESET time, τR, can be described as:
$${\tau }_{{\rm{R}}}\left(V\,\right)={\tau }_{{\rm{R}}0}\exp \left[\frac{{E}_{{\rm{a}}}}{{K}_{{\rm{B}}}(T+{R}_{\mathrm{TH}}P)}\right]$$
(3)
where τR0 is the RESET scale prefactor, Ea is the activation energy, KB is the Boltzman constant, T is the external temperature and RTH is the thermal resistance. The thermal resistance has been described in the literature in terms of two parallel paths for heat evacuation51. The longitudinal thermal resistance, RL, corresponding to heat transport along the channel (related to the electrical conductivity) and the transverse resistance, RT, associated with heat transport towards the surrounding material. The latter is independent of the filament size to the first order, while RL is inversely proportional to the filament area, represented here by nch, which is proportional to the area. Thus, we can write \({R}_{{\rm{L}}}={K}_{{\rm{L}}}/{n}_{\mathrm{ch}}\), where KL is a constant. The total thermal resistance is given by the parallel combination of RL and RT, so that \({R}_{\mathrm{TH}}=\left({K}_{{\rm{L}}}{R}_{{\rm{T}}}\right)/\left({n}_{\mathrm{ch}}{R}_{{\rm{T}}}+{K}_{{\rm{L}}}\right)\). It is worth remarking that we included only description of thermal dissipation with a phenomenological approach based on macroscopic parameters such as thermal resistances. While in principle quantum thermal effects cannot be ruled out, experimental works pointed out that these effects only become not negligible in the low-temperature regime52, that is, far away from the room temperature conditions of our work.
Because τS has a strong exponential dependence on voltage, it emerges that \({\tau }_{{\rm{S}}}\ll {\tau }_{{\rm{R}}}\) for positive voltages and \({\tau }_{{\rm{S}}}\gg {\tau }_{{\rm{R}}}\) for negative voltages. Because of this, we can separately consider the SET and RESET transitions with two separate differential equations. One for the SET:
$$\frac{{\rm{d}}{n}_{\mathrm{ch}}}{{\rm{d}}t}=\frac{{n}_{\max }-{n}_{\mathrm{ch}}}{{\tau }_{{\rm{S}}}}$$
(4)
And one for the RESET:
$$\frac{{\rm{d}}{n}_{\mathrm{ch}}}{{\rm{d}}t}=-\frac{{n}_{\mathrm{ch}}}{{\tau }_{{\rm{R}}}}$$
(5)
As far as the current is concerned, we have considered:
$$I\left(V\,\right)=\frac{{n}_{\mathrm{ch}}{G}_{0}}{1+{n}_{\mathrm{ch}}{G}_{0}{R}_{{\rm{S}}}}V+{I}_{{\rm{B}}}\,\sinh \left[\eta \left(V-I{R}_{{\rm{S}}}\right)\right]$$
(6)
where η is a shape parameter related to the potential barrier at the constriction when there are no conducting channels. The first term corresponds to the conduction through the nch channels, and the second to the background tunnelling regime, that is, when the filament has a gap. Although the considered voltage dependence of the background current can be discussed, this is not relevant to our work because we focus on situations where there is at least one conducting channel with a conductance which is generally much larger than that of the background. Finally, notice that nch couples the current and memory equations.
For the generation of random events, we follow an ‘on-the-fly’ method. If the number of events (conductance jumps) is n(t), the event generation rate is \(\lambda \left(t\right)={\rm{d}}n\left(t\right)/{\rm{d}}t\). During the SET transition, \({n}_{\mathrm{ch}}=n\left(t\right)\) so that \(\lambda \left(t\right)={\rm{d}}{n}_{\mathrm{ch}}/{\rm{d}}t\), while during RESET \({n}_{\mathrm{ch}}={n}_{\max }-n\left(t\right)\), so that \(\lambda \left(t\right)=-{\rm{d}}{n}_{\mathrm{ch}}/{\rm{d}}t\). Thus, the event generation rates can be obtained from equations (4) and (5) so that \({\lambda }_{{\rm{S}}}={(n}_{\max }-{n}_{\mathrm{ch}})/{\tau }_{{\rm{S}}}\) and \({\lambda }_{{\rm{R}}}={n}_{\mathrm{ch}}/{\tau }_{{\rm{R}}}\) during SET and RESET, respectively. Since \({n}_{\max } > {n}_{\mathrm{ch}}\) at any time, both generation rates are always positive as they must be. For the RESET transition, we will depart from an initial number of channels, ninit, which are the ones generated during the previous SET transition.
The events are generated with a random number u uniformly distributed between 0 and 1 along the simulation time. The simulation time is discretized in steps Δt which are small enough so that λ(t) can be assumed to be constant during Δt. It can be shown that under these conditions, the random time to a subsequent event at time t is \(\Delta {t}_{{\rm{u}}}=-\mathrm{ln}(u)/\lambda (t)\). During the simulation, if \(\Delta {t}_{{\rm{u}}} an event is generated at time t, otherwise, the event is rejected. Details on modelling are discussed in Supplementary Section 9.
Interlaboratory comparison
An interlaboratory comparison involving six participants was carried out for the electrical characterization of quantum conductance levels in memristive devices, with the aim of testing the intrinsic standard of electrical conductance (or resistance) and for evaluating laboratory-to-laboratory variability. For this purpose, samples assumed to be identical were distributed among participants and a common measurement protocol was defined. The participants were the following institutions: Istituto Nazionale di Ricerca Metrologica (Italian Institute of Metrology, NMI 1), Instituto Português da Qualidade (Portuguese Institute of Metrology, NMI 2), Turkiye Bilimsel ve Teknolojik Arastirma Kurumu (Turkish Institute of Metrology, NMI 3), Forschungszentrum Juelich GmbH (LAB 1), Fundación IMDEA Nanociencia (LAB 2) and Politecnico di Torino (LAB 3).
Measurement protocol
The equivalence of the measurements across the different laboratories was ensured by establishing and agreeing a measurement protocol that defines standardized measurement conditions to program, accept and stabilize the quantum conductance level, and defines the methodology to measure its conductance value under steady conditions (an example of the device programming methodology is reported in Extended Data Fig. 1). The generation of the quantum conductance states is achieved by running sequential SET/RESET cycles where an applied voltage to the two terminals of the device is swept between +1.5 V and −0.9 V. The positive part of the sweep (SET cycle) has a sweep rate of 96 mV s−1 (voltage steps of 50 mV). The negative sweep (RESET cycle) has a slower sweep rate of 2 mV s−1 (voltage steps of 1 mV). The current compliance was established as 500 µA and 10 mA for the positive and negative cycles, respectively. The voltage at the terminals of the device and the current that flows through it are continuously measured over SET/RESET cycles, and the corresponding conductance state is obtained for each applied voltage step. The formation of the quantum conductance steps during the RESET is continuously verified and a criterion to detect and accept G1 and G2 conductance states related to G0 and 2G0 quantum values, respectively, was established. If the last five consecutive measurements of the conductance state lay within either G0 ± 0.5G0 or 2G0 ± 0.5G0 (censoring interval), the sweep RESET cycle is interrupted, and a continuous read voltage of 10 mV is applied. The measurement of the step conductance value starts under this fixed applied control voltage and continues as long as it remains in the intervals [0.5G0; 1.5G0] or [1.5G0; 2.5G0]. The measurements were made at room temperature and under normal environmental conditions. The equipment used was a source meter (different equipment was used by the participants, as detailed in Supplementary Section 14) in autorange mode. The above-described methodology makes it possible to deal with the stochasticity of the conductive filament formation process establishing an initial limit to the variability around the nominal values of the desired quantum conductance steps (the validation of this programming methodology is discussed in Supplementary Section 15). Note that all measurements not strictly following the established comparison protocol were not considered for the interlaboratory comparison.
Evaluation of results and uncertainty budget
The evaluation of the average value and the variability of the programmed quantum steps was made from the observation of the measurements taken under repeatability and reproducibility conditions (described in appendix 2 of ref. 4). Here, repeatability conditions are understood as measurements of a specific device taken consecutively, while reproducibility is considered as the variability of the measurements taken from cycle-to-cycle operation of a specific device and programmed state as well as from device-to-device operations.
For each participant, j, the arithmetic mean and the experimental standard deviation were calculated for each series i of ni values:
$${\bar{G}}_{j,i}=\frac{1}{{n}_{i}}\mathop{\sum }\limits_{a=1}^{{n}_{i}}{G}_{i,a}$$
(7)
$${s}_{j,i}=\sqrt{\frac{1}{{n}_{i}-1}\mathop{\sum }\limits_{a=1}^{{n}_{i}}{\left({G}_{i,a}-{\bar{G}}_{j,i}\right)}^{2}}$$
(8)
The s.d. given by equation (8) is an estimate of the repeatability53,54 associated with series i of a programmed quantum conductance state measured by participant j. Only series with a minimum of 30 consecutive values and limited to a maximum of 100 values were considered as a fixed condition in this data evaluation (Supplementary Section 13). As each participant measured Nj series and there are series with different numbers of values, a polled standard deviation55 \({s}_{{\rm{p}},\,j}^{2}\) is calculated based on the following equation for its variance:
$${s}_{{\rm{p}},\,j}^{2}=\frac{{\sum }_{i=1}^{{N}_{j}}\left({n}_{i}-1\right)\times {s}_{j,i}^{2}}{{\sum }_{i=1}^{{N}_{j}}\left({n}_{i}-1\right)}$$
(9)
\({s}_{{\rm{p}},\,\,j}\) is therefore a weighted average of the \({N}_{j}\) s.d. where the number of degrees of freedom \(\left({n}_{i}-1\right)\) is the weight of each series.
For each participant, an average of the mean values obtained from the \({N}_{j}\) series and the experimental s.d. is calculated as:
$${\bar{\bar{G}}}_{j}=\frac{1}{{N}_{j}}\mathop{\sum }\limits_{i=1}^{{N}_{j}}{\bar{G}}_{j,i}$$
(10)
$${S}_{j}=\sqrt{\frac{1}{{N}_{i}-1}\mathop{\sum }\limits_{i=1}^{{N}_{j}}{\left({\bar{G}}_{j,i}-{\bar{\bar{G}}}_{j}\right)}^{2}}$$
(11)
The evaluation of the reproducibility of the programmed quantum conductance steps was based on the s.d.53,54 given by equation (11). Because the values obtained by each participant for each step are from different cycles and different devices, the reproducibility obtained is the result of cycle-to-cycle and device-to-device variability.
The measurement of quantum conductance states associated with each participant is expressed by the following measurement equation:
$${{G}_{j}={\overline{\bar{G}}}_{j}+{S}_{j}+{s}_{{\rm{p}},\,j}+e}_{j}$$
(12)
where \({\bar{\bar{G}}}_{j}\) is the mean value calculated by participant j, Sj is the related experimental s.d. according to equations (10) and (11), sp, j is the repeatability of the measurements according to equation (9), and ej is the error related to the accuracy of the measurement equipment used. It is assumed that these input variables are statistically random variables where Sj, sp, j and ej have an expectation value equal to zero and a s.d. estimated based on the experimental values presented before (Sj and sp, j) and in the manufacturing specifications of the equipment used (for ej). Note that random effects, including cycle-to-cycle variability but also variations related to small variations in the room temperature, humidity levels or even small fluctuations from the measurement set-up, are included in the estimation of the uncertainty component of the quantities Sj and sp, j, even if each specific contribution has not been disentangled.
The measuring uncertainty of Gj can be estimated by applying the law of propagation of uncertainties55 to equation (12):
$${u}^{2}\left({G}_{j}\right)={u}^{2}\left({S}_{j}\right)+{u}^{2}\left({s}_{{\rm{p}},\,j}\right)+{u}^{2}\left({e}_{j}\right)$$
(13)
where \({u}^{2}\left(x\right)\) is the variance (square of standard uncertainty) associated with the variable x and u2(Gj) is the square of the combined uncertainty of Gj.
The standard uncertainties of Sj and sp, j are estimated by the corresponding s.d. of the mean:
$$u\left({S}_{j}\right)=\frac{1}{\sqrt{{N}_{j}}}{S}_{j}$$
(14)
$$u\left({s}_{{\rm{p}},\,j}\right)=\frac{1}{\sqrt{{\sum }_{i=1}^{{N}_{j}}\left({n}_{i}\right)/{N}_{j}}}{s}_{{\rm{p}},\,j}$$
(15)
The relative standard uncertainty of ej is calculated from the combined relative uncertainty of the measurement of the voltage, \({u}_{{\rm{r}}}^{2}\left(U\right)\), and current, \({u}_{{\rm{r}}}^{2}\left(I\right)\):
$${u}_{{\rm{r}}}\left(e\right)=\sqrt{{u}_{{\rm{r}}}^{2}\left(U\,\right)+{u}_{{\rm{r}}}^{2}\left(I\,\right)}$$
(16)
The relative uncertainties of the measured voltage U and current I are estimated assuming a rectangular probability distribution for the voltage and the current measuring error with the plus/minus limits given by the manufacturing specifications of the equipment, usually identified as ‘accuracy’ (Supplementary Section 14):
$${u}_{{\rm{r}}}\left(U\,\right)=\frac{1}{\sqrt{3}}\frac{{U}_{\mathrm{accuracy}}}{U}$$
(17)
$${u}_{{\rm{r}}}\left(I\,\right)=\frac{1}{\sqrt{3}}\frac{{I}_{\mathrm{accuracy}}}{I}$$
(18)
Following the international recommendation to express the final measuring uncertainty with a coverage probability of approximately 95%56,57, the expanded uncertainty U(Gj) is calculated following the equation:
$$U\left({G}_{j}\right)=k\times u\left({G}_{j}\right)$$
(19)
where k is the coverage factor calculated according to annex G of ref. 55.
Evaluation of consensus value
The evaluation of the results achieved by the participants was done by comparing individual results with a consensus value36,37. The consensus value is established based on all results from the participants37, using a weighted average of their values35:
$${G}_{\mathrm{cons}}=\left(\sum _{\,j=1}^{6}{w}_{j}\times {G}_{j}\right)\left/\left(\sum _{\,j=1}^{6}{w}_{j}\right)\right.$$
(20)
where the weighting factors are given by:
$${w}_{j}=1/{u}^{2}\left({G}_{j}\right)$$
(21)
The combined uncertainty of the consensus value is estimated based on the participant uncertainties as follows:
$$u\left({G}_{\mathrm{cons}}\right)=\sqrt{1\left/\sum _{j=1}^{6}\right.{w}_{j}}$$
(22)
And the related expanded uncertainty is given assuming a coverage factor k = 2 (ref. 35):
$$U\left({G}_{\mathrm{cons}}\right)=2\times u\left({G}_{\mathrm{cons}}\right)$$
(23)
To identify an overall consistency of the results produced by this approach, a chi-square test was applied to the input values35:
$${\chi }_{\mathrm{obs}}^{2}=\mathop{\sum }\limits_{j=1}^{n}\left[{\left({G}_{j}-{G}_{\mathrm{cons}}\right)}^{2}/{u}^{2}\left({G}_{j}\right)\right]$$
(24)
The result of the test is considered to fail if \(\Pr \left\{{\chi }^{2}\left(\nu \right) > {\chi }_{\mathrm{obs}}^{2}\right\} where Pr is the ‘probability of’, \({{\chi }}^{2}\left(\nu \right)\) is the expected theoretical value of a chi-squared distribution for \(\nu\), and \(\nu\) is the degrees of freedom, which is the number of input values n minus 1 (in this case, 5). If the consistency check does not fail, then Gcons can be accepted as the consensus value and U(Gcons) can be accepted as its expanded uncertainty. Values obtained for the interlaboratory comparison were \({\chi }_{\mathrm{obs}}^{2}=6.3\) and \({{\chi }}^{2}\left(5\right)=11.1\). As \({\chi }_{\mathrm{obs}}^{2}\le {\chi }^{2}\left(5;0.05\right)\), the consistency of the participant’s values and the calculated consensus value was demonstrated, thus the obtained Gcons is the consensus value and U(Gcons) is its expanded uncertainty.
To qualify the result of each participant related to the consensus value, the normalized error35,37, En, j, was calculated by:
$${E}_{{\rm{n}},\,j}=\left({G}_{j}-{G}_{\mathrm{cons}}\right)/\sqrt{{U}^{2}\left({G}_{j}\right){-U}^{2}\left({G}_{\mathrm{cons}}\right)}$$
(25)
The value of En, j has the following meaning: if |En, j | ≤ 1.0, the result is consistent (passed); if |En, j | > 1.0, the result is inconsistent (failed). For all participants, results were observed to be consistent with the established consensus value. Based on statistical analysis, higher values of |En, j | (even if always ≤1.0) cannot be ascribed to eventual systematic errors affecting the measurement that are not being adequately corrected or considered in the evaluation of measurement uncertainty.

