A Study of the Electroforming Process in 1T1R Memory Arrays

For reproducible resistive switching in memristive devices, electroforming is a crucial process. However, a deeper understanding of the electroforming process is still lacking due to unavailability of a proper simulation tool. Here, we propose a physics-based compact model for the electroforming of valence change mechanism (VCM) memristive devices. The developed JART VCM Forming model is experimentally validated with the ZrOx-based memristive device. Furthermore, the electroforming process in different 1T1R memristive arrays is simulated with this model. The study shows that the electrical characteristics of each device in the array after the forming process are influenced by word/bit line series resistance. In addition, control effects depending on the channel width and applied gate voltage of transistor in the 1T1R cell are also investigated with the compact model simulation.


I. INTRODUCTION
V ALENCE change mechanism (VCM)-based memristive device is one of the most attractive candidate for nextgeneration memory because of its energy-efficient, fast, scalable characteristics. The structure of the VCM device cell consists of a simple metal/oxide/metal stack, which makes it suitable for nm-scale electric device integration processes. In various type of oxides, such as HfO 2 , TiO 2 , Ta 2 O 5 , or SrTiO 3 , memristive switching phenomena were shown [1]- [3], with different combinations of metal electrode [4]. Typically, a pair of one reactive metal with a low work function, reactive electrode (also called ohmic electrode), and one inert metal with a high work function (also called electronically active electrode) are usually chosen. The different chemical reactivity of the electrodes leads to the partial reduction of switching oxide layer, which takes place dominantly at the interface with an (d) Corresponding forming voltage as a function of the oxide thickness compared to measurement data of a ZrO x -based VCM device (orange) VCM cell [10]. For the simulation of the thickness variation, v 0 = 2 · 10 12 Hz was used instead of the listed value in Table I. ohmic electrode due to an anodic oxidation [5]. This is the driving force for the electroforming process and leads to a nonuniform distribution of oxygen defects during the electroforming [6]- [9]. To understand the main kinetics of this process and the corresponding change in the electrical characteristics of the VCM cell, a physically well-defined simulation compact model is required. Fig. 1(a) shows the detailed mechanism of the electroforming and switching process in a filamentary VCM cell. The process is induced by the migration of oxygen vacancies within a filamentary region provided that an electric field across the VCM cell. The filament is divided into a plug region (nearby the surface between the ohmic electrode and the oxide) and a disc region (nearby the surface between active electrode and oxide). To trigger the electroforming process from the pristine state of the VCM cell, a positive potential is applied at the ohmic electrode (or a negative potential at the active electrode). Due to an anodic oxidation process, oxygen ions are extracted from the oxide and partially oxidize This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ the ohmic electrode [11]. The partially oxidized electrode will be still conducting. The anodic oxidation leaves behind two positively charged oxygen vacancies in the oxide, which will move toward the disc region with the applied electric field. Because the vacancies act as mobile donors, the resistance of the disc/plug region decreases with an increased oxygen vacancy concentration creating a conducting path. When a high current starts to flow through the filament region, the electroforming process is completed. The successive switching process is controlled by the polarity of applied voltage inducing a partial increase/decrease of vacancy concentration in the disc. The electrical characteristics of the VCM cell are highly influenced by the electroforming process. For example, the total number of oxygen vacancies in the filament is varied by the current compliance (I cc ) and overshoot during the forming process [6], [12].
While various type of compact models for the VCM switching mechanism has been proposed [12]- [15], compact models for the electroforming process have been rarely reported. We have recently proposed a compact model for the electroforming process that includes the concept of oxygen exchange between the plug region and the ohmic electrode. The equivalent circuit diagram of the compact model is shown in Fig. 1(b) [16].
For memristive-based neuromorphic computing applications, such as the convolutional neural network (CNN), an 1T1R array is required for vector matrix multiplication (VMM) [17], [18]. The VCM-based memristive devices offer multilevel programming capability, structural simplicity, and energy efficiency and have attracted attention for the kernel node device. For this application, the electrical characteristics after the forming process in the 1T1R array should be well understood for avoiding the poor device yield and nonuniformity. Furthermore, one must investigate various type of the array configuration, because the signal path of VMM-based learning process could be different with organization of the 1T1R array [19], [20].
In this work, a physics-based compact model for the electroforming process is applied to different 1T1R memristive crossbar array configuration for studying the electrical characteristics at each node. Especially, the influence of the parasitic line series resistance in different array configurations is investigated. Furthermore, the influence of the transistor parameters in the 1T1R cell is analyzed. The simulation results obtained with the VCM compact model in the 1T1R array give a guideline for the parameter design of memristive crossbar arrays for neuromorphic computing application.

