Filamentary conduction models for RRAMs rely on the idea that the current magnitude is somehow related to the presence of a potential barrier or alternatively to a gap along the CF. For example, the QPC model [40,54–56] considers that the current flow through a single nanosized filamentary structure is controlled by the narrowest section along the filament. In terms of the physics of mesoscopic devices, this section has associated a transmission probability like any other barrier or scatterer. Following this idea, Miranda [34] and Patterson et al. [35] proposed a TE based on the similarity of the QPC model expression for the I-V curve with two identical opposite-biased diodes in series with a resistor, as shown in the inset of Figure 2A. Such a TE allows a progressive transition between an exponential (HRS) and a linear (LRS) curve by simply changing a single parameter of the model 0<λ < 1 called the memory state of the device [34].

According to the QPC model, the equation for the I-V characteristic of a memdiode M is given by Eq. 1: I = I 0 ( λ ) [ e β α ( λ ) ( V − I R S ( λ ) ) − e − ( 1 − β ) α ( λ ) ( V − I R S ( λ ) ) ] (1) where I 0 ( λ ) = I m i n ( 1 − λ ) + λ I m a x is the current amplitude factor, α a parameter related to the particular features of the conduction mechanism, and R _S a series resistance. Eq. 1 is the result of considering an inverted parabolic potential barrier with height Φ as the tunneling scatterer. Both α and R _S can have the same dependence on λ as I 0 . I m i n and I m a x are the minimum and maximum values of the current amplitude, respectively, and 0<β < 1 expresses the asymmetry in the potential at the two ends of the CF (in this work β = 1/2 is assumed). For a monomode ballistic conductor, I 0 = G 0 / α ⋅ e x p ( − α Φ / e ) where G 0 = 2 e 2 / h is the quantum conductance unit (e the electron charge and h the Planck’s constant). As I 0 increases in Eq. 1 because of the reduction of Φ (widening of the cross-section area of the CF), the I-V curve changes its shape from exponential (HRS) to linear (LRS) as experimentally observed for this kind of devices. This occurs because of the voltage drop across R _S and the reduction of the effective applied voltage to the constriction. Notice that Eq. 1 is equivalent to two opposite-biased ideal diodes (adding a + 1 and a − 1 in the right-hand side of the equation) for which the model was named memdiode. Throughout this work, “memristor” and “memdiode” are used interchangeably to refer to the same device.

From the modelling viewpoint, the ME is normally a voltage or a current driven differential equation [47] that usually relies on some non-linear dependencies as well as on the so-called window function. Although the introduction of a window function in the ME acting as a boundary condition for λ [32] represented a breakthrough in the modelling of the SET/RESET transitions, the approach was demonstrated not to be exempted from serious mathematical drawbacks [33,57]. The equation discussed here complies with a number of requirements such as dimensional homogeneity, dynamic balance and levelling-off behavior for large opposite biases, with end values 0 and 1 for HRS and LRS, respectively. Physically, the ME expresses the ion/vacancy movement caused by the external applied field. Among the possible candidates, the simplest first order differential equation is given by Eq. 2, where τ S , R are characteristic times linked to the SET and RESET transients (for V > 0 and V < 0 , respectively). d λ d t = 1 − λ τ S ( V ) − λ τ R ( V ) (2) τ S , R ( V ) = τ 0 S , 0 R e − V V 0 S , 0 R (3)

In this work, we assume an exponential dependence of τ S , R on the applied voltage, as expressed by Eq. 3 and Figure 2B. V 0 S , 0 R and τ 0 S , 0 R are fitting constants ( V 0 S > 0 for SET and V 0 R < 0 for RESET). Although Eq. 2 has already been used in connection with physical parameters of the CF [58–61], it was first employed as a behavioral memory equation in [62]. A central feature of Eq. 2 is that it can be represented as an RC circuit with voltage-dependent resistors as illustrated in Figure 2C. In this representation, λ corresponds to the voltage drop across the capacitor C H , whose initial value is λ 0 . This representation substitutes the voltage-controlled memory subcircuit considered in many previous memristor models [33].

