Versatile Filamentary Resistive Switching Model
Iosif-Angelos Fyrigos, Vasileios Ntinas, Georgios Ch. Sirakoulis, Panagiotis Dimitrakis, Ioannis G. Karafyllidis
11 Versatile Filamentary Resistive Switching Model
Iosif-Angelos Fyrigos,
Member, IEEE,
Vasileios Ntinas,
Member, IEEE,
Georgios Ch. Sirakoulis,
Member, IEEE,
Panagiotis Dimitrakis,
Senior Member, IEEE, and Ioannis G. Karafyllidis
Abstract —Memristors as emergent nano-electronic deviceshave been successfully fabricated and used in non-conventionaland neuromorphic computing systems in the last years. Severalbehavioral or physical based models have been developed toexplain their operation and to optimize their fabrication pa-rameters. All existing memristor models are trade-offs betweenaccuracy, universality and realism, but, to the best of ourknowledge, none of them is purely characterized as quantummechanical, despite the fact that quantum mechanical processesare a major part of the memristor operation. In this paper,we employ quantum mechanical methods to develop a completeand accurate filamentary model for the resistance variationduring memristor’s operating cycle. More specifically, we applyquantum walks to model and compute the motion of atomsforming the filament, tight-binding Hamiltonians to capture thefilament structure and the Non-Equilibrium Green’s Function(NEGF) method to compute the conductance of the device.Furthermore, we proceeded with the parallelization of the overallmodel through Graphical Processing Units (GPUs) to accel-erate our computations and enhance the models performanceadequately. Our simulation results successfully reproduce theresistive switching characteristics of memristors devices, matchwith existing fabricated devices experimental data, prove theefficacy and robustness of the proposed model in terms of multi-parameterization, and provide a new and useful insight into itsoperation.
Index Terms —Filamentary Resistive Switching, Memristor,Modeling, Non-Equilibrium Green’s Function, Quantum Walks.
I. I
NTRODUCTION
The main resistance switching (RS) mechanisms in metal-oxide-metal (MOM) memristive devices are very well de-scribed in [1], [2] and can be categorized [3] as follows: 1.electrochemical metallization bridge (EMB), 2. metal oxide- bipolar filamentary (MOxBF), 3. metal oxide- unipolarfilamentary (MOxUF) and 4. metal oxide - bipolar nonfila-mentary (MOxBN). A more compact classification based onthe physical mechanisms responsible for the nanoionic natureof the resistive switching [4] includes: the ElectrochemicalMemories (ECM), the Valence Change Memories (VCM) andthe Thermochemical Memories (TCM) [5], [6]. For the sakeof comparison, ECM and EMB are the same, VCM comprisesMOxBF and MOxBN while TCM is similar to MOxUF. Forall these mechanisms the macroscopic image of the resistanceswitching is the formation of a conductive filament (CF) thatshorts the two electrodes in MOM memory cells, except the
I.-A. Fyrigos, V. Ntinas, G. Ch. Sirakoulis and I.G. Karafyllidis are with theDepartment of Electrical and Computer Engineering, Democritus Universityof Thrace, Xanthi, Greece.V. Ntinas is also with the Department of Electronics Engineering, Univer-sitat Polit´ecnica de Catalunya, Barcelona, Spain.P. Dimitrakis is with Institute of Nanoscience and Nanotechnology, NCSR“Demokritos”, Athens, Greece. case of MOxBN where the transformation of the metal-oxidecontact into Schottky occurs. Obviously, MOxBN is the onlyone where RS depends on the area of the electrodes. Theconductive filament may be attributed to metal atoms due todrift/diffusion of metal cations from one electrode to the other(ECM) or to oxygen vacant sites (metal atoms in the oxidebulk with unsaturated bonds) due to drift/diffusion of oxygenvacancies which can be assumed as anions.On the other hand, several models have been developed todescribe and simulate the operation of resistance switchingdevices according to the above physical mechanisms. One ofthe main goals of these models is to describe the phenomenonboth quantitatively and qualitatively. Ion drift and diffusionequations are used to simulate the kinetics of ions in the insu-lator which should be considered in their time and temperaturedependent forms. The most common approach to solve thesystem of these equations is the finite element analysis (FEA)and kinetic Monte Carlo (kMC) [7], [8]. In order to tackle thephysics and the computation complexity of the microscopicmodels, other macroscopic models have been adopted [9],[10], [11], [12], [13]. Moreover, various heuristically derivedwindow functions and mathematical techniques have beendeveloped to address the reduced accuracy of the memristormodel near the boundaries of devices kinetics [14], [15], [16],[17], [18]. In parallel, phenomenological models based on theexperimental measurements of the RS device behavior havebeen developed in an efficient way to enable the memristorcircuit design [19], [20], [21], compensating accuracy lossand generality with efficient circuit simulations. In brief, themicroscopic models are necessary for use in manufacturersprocess development kit (PDK) and macroscopic models forcircuit simulations.All of the microscopic models described above do not usequantum mechanical methods to describe and simulate thememristor operation, although both the ionic kinetics andconduction phenomena in the nanoscale can be describedaccurately using quantum mechanical methods [22]. In thiswork, a filamentary based nanoscale memristor model forECM has been developed. The model assumes a standardmetal-insulator-metal RS device in 2D. The ion kinetics in thedielectric that results in filament formation and deformationare modeled using quantum walks, which incorporate theapplied potential distribution in the filament.The structure of the filament in the dielectric is describedusing tight-binding Hamiltonians. The electron transport in thenanoscale filament/dielectric system is modeled using Non-equilibrium Green’s Functions (NEGF) [22]. Tight-bindingHamiltonians combined with NEGF are used to compute theconductance of the nanoscale memristor as a function ofthe applied voltage. The process of the filament formation a r X i v : . [ c s . ET ] A ug and dissolution in the solid dielectric is modeled by dividingthe process into computational time steps. At each step thenew filament structure is simulated using quantum walks andthe Hamiltonian describing the system is constructed. ThisHamiltonian enters the NEGF method to compute the deviceconductance and current for this step. Our results reproducesuccessfully the I-V characteristics of RS devices and showthat the conductance quantization and the corresponding mul-tiple RS states become more prominent as the solid dielectricthickness decreases.The model is fully parametrized and ready for calibrationusing experimental data. As a proof of concept, experimentalresults of a fabricated metal-oxide-metal memristive devicewere used to prove the quantitative and qualitative agreementof the presented simulation results with the fabricated device’sones. Specifically, we fabricated and measured planar MIMdevices made of 30nm Cu (active electrode) on a 40nm thickSiO deposited by sputtering and 100nm thick Pt (counterelectrode). The whole MIM structure was sited on a Ti(5nm)/thermal-SiO (300nm)/n-Si (P-doped, 1-5 Ω cm). Current -voltage sweeps were performed in order to investigate the SET and
RESET switching processes. The experimental re-sults suggested a clear bipolar resistance switching operationwith
SET and
RESET voltages 4.5V and -4.5V, respectively.The Cu/SiO is a well-known combination to achieve ECM(MOxBF) ReRAM devices [23]. The resistance switching isattributed to the oxidation of Cu atoms lying at the interfacewith the oxide and the resulted Cu ions are drift and diffuseacross the bulk oxide layer towards the counter electrode. Thisprocess has been experimentally proven and modeled by manyinvestigators [6], [24], [25], [26], [27]. Furthermore, our modelcan reproduce successfully the resistive switching characteris-tics of MOxBF RS devices under the appropriate modificationof specific parameters. Such a quantum mechanical model canbe used to simulate the operation of nanoscale memristors andprovides a new and useful insight into their operation.II. M ODEL S TRUCTURE AND O UTLINE
During the memristor operation a filament is progressivelyformed in the dielectric, between the active electrode andthe counter electrodes. The process of filament formation issimulated in discrete computational time steps. At each stepthe applied voltage at the device electrodes is increased by aconstant value that is the voltage step of the applied voltagesweep. In this way, the dependence of the filament formationprocess on the voltage sweep rate can be investigated. Figure1 shows the progress of filament formation in the memristorat two successive computational time steps.The filament structure at computational time step s is shownin the upper part of Fig. 1. The simulation for computationaltime step s +1 involves the following computation steps, whichare briefly provided here and will be thoroughly explained laterthrough the corresponding sections:1) The structure of the filament and the potential distri-bution in the filament at step s , are the inputs for thecomputations for computational step s + 1 . Fig. 1. Filament formation in the dielectric at two successive computationaltime steps, s and s + 1 .
2) The filament structure and potential distribution in thefilament at computational time step s , enters the quan-tum walk algorithm, in a manner to be explained lateron. The quantum walk algorithm computes the kineticsof the metal ions and outputs the filament structure atcomputational time step s + 1 , shown in the lower partof Fig. 1.3) The applied voltage is increased by a constant valueand the potential distribution in the filament structure atcomputational time step s + 1 is computed.4) The tight-binding Hamiltonian for the dielectric-filament system at computational time step s + 1 iscomputed. The computation of the Hamiltonian will beexplained later on.5) The potential distribution at computational time step s +1 and the Hamiltonian enters the NEGF method, in amanner explained later on, and the conductance of thememristor device at time step s + 1 is computed.6) The potential distribution and structure of the filamentcomputed for computational time step s + 1 , are used asinputs for the next computational time step, s + 2 .As the applied voltage increases, the filament progressivelyis formed until it reaches the active electrode. After that,the applied voltage is decreased, the metallic ions dissolvein the dielectric and move towards the active electrode, sothe filament deforms. Again the quantum walk is used tosimulate the kinetics of the filament deformation. The modeldescribed briefly above computes the memristor conductanceas a function of the applied voltage. In the next sections,as already mentioned, the application of the quantum walkmodel and the NEGF method in the specific problem will bedescribed analytically. III. C
ONDUCTANCE C ALCULATION OF D EVICE U TILIZING
NEGF M
ETHOD
The memristor dielectric is divided into a matrix of atomic-size square cells. During the memristor operation these cellsare occupied either by atoms of the dielectric or metallic ions,diffusing from the active electrode. Electrons are transportedfrom one electrode to the metallic ions in the dielectricand from there to the other electrode. In the tight-bindingapproach electrons are transported stepwise from a cell to itsneighboring cells at each computational time step. This processis described by the tight-binding Hamiltonian, H D : H D = (cid:80) i,j k,l =+1 (cid:80) k,l = − ( − t i + k,j + l )ˆ c † i + k,j + l ˆ c i + k,j + l ++ (cid:80) i,j k,l =+1 (cid:80) k,l = − ˆ c † i + k,j + l ˆ c † i + k,j + l V ˆ c i + k,j + l ˆ c i + k,j + l (1)In the general case, an electron is located at the ( i, j ) cell.The electron motion is described using fermion annihilationand creation operators, namely c and c † , which annihilate andcreate electrons in the ( i, j ) cell and its neighbors ( i +1 , j ) , ( i − , j ) , ( i, j + 1) and ( i, j − . Electron motion is associatedwith a kinetic energy t i,j , which is the overlap integral. Theoverlap integral between non-neighboring cells is zero. Theoverlap integral also describes the effect of temperature onelectron mobility. The second term of (1) is the retentionenergy associated with the tendency of electrons to remainin the current cell, the potential in which is described by V .The conductance of the memristor for various electronenergies, G ( E ) , is given by the matrix equation [28]: G ( E ) = 2 q ¯ h T race [Γ G D Γ G A ] , (2) G D is the retarded Green’s function, and G A is the advancedGreen’s function. Γ and Γ are the broadening functions ofthe two contacts: Γ = i [Σ − Σ +1 ] , Γ = i [Σ − Σ +2 ] . (3)The contacts are described by their self-energies, Σ and Σ ,respectively. The retarded Green’s function is given by [29]: G D = [ E − H D − Σ − Σ ] − , (4)while the advanced Green’s function is: G A = ( G D ) + . (5)Fig. 2 depicts the outline of NEGF method applied to thememristor device. The device is described by the retardedGreen’s function, which determines its conductance through(5) and (2). The tight-binding Hamiltonian H D , describes theelectron transport in the dielectric-metallic ion system, the self-energy Σ the electron transport between the active electrodeand the dielectric, and the self-energy Σ the electron transportbetween the counter electrode and the dielectric. Fig. 2. The parts of the matrix eq. 4 for retarded Green’s function G D thatdescribe the active electrode, the dielectric containing metallic ions and thecounter electrode, respectively. A. GPU Acceleration of NEGF
One of the major issues resulting from the utilizationof the NEGF method in the application to the memristiveconductance calculations is the request for high performancecomputing resources. In order to tackle this issue, Graphicsprocessing units (GPUs) originally designed for computervideo cards but recently emerged as the most powerful chipin a high-performance workstation, are exploited to boost theperformance of the NEGF method. Unlike multicore commer-cial CPU architectures, GPU architectures contains thousandof cores capable of running numerous threads in parallel. Asa result, GPU by offering an almost ceaselessly increasednumber of cores couple with a high memory bandwidth,provides incredible computational resources that enable it as asuitable candidate to accelerate NEGF method’s computation.More specifically, NEGF method, when employed to com-pute the conductance of nanodevices, is a computationaldemanding technique and conductance becomes harder to becomputed as the number of atoms that constitute the deviceincreases. NEGF method parallelization suitability arises fromthe fact that it contains operations between large matriceswhich grow exponentially in size as the number of simulatedatoms increases. Moreover, the different energy levels ( E )for which the conductance ( G ( E ) ) has to be computed arecompletely independent from each other, meaning that thecomputation of each energy level can be allocated to a differentmultiprocessor of the GPU. Therefore, NEGF method wasaccelerated through GPU parallelization. More specifically,two different algorithms were developed, one for GPUs thatare optimized for single precision operations and another, moreprecise, for GPUs that are optimized for double precisionoperations, each one suitable for different hardware implemen-tations of GPUs. The results of the comparison as depictedin Fig. 3, between a up-to-date CPU and a correspondingGPU devices, are very promising showcasing that the NEGFmethod scales almost linearly in GPU while in the CPU thereis an exponential increase of execution time. As a result, theattainable speedup factor of the GPU compared to the CPUwas more than × . Grid Size E x e c u t i on T i m e ( s ) CPUGPU
Fig. 3. Comparison between CPU (i9-9900K) and GPU (2060 Super) duringexecution of NEGF method for different number of atomic grids
IV. M
ODELLING OF F ILAMENTARY G ROWTH S UB - SYSTEMUTILIZING Q UANTUM W ALKER ’ S PROBABILITYDISTRIBUTION
Due to the nanoscale dimensions of the memristors, thefabricated devices are governed by several quantum phenom-ena (quantum tunneling, Poole Frenkel, etc.). To integrate thequantum physics into the modelling of the filament, quantumwalks [30], the quantum mechanical counterpart of classicalrandom walks where utilized. The probability of each atomicslot to be occupied by a Cu or Di atom, arise from themotion of the quantum walker. More specifically, in general,any quantum walker is characterized by its position in the 2Dspace at any given moment as well as by its quantum state.In our realization the quantum state (coin state) of the walkeris described using 2 qubits. | c (cid:105) = a | (cid:105) + b | (cid:105) + c | (cid:105) + d | (cid:105) , (6) | a | + | b | + | c | + | d | = 1 , (7)where a , b , c and d are the probability amplitudes, the absolutesquare of which gives the probability of the walker to be inthe corresponding state.Coin state | c (cid:105) determines the state of the quantum walkerin a specific location. Every position of the 2D grid has astate vector. In the presented case, Fig. 4 shows the 2D squaregrid into which the dielectric-filament system is divided.Yellow cells represent dielectric atoms and blue cells representmetallic ion atoms. In discrete 2D quantum walk, the walker(metallic ion atom) moves on the lattice formed by connectingthe square cell centers, as shown in Fig. 4 by red dashed lines.The state of the quantum walker at computational time step s , namely | S s (cid:105) , is its position on the lattice at this step: | S t (cid:105) = | i, j (cid:105) . (8)The Hilbert space of the discrete quantum walk comprisestwo subspaces, the location subspace H L , which is spannedby the basis: Fig. 4. The grid of square cells and the quantum walk lattice. Yellow cells areoccupied by dielectric atoms while blue cells by metallic ions. The quantumwalk lattice is formed by connecting the centers of the square cells. | S (cid:105) = | − n, n (cid:105) , L | i − , j (cid:105) , | i, j (cid:105) , | i + 1 , j (cid:105) , | i, j + 1 (cid:105) , L | n, n (cid:105) , (9)and a four-dimensional coin subspace, H C , which is spannedby the four coin basis states | (cid:105) , | (cid:105) , | (cid:105) and | (cid:105) . TheHilbert space, namely H , of the quantum walk is the tensorproduct ⊗ of these two spaces: H = H L ⊗ H C . (10)The state of the quantum walker at computational time step s is the tensor product of its location and its coin state: | S t (cid:105) = | location at time s (cid:105) ⊗ | coin at time s (cid:105) . (11)For example, the quantum walker’s state at location ( i, j ) withcoin in state | (cid:105) is | i, j (cid:105) ⊗ | , (cid:105) . During the evolution of thequantum walk, the new state of the quantum walker at the nextcomputational time step, s + 1 is given by: | S s +1 (cid:105) = U ( s ) | S s (cid:105) , (12)where U ( s ) is a unitary operator that drives the quantum walkand it is the product of the following operators: U ( s ) = P ( s ) C ( s ) exp (cid:16) i ¯ h V ( i, j, t ) (cid:17) . (13)In the above equation the symbol i , which is divided by thePlank constant ¯ h is the imaginary unit, while V ( i, j, s ) is thepotential at the ( i, j ) cell at step s . The quantum walk modelsthe kinetics of metallic ions in the following way: First thecoin operator, C ( s ) , acts on the coin state, then the quantumwalker (metallic ion) moves to one of the neighboring cellsas a result of the shift operator, P ( s ) , acting on the previouslocation state and finally the product of these two operators ismultiplied by the exponent of (13), through which the potentialdistribution enters the computation. The coin operator C ( s ) isthe tensor product of two Hadamard quantum gates and by Filament evolution
Quantum walk probability distribution (a)
Filament evolution
Quantum walk probability distribution (b)
Filament evolution
Quantum walk probability distribution (c)
Filament evolution
Quantum walk probability distribution (d)Fig. 5. Filament evolution due to Quantum Walker’s probability distribution. (a) Initial filament growth with positive applied voltage. (b) Final filamentformation. (c) Degradation of the filament due to the applied negative voltage. (d) Termination of the filament acting on the current coin state brings it to a superposition ofthe four basis coin states: C ( s ) | coin at time s (cid:105) = c R | (cid:105) + c L | (cid:105) + c U | (cid:105) + c D | (cid:105) . (14)In the above equation, c R , c L , c U and c D are the probabilityamplitudes, the absolute square of which gives the probabilityof the walker to move to the right, left, up and down,respectively. The shift operator acts on the location state andmoves the quantum walker to one of the neighboring cellswith probability determined by the probability amplitudes ofeq. (14) as follows: P ( t ) = i (cid:88) j = − n | i + 1 , j (cid:105) , (cid:104) i, j | ⊗ | (cid:105)(cid:104) | + | i − , j (cid:105) , (cid:104) i, j | ⊗ | (cid:105)(cid:104) | + | i, j + 1 (cid:105) , (cid:104) i, j | ⊗ | (cid:105)(cid:104) | + | i, j − (cid:105) , (cid:104) i, j | ⊗ | (cid:105)(cid:104) | . (15)To compute the next state of the quantum walk operator,(13) is applied to all cells of the 2D grid. The exponentialfactor containing the potential determines the constructive ofdestructive interference of the quantum walkers arriving ateach cell, and the new distribution, at computational time step s + 1 , of the metallic ions into the dielectric is determined.V. S IMULATION R ESULTS OF Q UANTUM -W ALK BASED F ILAMENT FORMATION
By executing the quantum walk algorithm the probabilitiesof the walker positions start to spread. In every computational time step of the evolution of the walker, a thresholding checkis performed. The grid’s locations whose probabilities surpassthis threshold are occupied by a metallic ion. Consequently,a filament based on the evolution of the quantum walker isformed like the one in Fig. 5(a,b). During negative appliedvoltage, two quantum walkers are initialized to the upper sidesof the grid to simulate the degradation of the filament as shownin Figs. 5(c,d), respectively. More specifically, due to the nega-tive applied voltage, the probability of the backwards movingquantum walkers determines the probability of the metallicions to move away from the previously formed filament. InFig. 5(c) it can be observed that while quantum walker’sprobability distribution surrounds the filament, metallic ionsstart to get removed and the filament becomes thinner.In Fig. 6 the logarithmic I-V curve for a different numberof dielectric Retention Energies, namely
Eodi , and in Fig. 7of metal Retention Energies, namely Em , are displayed,respectively. By examining Figs. 6, 7 it can be observedthat when the dielectric retention energy is increased, the R off of the device becomes higher (vertically widening theI-V), while increasing the Metal retention energy both R off and R on increases (shifting the whole I-V down). These twovariables can manipulate the electron motion and therefore,act as fitting parameters to regulate the R off and R on statesof the fabricated devices.Finally, as a proof of concept of the proposed model abilitiesto simulate the operation of a filamentary resistive switchingmemristive device, a fitting of the presented model with afabricated device has been conducted. As mentioned before,we fabricated and measured planar MIM devices made of30nm Cu (active electrode) on a 40nm thick SiO deposited by -5 0 5 Voltage (V) -10 -5 C u rr en t ( A ) Eodi =10 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Eodi =10.4 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Eodi =10.8 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Eodi =11.2 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Eodi =11.6
Fig. 6. Logarithmic I-V for a range of dielectric’s retention energy (
Eodi ) between 10eV and 11.6eV utilizing Quantum Walk model -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Em =7.1 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Em =7.5 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Em =7.9 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Em =8.3 -5 0 5
Voltage (V) -10 -5 C u rr en t ( A ) Em =8.7
Fig. 7. Logarithmic I-V for a range of metal’s retention energy ( Em ) between 7.1eV and 8.7eV utilizing Quantum Walk model -6 -4 -2 0 2 4 6 Voltage (V) -10 -5 C u rr en t ( A ) SimulationExperimental
Fig. 8. Comparison between experimental data and simulation output. sputtering and 100nm thick Pt (counter electrode). Current -voltage staircase sweeps were performed in order to investigatethe
SET and
RESET switching processes. The experimentalresults suggested a clear bipolar resistance switching operationwith
SET and
RESET voltages 4.5V and -4.5V, respectively.The Cu/SiO is a well-known combination to achieve ECM(MOxBF) ReRAM devices [23]. The resistance switching isattributed to the oxidation of Cu atoms lying at the interfacewith the oxide and the resulted Cu ions are drifted and diffuseacross the bulk oxide layer towards the counter electrode.More specifically, in Fig. 8, the comparison between the modeland the experimental data can be thoroughly examined. Theoscillation of the current during negative applied voltage isattributed to the scattering of the electrons due to disturbanceof the material’s lattice uniformity. This phenomenon is notpresent to the experimental data because the measuring setuptakes multiple samples of the same voltage and, consequently,leading to a smoother I-V characteristic. Specifically, the I-V measurements were experimentally recorded using a semi-conductor parameter analyzer HP4155A with medium integra-tion time, that means the current was averaged over 1msec.A valuable further analysis that should be conducted isthe time analysis of the filament’s evolution. Evaluating thetime interval of the filaments evolution parameter ( dt ) isrequired to specify the timing and frequency derived for themodel, taking under consideration some physical variables incorrespondence to the simulations carried out, such as thetemperature distribution of the device, the electron mobility aswell as the ion concentration. The ion mobility correspondsto the kinetic energy t found in (1), the ion concentrationis linked to the number of atoms (size of the Hamiltonian)that constitute the device, as well as, to the quantum walker’sdistribution which is affected by the applied voltage and theinternal temperature of the device is linked to the overlapintegral. VI. C ONCLUSIONS
In this work a novel quantum based model for filamentaryresistive switching nanodevices, like MOxBF RS nanodevices,was presented. A quantum random walk based algorithm wasdeveloped to model the evolution of the conductive filament,while in each step of the filaments evolution the conductanceof the device is calculated by the GPU accelerated NEGFmethod. The simulation results provide clear identifications ofthe robustness of the model for a number of various retentionenergies for both dielectric and metal. Furthermore, as asingle proof of concept, experimental data were used froma fabricated resistive switching device, namely planar MIMdevices made of 30nm Cu (active electrode) on a 40nm thickSiO deposited by sputtering and 100nm thick Pt (counterelectrode) and were compared with the model’s output whenfitted to much the experimental data. The presented resultsshowed a fair qualitative and quantitative agreement betweenthe experimental and simulated data. As a further step, moreparameters like temperature distribution, electron mobility andion concentration of the under study devices will be consid-ered. Finally, the GPU acceleration will be improved utilizingcuBLAS and cuSPARSE libraries to embed the efficiency ofthe sparse matrices into the GPU.A CKNOWLEDGMENT
We gratefully acknowledge the Greece-Russia bilat-eral joint research project MEM-Q (proj.no./MIS T4 ∆ P Ω -00030/5021467) supported by GSRT and funded by Nationaland European funds. R EFERENCES[1] L. Chua, G. C. Sirakoulis, and A. Adamatzky,
Handbook of MemristorNetworks . Springer International Publishing, 2019.[2] P. Dimitrakis and E. Valov,
Metal Oxides for Non-volatile Memory .Elsevier, 2021.[3] A. H. Edwards, H. J. Barnaby, K. A. Campbell, M. N. Kozicki, W. Liu,and M. J. Marinella, “Reconfigurable memristive device technologies,”
Proceedings of the IEEE , vol. 103, no. 7, pp. 1004–1033, 2015.[4] M. Lanza, H.-S. P. Wong, E. Pop, D. Ielmini, D. Strukov, B. C.Regan, L. Larcher, M. A. Villena, J. J. Yang, L. Goux, A. Belmonte,Y. Yang, F. M. Puglisi, J. Kang, B. Magyari-Kpe, E. Yalon, A. Kenyon,M. Buckwell, A. Mehonic, A. Shluger, H. Li, T.-H. Hou, B. Hudec,D. Akinwande, R. Ge, S. Ambrogio, J. B. Roldan, E. Miranda, J. Sue,K. L. Pey, X. Wu, N. Raghavan, E. Wu, W. D. Lu, G. Navarro, W. Zhang,H. Wu, R. Li, A. Holleitner, U. Wurstbauer, M. C. Lemme, M. Liu,S. Long, Q. Liu, H. Lv, A. Padovani, P. Pavan, I. Valov, X. Jing,T. Han, K. Zhu, S. Chen, F. Hui, and Y. Shi, “Recommended methodsto study resistive switching devices,”
Advanced Electronic Materials ,vol. 5, no. 1, p. 1800143, 2019.[5] I. Valov, R. Waser, J. R. Jameson, and M. N. Kozicki, “Electrochemicalmetallization memories—fundamentals, applications, prospects,”
Nan-otechnology , vol. 22, no. 25, p. 254003, may 2011.[6] I. Valov, “Redox-based resistive switching memories (rerams): Electro-chemical systems at the atomic scale,”
ChemElectroChem , vol. 1, 012014.[7] S. Menzel, P. Kaupmann, and R. Waser, “Understanding filamentarygrowth in electrochemical metallization memory cells using kineticmonte carlo simulations,”
Nanoscale , vol. 7, pp. 12 673–12 681, 2015.[8] S. Menzel, A. Siemon, A. Ascoli, and R. Tetzlaff, “Requirements andchallenges for modelling redox-based memristive devices,” in , 2018, pp. 1–5.[9] D. B. Strukov, G. S. Snider, D. R. Stewart, and R. S. Williams, “Themissing memristor found,” nature , vol. 453, no. 7191, p. 80, 2008. [10] M. D. Pickett, D. B. Strukov, J. L. Borghetti, J. J. Yang, G. S. Snider,D. R. Stewart, and R. S. Williams, “Switching dynamics in titaniumdioxide memristive devices,”
Journal of Applied Physics , vol. 106, no. 7,p. 074508, 2009.[11] R. Waser, R. Dittmann, G. Staikov, and K. Szot, “Redox-based resistiveswitching memories–nanoionic mechanisms, prospects, and challenges,”
Advanced materials , vol. 21, no. 25-26, pp. 2632–2663, 2009.[12] T. Chang, S.-H. Jo, K.-H. Kim, P. Sheridan, S. Gaba, and W. Lu,“Synaptic behaviors and modeling of a metal oxide memristive device,”
Applied physics A , vol. 102, no. 4, pp. 857–863, 2011.[13] Z. Jiang, Y. Wu, S. Yu, L. Yang, K. Song, Z. Karim, and H.-S. P. Wong,“A compact model for metal–oxide resistive random access memorywith experiment verification,”
IEEE Transactions on Electron Devices ,vol. 63, no. 5, pp. 1884–1892, 2016.[14] F. Corinto and A. Ascoli, “A boundary condition-based approach to themodeling of memristor nanostructures,”
IEEE Transactions on Circuitsand Systems I: Regular Papers , vol. 59, no. 11, pp. 2713–2726, 2012.[15] Y. N. Joglekar and S. J. Wolf, “The elusive memristor: properties ofbasic electrical circuits,”
European Journal of Physics , vol. 30, no. 4, p.661, 2009.[16] Z. Biolek, D. Biolek, and V. Biolkova, “Spice model of memristor withnonlinear dopant drift.”
Radioengineering , vol. 18, no. 2, 2009.[17] T. Prodromakis, B. P. Peh, C. Papavassiliou, and C. Toumazou, “A versa-tile memristor model with nonlinear dopant kinetics,”
IEEE transactionson electron devices , vol. 58, no. 9, pp. 3099–3105, 2011.[18] V. Ntinas, A. Ascoli, R. Tetzlaff, and G. C. Sirakoulis, “A completeanalytical solution for the on and off dynamic equations of a TaOmemristor,”
IEEE Transactions on Circuits and Systems II: ExpressBriefs , vol. 66, no. 4, pp. 682–686, 2019.[19] S. Kvatinsky, E. G. Friedman, A. Kolodny, and U. C. Weiser, “Team:Threshold adaptive memristor model,”
IEEE Transactions on Circuitsand Systems I: Regular Papers , vol. 60, no. 1, pp. 211–221, 2012.[20] S. Kvatinsky, M. Ramadan, E. G. Friedman, and A. Kolodny, “VTEAM:A general model for voltage-controlled memristors,”
IEEE Transactionson Circuits and Systems II: Express Briefs , vol. 62, no. 8, pp. 786–790,2015.[21] I. Vourkas and G. C. Sirakoulis, “A novel design and modeling paradigmfor memristor-based crossbar circuits,”
IEEE Transactions on Nanotech-nology , vol. 11, no. 6, pp. 1151–1159, 2012.[22] S. Datta, “Nanoscale device modeling: the greens function method,”
Superlattices and microstructures , vol. 28, no. 4, pp. 253–278, 2000.[23] S. Tappertzhofen, H. Mndelein, I. Valov, and R. Waser, “Nanoionictransport and electrochemical reactions in resistively switching silicondioxide,”
Nanoscale , vol. 4, pp. 3040–3043, 2012.[24] S. Menzel, S. Tappertzhofen, R. Waser, and I. Valov, “Switching kineticsof electrochemical metallization memory cells,”
Phys. Chem. Chem.Phys. , vol. 15, pp. 6945–6952, 2013.[25] I. Valov, “Interfacial interactions and their impact on redox-basedresistive switching memories (ReRAMs),”
Semiconductor Science andTechnology , vol. 32, no. 9, p. 093006, Aug 2017.[26] I. Valov and M. N. Kozicki, “Cation-based resistance change memory,”
Journal of Physics D: Applied Physics , vol. 46, no. 7, p. 074005, Feb2013.[27] A. Mehonic, A. L. Shluger, D. Gao, I. Valov, E. Miranda, D. Ielmini,A. Bricalli, E. Ambrosi, C. Li, J. J. Yang, Q. Xia, and A. J. Kenyon,“Silicon oxide (siox): A promising material for resistance switching?”
Advanced Materials , vol. 30, no. 43, p. 1801187, 2018.[28] S. Moysidis, I. Karafyllidis, and P. Dimitrakis, “Design and simulationof graphene majority gate without back-gating,”
Engineering ResearchExpress , vol. 1, 08 2019.[29] S. Moysidis, I. G. Karafyllidis, and P. Dimitrakis, “Graphene logicgates,”
IEEE Transactions on Nanotechnology , vol. 17, no. 4, pp. 852–859, 2018.[30] Y. Aharonov, L. Davidovich, and N. Zagury, “Quantum random walks,”