II. ELECTROFORMING MODEL AND ITS VALIDATION
For modeling the electroforming process, a similar approach as in the Jülich Aachen resistive switching tool (JART) VCM models proposed in [13] is applied. The electroforming model was originally reported in [16], but the full equation system has been not well explained, yet. Thus, we here report a comprehensive model for the forming process. Each layer component in a memristive cell should be included in the equivalent circuit of the model. Fig. 1(b) shows the equivalent circuit diagram of the VCM cell. It consists of the resistance of the two electrodes, and the electrical components of the filamentary region as outlined below. The voltage is applied to the active top electrode while the bottom ohmic electrode is grounded.
In the proposed model, the filament is considered as a cylindrical region and divided into two regimes: 1) the disc region with the length l disc and resistance R disc and 2) the plug region with the length l plug and resistance R plug . The cross-sectional area of the filament is given by A = π r 2 fill and length l cell is equivalent to the oxide layer length. The interface between AE/disc and OE/plug denoted as V schottky,AE and V schottky,OE describes the electrical transport at each of the two interfaces. The electrical model is identical to the previous JART VCM v2 model [13]. In contrast to the JART VCM v2, oxygen exchange processes between the electrodes and the oxide region are also considered. The oxygen vacancy concentration in the disc N disc and the plug N plug is defined as state variables in the JART VCM forming model, but also the oxygen concentration in the OE N O,OE and the AE N O,AE are state variables. In the following, we first describe the electrical model and then the ionic model constituting together the JART VCM Forming model.
For calculating the main electrical characteristics of the VCM cell, the current flow through the cell is defined with Kirchhoff's law The resistances of the two regions of the filament are calculated with In (2) and (3), e is the elementary charge, z Vo is the charge number of the oxygen vacancies, k B denotes the Boltzmann constant, T is the temperature, μ n0 is the temperatureindependent prefactor of the mobility, and W is an activation energy for the temperature-dependent mobility.
The current flow through the Schottky diode is defined using thermionic and thermionic field emission in forward and reverse bias, respectively. In the VCM cell, the Schottky diode at the AE/oxide interface and the OE/oxide interface is connected in the opposite directions. The calculation of the current flow considers with a forward/reverse bias of applied voltage [21] At the AE/oxide interface, a positive voltage equals a forward bias (V forward = V Schottky,AE , V reverse = −V Schottky,AE ). In contrast, at the OE/oxide interface, a positive voltage corresponds to a reverse bias (V forward = −V Schottky,OE , V reverse = V Schottky,OE ).
The thermionic emission current I TE at the Schottky interface within a forward bias voltage is calculated by In (5), A denotes the cross-sectional area. The effective barrier height eφ Bn and the effective Richardson constant A * are calculated according to [21] Here, N D denotes donor concentration in the oxide, in VCM cell oxygen vacancies take a role of double-charged donor, and it is applied within N D = 2N Vo . The oxygen vacancy concentration N Vo is defined as the disc vacancy concentration N disc at the AE interface and the plug vacancy concentration N plug at the OE interface. eφ Bn0 is the energy difference between the metal work function and oxide electron affinity, m eff is the effective mass of the carrier, and φ B is the permittivity related to the barrier lowering.
Under reverse bias voltage, the thermionic field emission current I TFE is calculated by [21] The parameters W 00 , W 0 , and ζ are described as The energy difference between the Fermi level and the conduction band edge on each electrode is denoted with φ n,AE/OE , calculated based on the Fermi-Dirac integral by Nilsson [22] For implementing the electroforming process, the extracted oxygen ions, which is stored in the OE, is also considered. Thus, three main state variables, N disc , N plug , and the oxygen ion concentration in the OE N O,OE , are used for simulating this process. Three coupled ordinary differential equations (ODEs) are calculated in the model for these three state variables. The change of N disc and N plug is modeled with three types of ionic currents: I ion denotes the flux of oxygen vacancies between the plug and the disc. The oxygen exchange between disc and active electrode and the oxygen exchange between plug and ohmic electrode are described by the ionic currents I ion,BV,AE and I ion,BV,OE . The corresponding ODEs are read The ionic current I ion between the plug and the disc is represented with the drift current I ion,drift and the diffusion current I ion,diffusion and corresponds to [13] In (16)-(21), a is the hopping distance, ν 0 is the attempt frequency, W A is the migration barrier, and N min and N max are the limits of the vacancy state variables. The minimum vacancy concentration N min is defined as one defect in the respective volume and calculated as listed in Table I. During electroforming, the oxygen exchange process happens at the OE due to anodic oxidation. The oxygen exchange reaction at the AE can be neglected in this case. The influence of I ion,BV,AE element is also neglected in the ODEs. The change of the oxygen concentration in the OE is determined by I ion,BV,OE , and its ODEs are given by In (22), l OE,eff denotes the active thickness of the OE in which the extracted oxygen is stored.
To define the ionic exchange current density at the interface of an electrode, the exchange reaction is described with The ionic exchange current density by oxidation/reduction of the oxide is modeled using a Butler-Volmer equation, J O,I denotes the current density for forward reaction, and J O,II for reverse direction as k 0 I (k 0 II ) is the reaction rate coefficient for the reduction (oxidation) of the oxide, α is the charge transfer coefficient, z O is the charge of an oxygen ion, and η is the overpotential in the Helmholtz layer of the ion transfer system. N O,ad(M) is the concentration of oxygen ions absorbed in the metal, and N O,ox is the concentration of oxygen on lattice sites in the oxide. The corresponding activation energy for reduction (oxidation) of the oxide is denoted as G I ( G II ) [23]. At the surface of each filament region, the total ionic current density is defined by the oxygen exchange from metal electrode (M) to the oxide (Ox), and J O,M→Ox is calculated with the sum of the oxidation and reduction ionic exchange current densities The oxygen exchange current at the interface between the OE and the oxide I ion,BV,OE in (15) is calculated by The oxygen concentration in the plug N O,plug results from the total oxygen concentration N O,oxide,max of the oxide material and the oxygen vacancy concentration in the plug to The limiting factors F limit,OE,I(II) are window functions constraining the state variables to their limits. They are defined according to Here, N O,OE,min is the minimum oxygen concentration in the OE and depends on the effective volume of the OE in which the oxygen is "stored" as listed in Table I. The minimum concentration N mO,OE,min is equal to the case that exactly one O atom in this volume.
The electric field E in (16) and (20) is defined according to The local temperature T is calculated based on the power dissipation in the filament due to the Joule heating where I is the electronic current, R th,eff is the effective thermal resistance, and T 0 is the initial temperature of ambient. The simulation model is implemented in MATLAB and Verilog-A, which runs within the CADENCE environment. Table I lists all simulations parameters used in this study.
To validate the electroforming compact model, simulations with parameters based on Pt/ZrO x /Ta/Pt device are performed and compared to the experimental data published by Kindsmüller et al. [10]. In this study, a dependence of the electroforming voltage on the thickness of the ZrO x layer is observed. It shows that the forming voltage increases linearly with the oxide thickness. The Pt/ZrO x interface is regarded as the AE/oxide interface, and the Ta/ZrO x interface is the OE/oxide interface. For the electroforming simulation, the same layer thicknesses as in the experiment are used. Likewise, a voltage sweep rate is also applied with 1 V/s in the simulation. Note that we chose another polarity convention than in the experimental study setting the Ta/Pt electrode to GND and applying the negative voltage to the active Pt electrode. The forming voltage is defined as the voltage point at which a current level of 1 μA is reached. This criterion is the same in experiment and simulation. Fig. 1(c) shows the simulated I-V characteristics (blue line) of the electroforming processes using a 5 nm-thick ZrO x layer in comparison to the experimental data (gray line). The simulation result shows a good agreement with the measurement data at voltages |V applied | > |−1 V|. At the other range of voltage (|V applied < | − 1 V|), the measurement data are within the resolution limit of the measurement setup leading to a plateaulike feature. At high voltages, the forming event appears as an abrupt increase of the current value.
We also studied the thickness dependence of the forming process in the simulation. The same range of the ZrO x layer thicknesses was used as in experiment [10]. The simulation results feature that the forming voltage |V Forming | is increased with oxide thickness. The simulated trend of the forming voltage plotted in Fig. 1(d) also has a good agreement with the experimental measurement data.