For a periodic ramped input signal, the combination of Eqs 1-3 results in I-V loops such as those illustrated in Figure 2A, which starts with HRS (λ = 0) and evolves as indicated by the blue arrows printed on top. The SET and RESET processes triggered at the transition voltages V _TS and V _TR point out the switching from HRS to LRS and vice versa, respectively. The ratio between the LRS and HRS currents measured at a fixed low bias is referred to as the resistive window of the device and when properly normalized to unity represents the memory state window. Note that as the signal frequency increases, V _TS and V _TR shift toward higher absolute values, as indicated by Eq. 4 [63]: V T S = V 0 S ⁡ ln ( R R ) + V 0 S ⁡ ln ( τ 0 S V 0 S ) (4) where RR is the input signal’s ramp rate [V(t) = RR⋅t]. A similar expression holds for V _TR . These SET and RESET voltage shifts lead in turn to the mutual approximation of the LRS and HRS curves and thus to the reduction of the hysteresis lobe area. This is in agreement with previously reported experimental measurements of the I-V loops as function of frequency, as shown in Figure 3A. The evolution of the SET and RESET voltages is plotted separately in Figure 3B and fit with Eq. 4. The separation between the lines is a consequence of the series resistance effect [63]. In physical terms, the collapse of the memory window is attributed to the incapability of ions/vacancies to follow the input signal [64,65]. From the electrical viewpoint this effect can be directly related to the shift of the switching thresholds [66–68].

For a given RR, the HRS (exponential) to LRS (linear) transition as λ is swept from 0 to 1 is detailed in the inset of Figure 2A (solid blue lines), altogether with some intermediate states (dashed blue lines). For comparison purposes, a linear model [41,69] (faded-thick grey lines) is also plotted. sinh()-based models are omitted here as they require the simultaneous fit of multiple parameters to mimic the smooth linear-exponential to linear transition or even separate expressions for the HRS and LRS regimes [70,71]. Although the two models coincide at low voltages and exhibit a clear linear dependence, significant discrepancies arise as the voltage increases. As can be seen, the linear model cannot capture the departure of the HRS curves at intermediate voltages. As such, if used to fit I _LRS and I _HRS (the currents in LRS and HRS, respectively) at a nominal read voltage (V _read ), the linear model leads to an overestimation of the device current when lower voltages are applied. On the contrary, the DMM can accurately describe both the HRS and LRS curves by solely changing a single parameter (λ) in the TE.

Concerning the switching dynamics for discontinuous signals, Figures 2D–G illustrate the effects of a sequence of equal amplitude voltage pulses (V _write = 0.6–1.3 V) and period (indicated as a combination of frequency ranging from 100 to 100 kHz and Duty Cycle—0.2 to 0.8 –) on the time required to reach a given current level. The parameter values for the ME are V 0 S , 0 R = 6.77 × 10 − 2 V and τ 0 S , 0 R = 8.48 × 10 3 sec. Of course, these values depend on the switching material under consideration. Such input stimuli replicate those employed in OLT or CLT approaches as proposed in previous works [23]. For each pulse, the current increases as a function of the voltage and time due to the so-called potentiation effect in memristors [72]. Knowing the memory state value increase as a function of the number of pulses N with amplitude V applied during a time Δt, λ N = ( λ 0 − 1 ) e − ( NΔ t / τ S ) + 1 (5) the time required to reach a target current level can be estimated from Eq. 6 as T w r i t e = − τ 0 S e − V V 0 S D C ln ( I L R S − I T I L R S − I H R S ) (6) where I _LRS e I _HRS are the maximum and minimum currents measured at the pulse’s amplitude V, respectively. DC is the duty cycle. It should be noted that Eq. 6 is only valid in the linear region of the HRS regime, i.e. for a voltage V lower than V _TS . Also, it is assumed that the hysteresis lobe does not collapse in the region of interest. Otherwise the formula is not valid unless the variation of I _HRS and I _LRS are considered. In this region, the write time increases exponentially as the pulses’ amplitude decreases, but remains insensitive to their frequency as shown in Figures 2D,E. In addition, notice that as V _TS shifts to higher values due to the increase of the input signal frequency (RR), so it does the range in which Eq. 6 holds valid: it goes from 0-∼0.8 to 0-∼1.2 V when the frequency increases from 100 to 100 kHz (see Figure 2D). Similarly, when sweeping the frequency of the pulses, Eq. 6 holds valid for frequencies above 100 Hz at V _write = 0.8 V (the amplitude of the pulses is below V _TS for any frequency in the range considered) and for frequencies above 100 kHz at V _write = 1.2 V. Rather than the signal frequency, the write time is slightly sensitive to the pulse’s width, which is expressed in Eq. 6. as a dependence on the signals’ DC as shown in Figures 2F,G. As it can be seen, when increasing DC, the write time down-shifts as each signal period causes a larger conductance variation.

