Solid-State Electronics 185 (2021) 108040 Contents lists available at ScienceDirect Solid State Electronics journal homepage: www.elsevier.com/locate/sse SPICE modeling of cycle-to-cycle variability in RRAM devices E. Salvador a, *, M.B. Gonzalez b, F. Campabadal b, J. Martin-Martinez a, R. Rodriguez a, E. Miranda a a Departament d’Enginyeria Electrònica, Universitat Autònoma de Barcelona, 08193 Cerdanyola del Valles, Spain b Institut de Microelectrònica de Barcelona, IMB-CNM, CSIC, 08193 Cerdanyola del Valles, Spain A R T I C L E I N F O A B S T R A C T Keywords: In this work, we investigated how to include uncorrelated cycle-to-cycle (C2C) variability in the LTSpice quasi- RRAM static memdiode model for RRAM devices. Variability in the I-V curves is first addressed through an in-depth Variability study of the experimental data using the fitdistrplus package for the R language. This provides a first approxi­ Modeling mation to the identification of the most suitable model parameter distributions. Next, the selected candidate LTSpice distributions are incorporated into the model script and used for carrying out Monte Carlo simulations. Finally, Fitdistrplus the experimental and simulated observables (set and reset currents, transition voltages, etc.) are statistically compared and the model estimands recalculated if it is necessary. Here, we put special emphasis on describing the main difficulties behind this seemingly simple procedure. 1. Introduction underlying connection among their average and dispersion values. This interdependence must be addressed in order to obtain realistic simulated Memristors or Resistive RAMs (RRAMs) are electroformed metal- curves. In this work, we adopt a simplistic approach which consists in oxide-metal devices that can alter their resistance states in a non- estimating first the most suitable distribution functions for the experi­ volatile fashion. They are currently considered not only for memory mental observables and then attributing these functions to the corre­ arrays but also for a plethora of applications including neuromorphic sponding most relevant model parameters in each part of the I-V curve. computing, logic circuits, cryptography, etc. [1–3]. Nevertheless, one of However, recall that the final simulated curve is the result of a combi­ the major drawbacks that this technology currently faces concerns with nation of nonlinear functions with random parameters. In order to its inherent variability, which is associated with the alternate formation obtain acceptable results under this scheme, we follow a kind of recur­ and destruction of the filamentary conducting structure at atomic scale sive approach in which the model parameters are at the end selected [4,5]. Since variability is at the heart of RRAM operation, inclusion of according to the statistical results they generate. this phenomenon in any compact model intended for circuit simulation environments would be of utmost importance. In this work, we explore 2. Devices and experimental results the introduction of uncorrelated cycle-to-cycle (C2C) variability in the quasi-static memdiode model (QMM) for RRAM devices [6]. Correlated HfO2-based MIM structures [7] were investigated in this work. The C2C variability involving time series dependence (simulation results for oxide thickness is 10 nm and the area of the devices 5x5μm2. The bottom cycle i depend on simulation results for cycle i-1) and chained parameter electrode is a 200 nm-thick W layer and the top electrode is a 200 nm- simulations involving multivariate regression (parameter value i thick TiN layer on top of a 10 nm-thick Ti layer acting as oxygen getter explicitly depends on parameter value j within cycle k) are the following material [8]. Simulations are compared with experimental data ob­ logical steps to investigate after this first attempt. However, before tained from applying 450 voltage sweeps to the devices described pre­ entering into the subject, it is important to understand how variability viously. The complete set of experimental I-V curves and the median affects the simulation results. Uncorrelated C2C variability does not curve are shown in Fig. 1.a. High (IHRS) and low (ILRS) resistance state mean that the parameter estimands can be independently chosen within currents are extracted at V = 0.2 V. The snapback (SB) correction is a certain cycle. Since the generated I-V curve is the result of a sequence calculated as VSB = Vapplied − Ri ⋅Imeasured , where Ri is a series resistance of operations in which multiple parameters intervene, there is an (internal). Ri is chosen so as to achieve a current increase at a constant * Corresponding author. E-mail address: emili.salvador@uab.cat (E. Salvador). https://doi.org/10.1016/j.sse.2021.108040 Available online 29 May 2021 0038-1101/© 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). E. Salvador et al. Solid State Electronics 185 (2021) 108040 Table 1 Goodness-of-fit statistics (Kolmogorov-Smirnov, Cramer-von Mises, Anderson- Darling) and criteria (Akaike’s and Bayesian) for a) VT , b) for VR and c) for RI . 3. QMM with variability Once the model parameter distributions are established, they are included in the header of the QMM model script as shown in Table 2. Fig. 1. a) Experimental I-V curves: 450 cycles and median curve. b) Voltage Some model parameters are assumed to be constant in order to avoid correction illustrating the snapback effect in a single cycle and parameter over-randomness. The LTSpice gauss function and its transformations specification. are used to generate the appropriate parameter values and variability. According to the QMM [2] the I-V characteristic reads: transition voltage VT (see Fig. 1.b) [9]. This correction is calculated for I(V) = I0 (λ)sinh{α(λ)[V − (RS (λ) + Ri )I ] } (1) all the cycles using a different Ri value. The RESET voltage VR is the voltage corresponding to the maximum current reached at negative bias. where I0 (λ) = Iomin + (Iomax − Iomin )⋅λ is the current amplitude factor, RS a The SB current ISB is found from the last data point before the occurrence series resistance, α a fitting parameter and Iomin and Iomax , minimum and of the SET jump. The parameter extraction was carried out using maximum current values, respectively. Equation (2) relates the memory MATLAB. The extracted parameters for the whole set of curves were state λ to the voltage across the filament’s constriction VC through the analyzed using the tools available in the fitdistrplus package for the R recursive hysteresis operator: language [10]. Different distributions (normal, lognormal, gamma, and { [ ← ]} Weibull) were fitted to the data and compared. The goodness of fit and λ(VC ) = min Γ− (VC ), max λ(VC ), Γ+ (VC ) (2) criteria (minimum values in Table 1) in combination with the quanti­ le–quantile (Q-Q) plots are illustrated in Fig. 2, which graphically where Γ+ and Γ− are the so-called ridge functions, which represent the compare the experimental distributions with the parametric models. (←) This information is used to determine the best candidate distributions. ion/vacancy movement for SET and RESET, respectively. λ VC is the Fig. 2.a, .b and .c show the density plots for VT, VR memory value a timestep before (dictated by the simulator). The model andRi corresponding to the four distributions investigated. According to includes other parameters for the fine-tuning of the simulated curves. this analysis, VT follows a lognormal distribution while VR and Ri are better described by a normal one. The same analysis was carried out for 4. Comparison between experimental and simulation results the rest of the above mentioned parameters resulting in: normal distri­ bution for ISB, and lognormal distributions for IHRS and ILRS. This prior The next step consists in performing Monte Carlo simulations using information will be used to define the distribution of the leading model the QMM model discussed in the previous Section. In Fig. 3, experi­ parameters in the device script. mental (grey) and simulated (red) curves are compared using three 2 E. Salvador et al. Solid State Electronics 185 (2021) 108040 Fig. 2. a), b) and c) show density plots for VT , VR and Ri , respectively, with normal, gamma, lognormal, and Weibull distributions. d), e) and f) show Q-Q plots after selecting the best fitting distribution. different plots for the I-V curves: a) linear–linear, b) log-linear, and c) specifically selected so as to compile information from the entire set of linear–linear (with SB correction). Although the fitting is not perfect, the the I-V curves in the different regimes: IHRS for the HRS curve, VT for the simulations reproduce the main features of the experimental curves in SET transition, ILRS for the LRS curve and VR for the RESET transition. all the cases. The most conspicuous difference occurs in LRS, mainly This compilation was performed without considering the SB correction because of a peculiarity of the experimental data. This will be discussed and for a direct evaluation of IHRS and ILRS . As shown in Fig. 4, simulation in detail next. results for IHRS , VT , and VR are consistent with the experimental data. In order to make a fair comparison, the statistical distributions of the However, ILRS exhibits a notorious discrepancy. The reason is clear, for four observables is considered in Fig. 4. These observables were the particular case under study, the experimental ILRS histogram presents 3 E. Salvador et al. Solid State Electronics 185 (2021) 108040 Table 2 LTSpice QMM script including variability in the most significant parameters taking into account the previously extracted best candidate distributions. In red the Gaussian distribution and in blue the lognormal distributions. two peaks which likely corresponds to two different configurations of the filamentary path or to competing filaments. The simulated curves are unable to capture this feature because of the use of a single distri­ bution function. This single distribution is responsible for the large central peak shown in Fig. 4.c. This inconsistency illustrates that caution should be exercised when unexpected deviations in the experimental data occur. Finally, Fig. 5 illustrates the importance of C2C correlations in the experimental curves [11]. Experimental and simulated data corre­ sponding to the time evolution of two parameters (high resistance state current IHRS and reset voltage VR) are compared. Notice that our approach is unable to generate a trend since the LTSpice simulations considered in this work are uncorrelated. These issues will be further investigated in the future. 5. Conclusions In this work, we explored how uncorrelated C2C variability can be included in the QMM model for resistive switching devices. The exper­ imental I-V curves obtained from HfO2-based devices were analyzed and the main parameters extracted. The best candidate distributions for the experimental parameters were determined using the fitdistrplus package from the R language and this information was included in the model Fig. 3. Experimental and simulated curves comparison. a) I-V curves, b) I-V script. We showed that the model results reproduce reasonably well the curves in logscale and c) I-V curves applying the SB voltage correction. 4 E. Salvador et al. Solid State Electronics 185 (2021) 108040 Fig. 4. Comparison of experimental and simulated parameter distributions: a)IHRS , b) VT , c)ILRS , and d)VR . main features exhibited by the experimental curves. Future in­ vestigations will include time series and chained-parameter analysis for correlated C2C variability simulations. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work was supported by the Spanish Ministry of Science, Inno­ vation and Universities through projects TEC2017-84321-C4-1-R, TEC2017-84321-C4-4-R, and PID2019-103869RB-C32. References [1] Yu S, et al. IEEE Solid State Circuits Magazine 2016;8(2):43–56. [2] Lee JS, Appl., et al. Phys. Rev. 2015;2(3):031303. [3] Jo Sung Hyun, et al. Nano Lett 2010;10(4):1293–301. [4] Chen A, et al. IEEE IRPS 2011;843–846. [5] Bargallo M, et al. IEEE Trans. Electron Devices 2016;63(8):3116–22. [6] Miranda E. IEEE TNano 2015;14:787–9. [7] Piccolboni G, et al. IEEE Electron Device Lett. 2016;37(6):721–3. [8] Poblador S, et al. Microelectron. Eng. 2018;187–188:148–53. [9] Karpov V, et al. Appl. Phys. Lett. 2015;109:1–5. [10] Delignette-Muller ML, et al. J. Stat. Softw. 2015;64:1–34. [11] Alonso FJ, et al. Chaos. Solitons and Fractals 2021;143:110461. Fig. 5. Time evolution of experimental and simulated parameters for: a) HRS current and b) reset voltage. 5 E. Salvador et al. Solid State Electronics 185 (2021) 108040 Emili Salvador Aguilera is a PhD student at the Electronical Javier Martin-Martinez received the M.S. degree in physics Engineering department from the Universitat Autònoma de from the Universidad de Zaragoza, Zaragoza, Spain, in 2004, Barcelona under the supervision of Enrique Alberto Miranda and the Ph.D. degree from the Universitat Autònoma de Bar­ and Rosana Rodríguez, from the same department. Prior to the celona (UAB), Bellaterra, Spain, in 2009. He was with the PhD, coursed my Bachelor degree in Nanoscience and Nano­ Università degli Studi di Padova, Padua, Italy, and IMEC, technology and Master degree in Advanced Nanoscience and Leuven, Belgium. He is currently an Associate Professor with Nanotechnology at the Universitat Autònoma de Barcelona. He UAB. His main research interests include the characterization carried out my Master thesis at the Institut Català de Nano­ and modeling of failure mechanisms in MOSFETs and also ciència i Nanotechnology (ICN2) working in the graphene RRAM characterization and modeling for neuromorphic transistors field. applications. Mireia Bargallo Gonzalez received the degree in physics from Rosana Rodriguez received the Ph. D. in Electrical Engineer­ the University of Barcelona, Barcelona, Spain, and the Ph.D. ing from Universitat Autonoma de Barcelona (UAB) in 2000. degree on the topic of stress analysis and defect characteriza­ Funded by the Fulbright program, she worked on devices and tion techniques of semiconductor materials and devices, from circuits reliability at the IBM Thomas J. Watson Research Katholieke Universiteit Leuven, Leuven, Belgium, in 2011. She Center (USA). Currently, she is associate professor at the UAB. pursued her Ph.D thesis with the Interuniversity Microelec­ Her research is focused on the variability and reliability of tronics Center (imec), Leuven, Belgium. In 2011, she joined the advanced CMOS devices. She is interested in the electrical Institut de Microelectrònica de Barcelona (IMB-CNM, CSIC). characterization and modeling of process-related and time­ Her current research interests include electrical characteriza­ –dependent variability sources as Random Telegraph Noise tion, modeling and applications of resistive switching devices. (RTN) and aging mechanisms as Bias Temperature instability (BTI) and Hot Carrier Injection (HCI). Her research includes the study of the variability impact on the performance of single devices and digital and analogical circuits. She is also interested in the characterization of resistive switching devices (memristors) and their application for non-volatile memories, computing and neuromorphic applications. Francesca Campabadal received the Ph.D. degree in physics from the Universitat Autònoma de Barcelona, Bellaterra, Spain, Enrique Miranda is Professor at the Universitat Autonoma de in 1986. She joined the Institut de Microelectrònica de Barce­ Barcelona (UAB), Spain. He has a PhD in Electronics Engi­ lona, Consejo Superior de Investigaciones Científicas, Barce­ neering from the UAB (1999) and a PhD in Physics from the lona, Spain, in 1987, where she is currently a Research Universidad de Buenos Aires, Argentina (2001). He received Professor. Her current research interests include the deposition numerous scholarships and awards including: RAMON y CAJAL of high-k dielectric layers, their electrical characteristics, and (UAB), DAAD (Technical University Hamburg-Harburg), the resistive switching phenomena in RRAM devices. MATSUMAE (Tokyo Institute of Technology, Japan), TAN CHIN TUAN (Nanyang Technological University, Singapore), WALTON award from Science Foundation Ireland (Tyndall National Institute), Distinguished Visitor Award (Royal Acad­ emy of Engineering, UK), CESAR MILSTEIN (CNEA, Argentina), Visiting Professorships from the Abdus Salam International Centre for Theoretical Physics, Slovak Academy of Sciences, Politecnico di Torino, Leverhulme Trust (University College London, UK), and Nokia Foundation (University of Turku, Finland). He serves as member of the Distinguished Lecturer program of the Electron Devices Society (EDS-IEEE) since 2001 and as Associate Editor of Microelectronics Reliability since 2003. He has authored and co-authored around 250 peer-review journal papers. 6