III. ELECTROFORMING SIMULATION IN 1T1R CELL
In the memristive 1T1R memory array, each node cell consists of one transistor and one memristor connected serially. To simulate the electroforming in the 1T1R single cell, the JART VCM Forming model and a BSIM 4 transistor model for bulk CMOS are applied. Fig. 2(a) shows an illustration of the simulated 1T1R single cell. The voltage source for electroforming V WL is connected with the OE of memristor via the word line resistance R WL . The AE of memristor is directly connected to the drain of transistor. The source of transistor is connected to GND via the bit line resistance R BL . The current flows through the word line resistance, the memristor, and the transistor channel over the bit line to GND. In addition, the voltage source for the transistor control V select is connected to the gate of transistor. The memristive devices are arranged in a way that they form with a negative voltage applied to the word line to be consistent with the voltage convention in the previous section. In reality, the memristive device in the 1T1R cell would be placed "upside down." Then, a positive voltage needs to be applied to the OE in order to form the cell.
In this simulation, we have chosen the 45-nm CMOS technology node, where the transistor width and length are kept to be same as the node. In the memristor model, the maximum vacancy concentration of N max listed in Table I is replaced by 8 · 10 27 m −3 . The transistor model includes a parasitic resistance factor, i.e., the junction resistance at the drain [24], [25]. Its value is smaller than the series resistance parameter R series , which can be considered as a contact resistance in the VCM model. The influence of the junction resistance between VCM cell and transistor is neglected, but is supposed to be lower than the segment resistance. Fig. 2(b) and (c) shows the simulated electrical characteristics of the memristor in a single 1T1R cell. In this simulation, each line series resistance for WL/BL was assumed to be 10 and voltage sweep rate was 1 V/s for both memristor OE and transistor gate electrode. The channel current was limited by a gate voltage of 0.5 V. During the electroforming phase when the forming voltage 1.5 V was applied for 6.5 s to the VCM cell, the transient current level suddenly increased after 5 s. The vacancy concentration within each layer of the filament also increased at the same time as shown in Fig. 2(c). While the disc concentration is almost constant at first, only the plug concentration increases. Then, the disc concentration starts increasing, but slower than plug concentration. After about 8 s, both concentrations increase abruptly also leading to the abrupt current increase in the I-V characteristic. This sudden increase is due to the positive feedback of Joule heating effect because of a thermal runaway similar to the case reported for the SET resistive switching transition from a high resistive to the low resistive status [26], [27]. With this simulation result, the electroforming process in the 1T1R single node cell is verified with the compact model simulation