The model discussed in the previous Sub-Sections was put under test by fitting experimental data extracted from different published works. In particular, Figures 4A–C show the results obtained for different RRAM structures including HfO₂ [73], Ta₂O₅/TaO_X [74] and TaO_X [75] structures, respectively. In all the cases the I-Vs were measured at room temperature and under voltage sweeps. The experimental data were fitted using the SPICE model depicted in Table 1 based on Eqs 1–3 and applying driving signals as described in the corresponding references. The fitting parameters are listed in each of the panels of Figure 4 as reference, as well as the details of the stack structure. It should be mentioned that the proposed DMM does not only provide a simple SPICE-compatible implementation for the resistive memory devices but also a versatile one, as it can accurately fit the I-V loops experimentally measured in different RRAM devices.

Moreover, we have tested the capability of the DMM to replicate the Long-Term Potentiation (LTP) and Depression (LTD) of memristive synapses used in neuromorphic hardware. This kind of evolutionary behavior is required to achieve gradual conductance changes upon pulse applications. Figure 4D shows the LTP and LTD measurements reported in [75] for TiN/(25 nm)TaO_X/Pt -based devices (the experimentally measured I-V loop is plotted in Figure 4C altogether with the corresponding fit with the DMM) by the application of 300 identical pulses of 1 V (LTP) followed by 300 identical pulses of -1.1 V (LTD). All pulses have the same width (100 ns, t _ON ) and they are applied every 20 msec During the time in between pulses, a low voltage (0.1 V) pulse is applied to read the conductance (memory) state of the device. These measurement conditions were replicated in the SPICE simulations shown in Figure 4C. The simulated LTP and LTD trends are superimposed to the measurements shown in Figure 4D, showing a good fit of the experimental trends. There is also a voltage acceleration of the LTP/LTD trends, as reported for instance for the TiN/(3-nm)HfO₂/Pt -based devices measured in [76], and which can be reproduced by the DMM as shown in Figure 4E. Finally, the versatility and capability of the DMM to reproduce the evolutionary behavior of memristive devices is shown (see Figure 4F) by also fitting the gradual current increase during LTP of Ag/ZnO/Pt based nanowire memristors [77], by considering the same stimuli (pulses with an amplitude of 2.5 V and a duration of roughly 2 msec.). Thereby, the DMM is suitable to model the response of RRAM devices with a large number of incrementally accessible conductance states.

First, a software-based SLP of size n ² × m (therefore with n ²·m synapses) is created and trained using a given dataset formed by n × n px. images (as those represented in Figure 5B) distributed into m classes. For the sake of simplicity, ex-situ supervised learning is assumed here with the Scaled Conjugate Gradient [78] (SCG) training algorithm. This algorithm provides a good trade-off between accuracy and learning time for the different datasets considered. The possible impact of the learning method was discussed elsewhere [14] where it was demonstrated that no significant statistical differences in terms of the inference accuracy occur. The training procedure generates a n ² × m weight matrix W M ∈ ℝ , which can be further decomposed as W M = W M + − W M − based on Eq. 7 and Eq. 8, as proposed in the literature [79]: w M i , j + = { w M i , j , w M i , j > 0 0 , w M i , j ≤ 0 (7) w M i , j − = { 0 , w M i , j > 0 − w M i , j , w M i , j ≤ 0 (8) with both W M + and W M − comprising only positive elements. This allows rendering both positive and negative synaptic weights in W M as well as doubling the dynamic range and reducing the noise and variability susceptibility [80].

In the next step, the conductance matrices G M + and G M − (each sized n ² × m, thus totaling 2·n ²·m synapses in the SLP) to be mapped into the CPAs are calculated by the linear transformation [22,81]: G M + ( − ) = ( G m a x − G m i n ) W M N o r m + ( − ) + G m i n (9) where [ G m i n ,   G m a x ] is a selected conductance range for a linear computation in matrix-vector calculations. For simplicity, we consider G m a x = G L R S = 1 / R O N and G m i n = G H R S = 1 / R O F F , with R _ON and R _OFF defined by the I-V characteristic of the memdiode model play evaluated at a given V _read . Then W M N o r m + ( − ) is the positive (negative) weight matrix normalized within the range [0,1] by one of the following normalization methods NM ₁ -NM ₂ , indicated by Eq. 10 and Eq. 11: M N 1 :   W M N o r m + ( − ) = W M + ( − ) m a x { a b s { W M } } (10) M N 2 :   w M i , j = { 1 , w M i , j > μ W M +   n σ W M w M i , j μ W M −   n σ W M , μ W M −   n σ W M < w M i , j < 0 w M i , j μ W M +   n σ W M , 0 < w M i , j < μ W M +   n σ W M − 1 , w M i , j < μ W M −   n σ W M (11) m a x { W M } and m i n { W M } are the maximum and minimum synaptic weight values in W M and μ W M , σ W M are the mean and standard deviation of the synaptic weights in W M , respectively. In this way, the W M + and W M − matrices are converted to conductance values within the range [ G H R S ,   G L R S ] .

The subsequent sub-routines write down the circuit netlist for the dual-n ² × m memdiode CPA-based SLP, adding the parasitic line resistance (R _L ), connection scheme, and control logic necessary for both the CLT and inference phase. Two approaches were followed in order to improve the voltage effectively delivered to the RRAM cells as shown in the simplified equivalent circuit depicted in Figure 5C: i) A Dual Side Connection (DSC) scheme and ii) the partitioning of the G M + and G M − matrices into NP (number of partitions) smaller arrays [14,69,80]. Despite the increased peripheral circuitry complexity, DSC improves the voltage delivery to each synapse [13] by connecting both wordline terminals to the input stimuli. Similarly, the small size of each partition helps reducing the parasitic voltage drop along the line interconnections. Exploding the integrability of the CPA with CMOS circuitry, the vertical interconnects linking the outputs of each partition may be placed under the CPAs, as well as the analogue sensing electronics, thus minimizing the area overhead of the partitioned architecture [80].

To account for both the inference and write phases, the present implementation follows a reconfigurable approach, in which the partitioned CPA-based SLP is alternatively connected to the input/output signals (inference phase, Figure 5C) or to the writing stimuli (write phase, Figure 5D). The analogue electronics required for the inference phase has been described elsewhere in the literature [14,15]. However, the circuitry needed for the write phase requires a more complex circuital implementation, as the input stimuli are not passed simultaneously to all the CPA inputs, but sequentially. Thereby, as shown in Figure 5D, this circuital arrangement involves a CPA address block, the Row/Column address decoders, the Row/Column selectors, and the Write Acknowledge block. All of them have been designed assuming a 130 nm CMOS process from GlobalFoundries.

The CPA address block is basically an asynchronous counter that produces n ² /NP m output pulses (the number of memdiodes in each partition), which propagate to both the Row and Column decoders. Both decoders consist of another ascending asynchronous counter that counts up to n ² /NP (Row decoder) and m (column decoder) pulses, and outputs a binary code indicating the row and column being addressed, respectively. Then, these binary codes are simultaneously propagated to all the Row/Column selectors. The latter are an array of analogue switches that connect the input node of each CPA Row to V _write or V _read (for addressing that particular Row during the write procedure), V _DD /2 (if another row is being addressed) or to V _input (when the ANN is operating in the inference phase). The column selector is a similar array that connects the columns output nodes to a sensing amplifier (SA, a transimpedance amplifier coupled to a comparator) if that particular column is being addressed, or V _DD /2 (if other column is being addressed). Finally, the target current level to be detected by the SA to determine the completion of the write process of a given element g _i,j in the G M + and G M − matrices is calculated based on the nominal V _read and conductance g _i,j . Once the circuit netlist is written, it is passed to a SPICE simulator which evaluates the voltage and current distributions in the CPA circuit while it is being programmed, and then, as the input images are processed and classified. Finally, the resulting waveforms are passed back to the MATLAB routine for evaluation and metrics extraction.

During the write operation each memdiode in the CPAs (M _i,j ) is individually addressed and supplied with a train of alternating read and write pulses of amplitude V _read and V _write , respectively. This causes a gradual increment of the memdiode conductance until a target value is reached. Such addressing procedure is performed following the V _DD /2 approach as it minimizes the line disturbance [41]. Within this writing method, the non-addressed rows are set to a constant voltage value V _DD /2. Similarly, the output node of the column of the addressed memdiode is grounded through the SA, which measures the current flowing out from this column (the other columns are at V _DD /2). The measured current is proportional to the applied voltage pulses and the memdiode conductance (g _i,j ) plus the parasitic wire resistance corresponding to the addressed device (M _i,j ). This allows to estimate the conductance of the addressed memdiode. This process is represented by the simplified equivalent circuit shown in the inset of Figure 6.

The write procedure for the addressed memdiode M _i,j begins by sensing the output current during the read pulse of magnitude V _read . In case this current is lower than a target value, a write pulse of voltage V _write is applied (V _write > V _read ), causing an increment of the M _i,j conductance (g _i,j ). Then a new read pulse is applied and the current sensed again. This process continues iteratively until the sensed current during the read pulse meets the target value. Once reached, the SA outputs a pulse that indicates the completion of the writing procedure for the addressed memdiode (M _i,j ), stopping the train of read/write pulses and preparing the following devices for the programming step.

It is worth noting that the partitioned architecture allows the simultaneous programming of the M _i,j memdiode in all partitions using a simpler control circuit. Let us assume that the devices to be programmed are the M _i,j memdiodes of a 2 × (n ² × m) RRAM-based SLP implemented with NP partitions (thus totalling 2·NP CPAs), such as the one presented in Figure 5D. In this case, the i ^th output of the row decoder (n ² /NP outputs) will be the only active output, as well as the j ^th output (m outputs) of the column decoder. Then these output vectors are passed to every Row/Column selector, which therefore simultaneously select the M _i,j memdiode in every CPA. This causes all the i ^th rows to be connected to a train of alternating read and write pulses and all the j ^th columns to be connected to the partition’s SA (each CPA partition has its own SA). All other rows and columns are connected to V _DD /2. The current flowing through each of the M _i,j memdiodes (and therefore out from the j ^th columns) is sensed by its associated SA until the target conductance value for that M _i,j memdiode is achieved. Then the associated SA propagates an acknowledge pulse (ACK) to the Write Acknowledge block. This block waits for the ACK pulses coming from the SAs of every partition. Once all ACK pulses are received, the i ^th ,j^th position of all CPAs is considered to be successfully written, and by the time the Write Acknowledge block receives the following system clock pulse, it instructs the CPA address block to address the M _i,j+1 memdiode and the write sequence starts again. This process continues until the CPA address block has addressed all the memdiode positions in the CPA partitions (n ² /NP × m positions).

The MNIST (Modified National Institute of Standards and Technology) dataset of handwritten digits was selected for the training and inference phases considered in this work. This dataset comprises a series of k input feature vectors [x(k) for sample k] and a target output vector [t(k), with 10 dimensions, each corresponding to one digit]. For the classification problem, t _c (k) =1 if sample k belongs to class c, otherwise t _c (k) = 0. The input feature vectors (n ² × 1) are the unrolled grayscale pixel values of thetwo-dimensionalimages (n × npx.). Pixel’s brightness is codified in 256 gray levels between 0 (fully OFF, black) and 1 (fully ON, white). The MNIST dataset contains 60,000 training images and 10,000 testing images, both in grayscale and with a 28 × 28 px. resolution [82]. A few examples of these images can be seen in the left panel of Figure 5B. The MNIST images were further down sampled using bicubic interpolation to 8 × 8 px. to allow smaller SLPs. This modification reduces the inference accuracy degradation due to the line resistance effect and speeds up the simulations (also reduces the image readability as shown in the right side of Figure 5B). Finally, the input stimuli are obtained by scaling the input feature vector by a voltage V _read . V _read is chosen such as to prevent altering the memdiode state during the inference phase simulation. In this way, during the inference process each of the test images is presented to the CPA as a vector of analog voltages in the range [0, V _read ].

All three metrics considered in this study show similar trends regardless of the normalization method, R _ON and R _OFF resistances and line resistance (R _L ). When addressing the CPA write time presented in Figure 7A, it can be seen that as the voltage of the write pulses decreases, the CPA write time slightly increases until reaching roughly 1.1 V where it increases abruptly. The reason behind this behavior is a consequence of the I-V characteristics of the memristors themselves, in particular the SET transitions. Above the 1.1 V threshold voltage, the reduction of the write pulses amplitude causes the conductance increase produced by the application of each pulse to decrease too, thereby requiring more pulses to program each memdiode in the CPA. However, when applying a voltage lower than the 1.1 V threshold, the voltage delivered to a number of cells in the CPA is not enough to alter its conductance due to a reduced read (or write) margin. RM (WM) is the fraction of the applied input voltage V _in (which may be V _read or V _write ) effectively delivered to the memory cells (V _cell ), i.e. V _cell /V _in .

This interpretation is in close agreement with the results reported in Figure 7B, where the SWV metric is represented in terms of the write voltage. Note that as the conductance change produced by each write pulse reduces with the pulse’s amplitude, the write procedure has higher control over the programmed conductance, and thereby the write error reduces to a minimum. Nonetheless, when crossing the 1.1 V threshold, a subtle rise in the SWV can be seen, coincident with the abrupt increase of the write time. This could be explained by considering the sneak-path effect. In this context, the programming of certain positions in the CPA is not terminated by the current being driven by the device located in that position, but by parasitic currents loops caused by already set devices located close to the input/output ports. A direct consequence of the write error increase is the degradation of the inference accuracy as shown in Figure 7C. Note that both below the 1.1 V threshold and above 1.6 V there is a reduction of the inference accuracy that almost coincides with the increasing writing errors during the memristor programming phase.

The matrix elements in W M are in the range [ m i n { W M } , m a x { W M } ] . To be mapped to a conductance magnitude in the range [ G H R S , G L R S ] , they must be first normalized to the range [ − 1,1 ] . Very often [22,81], such normalization is achieved dividing W M by the absolute value of the maximum element in W M (Normalization Method 1, NM-1, see Eq. 10). As expected, the normalized W M matrices W M N 1 preserve the exact same distribution and the m a x { W M } / m i n { W M } ratio. Interestingly, for the case of the MNIST images resized to 8 × 8 px. used in this work for benchmarking, ∼95% of the elements from W M N 1 falls within the range [−0.5,0.5]. Thereby, NM-1 does not exploit the entire dynamic range of the memdiodes as most of the devices will be set to a conductance value in the range [ G H R S , ( ( G L R S + G H R S ) / 2 ) ] . To make a better use of the entire memristor conductance range, an alternative approach (NM-2, see Eq. 11) is considered. In this context, an element w i , j ∈ W M has a probability P _i of being within the range μ W M ± i σ W M , where μ W M and σ W M are the mean and standard deviation of the values of W M . For i ranging from 1 to 4, ∼68.3%, ∼95.5%, ∼99.7 and ∼99.9% of the synaptic weights will be within such range, respectively [84]. Thus, values exceeding such limits are set equal to μ W M ± i σ W M and then W M is normalized to obtain W M N 2 . This normalization method allows not only a better exploitation of the memristor dynamic range but also minimizes the impact of stuck-at-ON faults, as explained in Ref. [85].

The impact of both NM-1 and NM-2 were compared in terms of the proposed metrics in Figure 7, with the CPA write time and SWV showing a factor 2 of difference between the two normalization methods. Given that, when considering NM-1, most of the conductance values are mapped close to G _HRS , a lower number of write pulses are required to meet the target conductance. This results in shorter writing times, as seen in Figure 7A. The lower mean conductance of the synaptic weights resulting from applying NM-1 also results in a lower SWV (See Figure 7B). As the conductance of the synapses increases, so does the current flowing from wordlines to bitlines. This not only causes a larger IR-drop in the interconnection lines, but also exacerbates the sneak-path effect too, resulting in a similar scenario to that depicted in section Simulation Results. A: the current measured by the sensing amplifier that terminates the writing cycle in a given {i,j} position of the CPA is in fact the result of multiple conductive paths. As a consequence, the writing cycle in the {i,j} position is terminated before reaching the target conductance, giving rise to a significant increase of the SWV metric. Finally, no significant changes occur in the inference accuracy of the programmed CPA (see Figure 7C), except for a higher sensitivity to the line resistance. This point will be thoroughly discussed in the next section.

Assuming conventional metal wires, the resistance of the interconnect wires (R _L ) between two adjacent cells in a 4F ² CPA structure can be estimated to be ∼4.5, 3, and 1.5 Ω for the 16, 22 and 32 nm technology nodes, respectively [86]. Nevertheless, for the novel technology nodes (10 nm and below) both surface and grain boundary scattering cause a size-dependent resistivity of the Cu wires [87–89] as the mean free path of electrons becomes comparable to the wire dimensions. These two effects are well-known and can be quantified using the Fuchs-Sondheimer (FS) [90] and the Mayadas-Shatzkes (MS) [91] models. They reveal that for highly scaled nodes (∼5 nm), R _L can be as large as ∼100 kΩ [69]. In this work, we have considered R _L values varying from 1 to 60 Ω.

Interestingly, all the three metrics (CPA Write Time, Inference Accuracy and SWV) show a much higher sensitivity to R _L when considering the NM-2. This can be attributed to the CPA’s RM. RM is jointly determined by R _L and the memdiode resistance (R _memd , which varies in between R _OFF and R _ON ) [14,41,69]. From a very basic analysis, each memdiode is part of a conductive path between the CPA’s input wordline i and the output bitline j. For a N×M CPA, the average parasitic resistance associated with such path is R _L [(N + M)/2 + 1] [14,37]. Within this simplified scenario, the V _cell /V _in ratio can be obtained from the voltage divider between R _memd and R _L [(N + M)/2 + 1]. This can explain the greater robustness of the NM-1 against R _L variations: the larger average value of R _memd minimizes the IR drop and as a consequence, the impact of R _L . On the other hand, for NM-2, the IR drop across R _L significantly limits the minimum write voltage, as for increasing values of R _L the voltage effectively delivered to the CPA cells reduces causing an exponential increase in the CPA write time. This effect also worsens the impact of the sneak-path problem, resulting in larger SWV values and lower inference accuracy (see Figures 7B,C, respectively).

As mentioned in Section Dynamic Memdiode Model, the I-V characteristic of the memdiode is sensible to the frequency of the input signal [63]. Thereby studying the impact of the write pulses frequency on the CPA write time, SWV and inference accuracy is necessary to complement the trends reported in Figure 7. In this regard, the aforementioned metrics were evaluated for different combinations of write pulses voltages and frequencies and for four different I-V loops, in which the abruptness of the SET transition is progressively increased. For the sake of simplicity, R _L is kept constant and equal to 10 Ω, as well as the amplitude of the verify pulses (0.3 V) and the read frequency and amplitude (100 kHz and 0.3 V).

In general, the simulation results presented in Figure 8 indicate that the best performance metrics are obtained for low voltage write pulses of high frequency (upper right corners of the SWV and accuracy metrics). This is a direct consequence of the smaller conductance change induced by each pulse when reducing their amplitude and width. On the contrary, in all the cases, the lower left corner is associated with the highest SWV and consequently with the lowest accuracy, given the lesser controllability of the conductance increase (pulses of high amplitude are applied during longer periods). Nevertheless, such a better controllability comes at the cost of a much higher writing time given that a larger number of pulses is required. This was already seen in the results reported in Figure 7.

Beyond the overall trends, the abruptness of the SET transition (see Figure 8A) also has a remarkable impact on the evaluation metrics. Note that as the SET transition becomes more abrupt, the SWV increases and consequently, the inference accuracy reduces. The changes can be regarded as an up-shift in the heat-maps presented in Figure 8, as no significant variations can be observed as a function of the frequency. This change in the heat-maps can be interpreted once again as a consequence of different conductance modifications caused by the pulses amplitude for the four different SET transitions considered: for a pulse with a given amplitude and width, a more progressive SET transition results in a smaller conductance increment, whereas the conductance modification produced for the same pulse in a memdiode with a steeper SET transition is larger.

Interestingly, apart from modifying the write performance of the conductance values, the frequency and voltage dependences of the memdiode also play a role during the inference phase, mainly because of the resistance shift occurring in the CPA memristors. Depending on the maximum voltage V _read used to translate the input pattern into a vector of analogue voltages, and the frequency at which such vectors are presented to the CPA, a resistance shift for the stored synaptic weights may occur (see Section Switching Dynamics). Notice that this is ultimately a consequence of the ionic movement inside the dielectric. This is presented in Figure 9A in terms of the SWV computed for the control parameter λ after having presented a series of images (x axis) and for different combinations of V _read values and input frequencies. The overall trend is that the SWV metric increases as more images are presented to the CPA and this can be minimized both by increasing the input frequency or reducing the read voltage. Such an increase in SWV causes inevitably a reduction of the inference accuracy as illustrated in the heat-maps in Figures 9B–E. Note that four cases with different SET abruptness were considered, and similarly to the results presented in Figure 8, the more abrupt the SET transition, the higher the sensitivity to the read voltage and the input frequency. As a corollary, it can be seen that V _read could be largely increased to improve, for instance, the Signal-to-Noise ratio as long as the input frequency is increased accordingly.

In this regard, it is crucial to consider the stability of the programmed resistive state of the cells. We have thereby studied in detail the distribution of synaptic weights and conductance matrices over time and reported the results in Figures 9F–T. This has been carried out by stopping the simulation and reading the value of the state variable λ for each memristor in the CPA (this value corresponds to the voltage of the node H in the model script). When the simulation is stopped, all the nodal voltages are stored and the λ-SWV computed with respect to the original mapped values. Then, the stored nodal voltages are assigned as initial conditions to the corresponding nodes and the SPICE simulation is resumed starting from the previous processed image. Note that we have considered accelerated stress conditions, V _read and Frequency values of 800 mV and 1 kHz. Under these circumstances the trends can be obtained in reduced times (the simulated time is 10 s for 10,000 images). Normal operation voltage (approx. 200 mV) and frequency (above 1 MHz) should ensure that the weights remain stable over a long period of time, otherwise, re-tuning of the synaptic weight shall be performed periodically. From Figures 9F–T, the statistical distribution of the programmed weights shift towards G _LRS as the number of images increases (i.e., over time). When considering the changes in the conductance matrices associated with these shifts in the statistical distribution of the synaptic weights, it can be seen that a certain pattern emerges, with the upper and lower region of the CPA arrangement showing no changes. Moreover, equidistant stripes of devices with no variations in middle of highly potentiated devices are observed. Such changes can be explained considering the particular features of the dataset analyzed in this work (MNIST). In the MNIST data-base, pixels located at the edges of the images are normally OFF. When unrolling the n×n images into n ² × 1 column vectors, the pixels from the left border of the image (normally unactive) occupies the first n positions of the n ² × 1 column vector, and those on the right border, the last n positions. Then the firsts and lasts n rows of the CPA arrangement are never stressed, so that they are always biased with approx. 0 V. Then, moving forward with the same reasoning, the upper and lower pixels of the middle columns of the n × n images are also normally unactive. When feeding the CPA arrangement with a n ² × 1 column vector, this results in equally spaced unbiased rows, which thereby suffers no variations of the weights stored in the corresponding memristors.

Finally, it is worth recalling that even for the case of high stability of the programmed memristor devices, a possible device-to-device (D2D) variability and the consequent programming error need to be considered. The intrinsic stochastic nature of the switching mechanism can induce a high degree of variability [92]. Such variability in the memristive structures also increases the possibility of a deviation of the weighted sum from the target conductance value [93]. The variability of the conductance states G _LRS and G _HRS across a matrix is largely influenced by the choice of the stack’s materials (e.g., single material HfO_X vs. bilayer HfO_X + TaO_X) [94,95], as well as by the device scaling. Extreme scaling seems to reduce variability likely because of a reduction of the area where the switching occurs [96]. Furthermore, the conductance value set during the CLT procedure also exhibits variability [19,93].

The normalized device variability can be expressed as σ/µ, where σ is the standard deviation and µ the mean value of the LRS and HRS conductance (G _LRS and G _HRS ) distributions. The DMM is suitable to study the device variability in terms of the G _LRS and G _HRS dispersion, as well as due to errors in the conductance programming, expressed as a variability of the DMM control parameter λ (see Figure 10A). To do so, we made a small amendment in the SPICE subcircuit representation of the DMM (see Table 1), and defined the parameters λ, I _min and I _max in terms of the gauss() SPICE function. For instance, for λ, we define H as param H = gauss(H_0,H_var,3), being H_0 and H_var two arguments passed to the subcircuit representing the nominal value of λ and its relative variability specified for 3 sigma, respectively. Similar variability is assumed for I _min and I _max parameters, which introduces dispersion in the G _LRS and G _HRS values. Because of the D2D variability among devices in the CPA, Monte Carlo (MC) simulations were performed, assuming different levels of variability for all devices. In each MC run, the characteristic (stored memory state for the resistive window) of each device was individually altered following a Gaussian distribution around the nominal value with increasing dispersion. Subsequently, the CPA was used to classify the images from the MNIST dataset. 10 MC runs were considered for each value of σ/µ.

In Figure 10B, the influence of the control parameter λ variability (σ _λ /µ _λ ranging from 0 to 30%) over the inference accuracy is studied for two different G _LRS /G _HRS ratios (10 and 100) achieved by two model plays (namely C1 and C2, shown in the inset of Figure 10B) and considering R _L = 10 Ω. The image size remains 8 × 8 px. and 4 partitions were used (that is, 4 CPAs of size 16 × 10 for each polarity of synaptic weights). For simplicity, no variability in the major I-V loop is considered (σ _GHRS = σ _GLRS = 0). Two trends are clearly observed. On one hand, the model play having a G _LRS /G _HRS ratio equal to or greater than 100 (C2) exhibits a very reduced sensitivity to λ variations (accuracy loss is below 5% for variabilities up to 30%). On the other hand, there is a sustained accuracy reduction for model play C1 over the same range of σ _λ /µ _λ . Given that the case of model play C2 is more robust to λ variations it was then selected to be thoroughly studied by considering the joint variability of G _LRS and G _HRS (σ _GLRS/µ _GLRS and σ _GHRS/µ _GHRS, respectively). As the variability is normally higher in HRS than in LRS [38], they were swept independently, resulting in the accuracy map illustrated in Figure 10C. The trend observed in Figure 10B is repeated in Figure 10C. Interestingly, σ _GHRS/µ _GHRS has a higher impact on the inference accuracy, likely due to a higher number of memdiodes mapped close to the G _HRS value.