IV. MEMRISTIVE 1T1R ARRAY CONFIGURATIONS
Three different types of 1T1R crossbar array with m word lines (WL) and n bit lines (BL) are feasible when the lines are only structured horizontally or vertically. Each array unit consists of a memristive VCM cells (1R), and a select transistor (1T) to prevent unintended sneak path currents and half-select issues [28]. The select lines (SLs) are connected to the transistor gates for activating the transistor of each 1T1R cell. Fig. 3 shows each type of 1T1R array configuration. In case of the typical-type array in Fig. 3(a), the voltage for the electroforming of VCM devices is applied to the vertical WL, and the current flows through the WL and the unit cell over the BL to GND, when the transistor channel at the targeted cell is conducting. The SL is connected to the gate of the transistors along horizontal lines. The operation of each transistor is controlled by the voltage on the SL. In contrast to the typical array, the vertical-type array in Fig. 3(b) has vertical SLs for each transistor gate, and the voltage for the electroforming is applied to the horizontal WLs. In case of the pseudocrossbar array in Fig. 3(c), both the WL and the SL are aligned in parallel in horizontal direction. The devices on the same WL can be operated at once. In this way, the cell operation on the line cannot be controlled individually. In the neuromorphic application system, an analog VMM with the crossbar array is presented by M ij is the synaptic weight of each array node, w i is the input read signal for mapping, and z i is the accumulative decipher signal with output.
In the 1T1R memristive device array, the analog conductance value of the memristor in each cell G ij reflects the synaptic weight, the read voltage applied to the WL v i is the input read signal, and the output current sensed at the BL I j is the output signal. The relation with these elements is calculated as the VMM Since the 1T1R array has different types of configurations, the most suited type of array should be chosen depending on the required array operations and applications. In case of the typical-crossbar array, the transistor in each cell can be opened for several clock cycles to achieve analog inputs, and the BL can be charged to different levels. This scheme requires a capacitive read out (voltage sensing). For a correct operation, the charging time of BL should be longer than the longest time one of the transistors remains open. Thus, it can make the computing time longer. In case of the vertical array, different analog voltages can be applied to the WL and the accumulated current can be read at the BL. The programming analog values are best performed on a single-cell level. However, it comes at the expense of an increased energy consumption to prevent write disturb. Otherwise, in case of the pseudocrossbar array, a full vector-matrix product can be performed in one signal cycle by reading the information on the BLs in parallel. During this operation, however, this array is a passive array, and sneak paths can make read failures. In addition, it is not clear if all BLs can be read in parallel anyway as the required read circuit might become too complex and large. Furthermore, the complex read circuit for the pseudocrossbar array can make the integration of neuromorphic system inefficient [17], [29].
In this research, the line series resistance of WLs and BLs is considered in the various array configuration. Each line segment resistance R seg is assumed to be identical for the BL and the WL, i.e., R BL,seg = R WL,seg = R seg . Moreover, a quadratic configuration of the array with n = m is considered. The total series resistance R ser on each cell can be different according to array type. The critical path as the longest way of current flow is colored as an orange line, and the initial path as the shortest way of current flow is colored as a sky blue line is drawn in Fig. 3(a)-(c). The R ser for the initial/critical path of each arrays are also calculated.
In the typical array, each cell on the same column has a uniform R ser , because it is calculated as R ser = (m − i) · R WL,seg + i · R BL,seg when the cell (i, j) on the ith row and the jth column is formed. R ser of critical path is R ser = m · R seg , same as the case of the initial path. Therefore, there is no change on R ser , and a minor effect of electroforming on the position of the cell is expected. Note that using different materials (and/or segment length) for WL and BL, the total line resistance will differ slightly, but the length of the critical path from input voltage to GND remains the same.
In case of the vertical array and pseudocrossbar array, the series resistance of each cell is calculated with R ser = j · R WL,seg + i · R BL,seg . R ser can change dramatically depending on position (i, j) in the array on each cell. The maximum R ser is found with the critical path at position (m, n) as R ser = (m + n) · R seg , and the minimum R ser is found with the initial path at position (1, 1) as R ser = 2 · R seg . With this difference of R ser , we can expect a larger effect on the characteristics after forming in the vertical and pseudocrossbar array than typical-type array depending on the array cell position. Fig. 3(d) and (e) shows an example of the simulated electrical characteristics distribution after the forming of each VCM cell in the vertical-type 1T1R array circuit. At point (1, 1) of the array, the cell has the highest value of current level after forming, because it has lowest series resistance. Especially, the vacancy concentrations after forming is also tilted with the line series resistance at each point. The results verify that the forming phenomena during the process are highly influenced by the line series resistance.
In the following, we translate the 3-D plots to 2-D box plots for better comparability. Each box in Figs

V. ELECTROFORMING SIMULATION OF 1T1R ARRAY
A 1024 × 1024 1T1R memristive cell array is considered in this simulation. To form the cell (i, j), in the typical array case, V select,i is applied on the ith SL to open the transistor channel, the forming voltage is applied to V WL,j on jth WL, and the BL is set to ground state GND. On the contrary, in the vertical array, V select,j is applied on the jth SL to activate transistor, and the forming voltage is applied to V WL,i on the ith WL. In case of the pseudocrossbar array, V select,i is applied on ith SL and  the forming voltage is applied to V WL,i on the ith WL. The signals on the other lines are set to 0 V. The voltage pulse for electroforming consists of the forming phase and the read phase, both have a rectangular shape. In the forming phase, the voltage potential applied to the WL increases to 1.5 V with a sweep rate of 1 V/s. Then, the voltage is kept constant for 6.5 s and declines again to zero. During this phase, the voltage of 0.5/1 V is applied to SL and the transistor gates. Two different voltage levels are used to investigate the impact on the programmed states after electroforming. In the read phase, the signal of 0.3 V is applied to WL to read the current level. In addition, the vacancy concentration in the disc after forming at each cell is evaluated in this phase. In contrast to the forming phase, the voltage of 1.5 V is applied on the SL to fully open the transistor channel in the cell during the read phase. In this simulation, 45/90 nm was chosen for the transistor channel width condition, and a constant length 45 nm was applied. In addition, maximum vacancy concentration N max in the JART VCM forming model was kept 8·10 27 m −3 , same as for the 1T1R single cell simulation. Fig. 4 shows the simulated electrical characteristics of each cell in a typical 1T1R array. The implemented R seg was chosen to 0.1, 0.5, 1, 10, and 50 for investigating the impact of line resistance on cell characteristics after the forming. The quantitative maximum/minimum value of simulated results on each condition is represented in Table II. Fig. 4 shows that the current level and the vacancy concentration in the disc N disc after the forming strongly decreases for higher R seg , and for lower V g . Especially, the meaning that N disc is changed with R seg and V g indicates that the forming process on each cell was significantly influenced by R ser . The distribution of the characteristic values is also getting larger with higher R seg . When the transistor width is 45 nm, the N disc value after forming is limited even though N max is set to 8 · 10 27 m −3 . This limit can be released with increasing transistor width to 90 nm, and N disc increases and approaches N max for low R seg and higher V g range.  Fig. 5 shows the characteristics in the vertical 1T1R array. Similar trends of characteristics with R seg and V g dependence are also shown in the vertical 1T1R array simulation. The quantitative maximum/minimum value of simulated results are represented in Table III. In contrast to the case for the typical 1T1R array, the distribution range of the characteristics is much larger, and also widens for the higher R seg . The vertical 1T1R array has a larger range of R ser depending on the array point, since the length of the critical paths varies. Thus, the characteristics of each cell after the forming are more influenced by R ser in the vertical-type array. N disc after the forming also remains below N max for 45-nm transistor channel width in the vertical 1T1R array. This, however, is inflated to N max for 90-nm transistor channel width. Fig. 6 shows the electrical characteristics in the pseudocrossbar array. The quantitative maximum/minimum value of simulated results is represented in Table IV. Also, in this case, the current level after electroforming and the vacancy concentration in disc decrease for higher R seg , and for lower V g . The trends of the characteristic values and the phenomena with limited N disc for 45-nm transistor channel width are also similar as for the vertical-type array. In addition, the distribution range of the characteristic parameters is much larger than the typical-and the vertical-type array. In integrated memory cells, BL and WL may consist of different materials. Moreover, the pitch in lateral and horizontal direction may be different as BL and WL could be integrated in different metal layers of the integrated stack. Thus, in general, different segment resistances for BL and WL are expected. To further study the influence of the line series resistance in the arrays, we performed simulations with an asymmetric line segment resistance for WL and BL. In the simulation, R seg,WL is fixed to 1 , and R seg,BL was chosen to 0.1 (i.e., R seg,WL > R seg,BL ) and 50 (i.e., R seg,WL < R seg,BL ). The transistor channel was 45 nm, and V g = 0.5 V was chosen for the simulation condition. Fig. 7(a) and (b) shows the simulated read current level after forming at each array point in the typical-type array. In Fig. 7(a) with R seg,BL = 0.1 , the influence by the line series resistance for WL should be dominant. The current level after forming is getting higher with the increasing point number of row i, because the series resistance in the typical-type array is identified with R ser = (m−i)·R WL,seg +i·R BL,seg . In contrast, in Fig. 7(b) with R seg,BL = 50 , the current level after forming is getting lower with the increasing point number of row i. Fig. 7(c) and (d) shows the simulated read current level after forming in the vertical-type array. In this case, the series resistance in the vertical-type array is identified with R ser = j · R WL,seg + i · R BL,seg . In Fig. 7(c), with R seg,BL = 0.1 , the influence by WL resistance is dominant, and the current level after forming is getting lower with the increasing point number of column j. In contrast, in Fig. 7(d) with  Fig. 7(e) and (f). The series resistance in the pseudocrossbar array is similar to the verticaltype array, and the influence of the line series resistance is comparable to the case of the vertical-type array.

VI. DISCUSSION
As the array size is getting larger for neuromorphic application, the impact of the line series resistance on each array cell becomes more severe [30]- [32]. The best way to minimize this impact is that the larger height and wider width of lines are used for array design. This might however not be an acceptable solution due to the scaling issue. Another approach is the control of the final state by the applied gate voltage or implemented transistor channel width.
Each type of array offers it advantages and disadvantages. As the critical path and the initial path are identical for the typical 1T1R array, it provides a better control on the forming state. This can be observed when comparing the scatter in the box plots. Even when the series resistance is clearly influencing the final characteristic, the variation over the array is smaller than in the other two array types. The smaller variation in the formed state potentially reduces the device-to-device variability across the array during switching [30], [33]. For a pure memory operation, this type of array is certainly a better choice, however, this suffers from the drawback discussed above during the VMM operations. An analog VMM can only be achieved by time coding of the input vector and measuring the charge over time. In this case, a constant read voltage is used, but the transistors are opened/closed for a certain number of cycles corresponding to the analog input value. This will cost more time. If the memory cell, however, shows a nonlinearity in the I-V characteristic using this scheme might be anyway better suited for the accuracy of the VMM operation. Otherwise, the currents would not scale linearly anymore with the analog read voltage. Additional issue with using analog voltage inputs is the read disturb. As the switching kinetics are highly nonlinear, the VCM cell will switch at any voltage, but at orders of magnitude different time scale. This means increasing the read voltage may lead to undesirable changes in the stored synaptic weight.
In the present study, we did not consider the inherent variability of the VCM cells [34], as the focus of this study was on the understanding of the interplay between the VCM cell and the array. Device-to-device variability of the VCM cells, however, is an important issue that needs to be considered. There are three ways of including the variability in the compact model. The first approach takes the similar way as in this article from Bengel et al. [35]. In this case, the parameters of the model are taken randomly from a truncated Gaussian distribution. Alternatively, one can consider corner cases, meaning the fastest and slowest forming device. As shown in [35], this approach already gives the expected range of switching variability. In another work from Wiefels et al. [36], it was shown how to include the read instability in the JART VCM compact model. In principle, the forming model could be extended in the same manner as in these papers. A third approach, including the variability in the model during the forming process, would be to include a noise source as additional term for the change of the state variables or the temperature as shown by Guan et al. [37]. The variability issue will be also included in the future version of the JART VCM forming compact model.

VII. CONCLUSION
In this research, a highly predictable electroforming compact model for the VCM-based memristive device was developed and applied for the study of 1T1R array application for various configuration. The compact model was fitted with the ZrO x -based VCM device parameters, and validated with its experimental measurement data. First, the detailed electroforming phenomena of the model were checked with the simulated electrical characteristics in the 1T1R single node cell.
With the simulation result in 1T1R array configurations, it was confirmed that the line series resistance for each array point takes a critical role to impact the cell characteristics after the electroforming. In case of a node point, which has the shortest path from voltage source to GND, the highest value of current level and vacancy concentration after the forming is observed. The distribution of the characteristic values scales with the value of line segment resistance. In addition, the impact of line series resistance can be mitigated with the different gate voltage or channel width condition applied on the transistor part of the 1T1R cell. This simulation study gives a guideline for design parameters in various type of the 1T1R array configuration. From the view of the variability of the formed state, our study shows that the typical 1T1R array structure is the best choice as the length of the current paths is always the same, meaning the series resistance due to wire resistance is always identical.