Transfer Modeling for 1T1R Crossbar Arrays with Line Resistances

Progress of neuromorphic computing and next-generation information storage technologies hinges on the development of emerging nonvolatile memory (eNVM) devices, which are typically organized employing the crossbar array architecture. To facilitate quantitative performance analysis of eNVM crossbar array architecture, this paper proposes a way to study the one-transistor-one-resistor (1T1R, R: eNVM devices) crossbar arrays based on matrix algebra method. The comparative analysis of 1T1R crossbar array modeling based on matrix algebra method and compact-model SPICE simulations verifies the accuracy of the proposed method, which can be directly used for static quantitative analysis and evaluation of 1T1R crossbar array performance. With the proposed method, the optimization of array operation schemes and current backflow issue are discussed. Our analysis indicates that the proposed method is capable of flexibly adjusting array parameters and consider the influence of line resistance on array operation, and can provide guidance for improving the sensing margin of the array through multi-parameter co-simulation. The proposed matrix algebra-based 1T1R crossbar array modeling method can bridge the gap between the accuracy and flexibility of the available methods.


Introduction
The rapid development of artificial intelligence (AI) technology represented by deep learning has made it possible for human beings to enter a highly intelligent era [1][2][3][4]. This technology requires high-bandwidth data storage and computing capabilities, and places high demands on hardware [5,6]. From a hardware perspective, the crossbar array of nonvolatile memory (NVM) devices has received extensive attention in recent years [7][8][9]. In a crossbar array, NVM devices can be integrated at the junctions of orthogonal access lines [9][10][11]. Due to the advance of emerging NVM (eNVM) devices such as phase-change random-access memory (PCRAM) [12,13], resistive RAM (RRAM) [14,15], etc. and their compatibility with CMOS technology, one-diode-one-resistor and one-transistor-one-resistor crossbar arrays (each junction of crossbar arrays integrating an eNVM device and a diode or transistor, abbreviated as 1D1R and 1T1R, respectively, R: eNVM devices) are promising architectural routes to realize future information storage or neuromorphic computing [8,[16][17][18][19].
Whether it is applied to information storage or neuromorphic computing scenarios, quantitative electrical analysis for 1D1R or 1T1R crossbar arrays is a crucial task, and many efforts have been made [20][21][22][23][24][25]. Most of these researches are based on SPICE-tool due to its highly optimized and compact features (high accuracy). However, in the case of a large-scale crossbar array, its inherent node analysis characteristics result in a slow calculation speed, which poses high demands on the computation capability. Moreover, SPICE-based circuit simulation lacks flexibility in adjusting several array parameters, and usually does not consider the influence of line resistance on array operation [7,23]. Therefore, two methods of crossbar array modeling based on numerical iterative are proposed, which are implemented adopting C and MATLAB codes respectively [23][24][25]. These two methods can easily adjust the array parameters and take into account the influence of line resistance on array operations. However, in available literatures that based on these two methods, only 1D1R crossbar arrays are considered, and the accuracy of these array modeling is not evaluated. Among these two numerical modeling methods, the matrix algebra-based method is easy to implement parallel computation and sparse matrix storage, which is more suitable for modeling crossbar arrays [23,25]. But so far, the matrix algebraic modeling and quantitative analysis for 1T1R crossbar arrays have not been reported.
Under the theoretical framework of [24], this paper implements the 1T1R crossbar array modeling based on the matrix algebra method on the basis of [25] through transfer modeling. In this work, we propose a method of dynamically fitting each 1T1R cell to improve the accuracy, and complete the array calculation by adopting a piecewise function call method according to the location of each array cell. Line resistances are considered in the modeling, and the matrix algebra calculation is implemented as per Kirchhoff and Ohm's law. The accuracy of the 1T1R crossbar array modeling in this paper is verified by comparing with SPICE-based circuit simulation. The matrix algebra-based 1T1R crossbar array modeling method proposed in this paper can give consideration to flexibility and accuracy, so as to provide flexible and reliable quantitative analysis guidance for research groups based on 1T1R crossbar array architecture.
To realize the transfer modeling from 1D1R crossbar array to 1T1R crossbar array, the main challenges are: (i) Compared with the two-terminal unipolar passive diode, the MOSFET is a three-terminal active device, which increases the difficulty of array configuration, and a number of its parameters are sensitive to voltage; (ii) The 1D1R cell can be represented by a mature model function (can be expressed by Lambert-W function [26]), while the 1T1R cell has a complex multi-region working mode, and there is no ready model function for reference. In response to the two points, we give a detailed modeling solution from 1T1R cell to array in the subsequent paper. Section II introduces the modeling and dynamic fitting method of 1T1R cell. Section III gives the modeling method and details of 1T1R crossbar arrays. The simulation results and analyses of a typical 1T1R crossbar array are included in Section IV. Section V presents a brief conclusion and discussion.  The 1T1R crossbar array is a popular array architecture, which is usually emplyed to implement memory and computing integration, with the R stands for eNVM devices [27][28][29]. Generally, the eNVM devices mainly include PCRAM and RRAM, and so on [9]. Figure 1 shows the schematic diagram of a typical 1T1R array architecture and the internal connections of a single 1T1R cell, in which the transistor is enhanced N-channel metal-oxide-semiconductor field-effect transistor (NMOSFET), and its drain (D) is connected to the bottom electrode (BE) of the eNVM device. The gate (G), source (S), and substrate (B) of the NMOSFET are connected to the word line (WL), source line (SL), and the lowest potential of circuit (generally grounding), respectively. The top electrode (TE) of the eNVM device is connected to the bit line (BL). In the array architecture of Fig. 1, the main role of each device is vividly illustrated by replacing the 1T1R cell with a switch in series with a variable resistor.
Compared with the diode selector in 1D1R array, the NMOSFET in 1T1R array possesses three terminals, which has better control initiative, but increases the difficulty of modeling. In 1T1R crossbar arrays, the transistor plays the role of a switch under ideal circumstances (its resistance is 0 when it is turned on, and its resistance is infinite when it is turned off), but in fact it cannot be simply equivalent to a switch for its many second-order effects directly or indirectly affect array's performance. To make matters worse, the impact of these second-order effects on large-scale arrays will be more pronounced. Therefore, to model a 1T1R crossbar array, the single 1T1R cell should be analyzed first.
The cell schematic modelled in this work is a controller NMOSFET in series with a binary resistor representing high-or low-resistance state (HRS, LRS), as shown in Fig. 2. For the transistor used as a controller in 1T1R crossbar arrays, the working region of it mainly includes subthreshold (for the convenience of analysis, the cutoff region is merged into the subthreshold region) and triode regions. As per the series connection mode in Fig. 2, one can get the drain current of the NMOSFET working in the subthreshold and triode region as where I0 is a parameter related to the process, VGS and VTH are respective gate-source voltage and threshold voltage, ζ is a non-ideal factor, VT is thermal voltage (typically 26mV), μn, Cox, W and L are respectively electronic mobility, oxide capacitance per unit area, and channel width and length, VBS is the potential difference between BL and SL, R is the resistance of the eNVM device. In (2), K=μnCox is called transconductance parameter, VDS is the drain-source voltage of NMOSFET and can be easily deduced equal to (VBS-IDR) according to Fig. 2. It should be noted that the subscript B in VBS does not refer to the substrate of NMOSFET. A fact is that the deviation appears when the drain current of the 1T1R cell computed by using the equations (1)(2) in MATLAB and simulated by using the mature Silterra 0.18μm process in Cadence under the same set of variables {VGS, VDS}. This is because several parameters in (1) and (2) are treated as constants in the calculation process, but the actual situation is that these parameters themselves are affected by VGS and/or VDS. For example, the non-ideal factor ζ is affected by VGS, K=μnCox is affected by VDS and VGS, and VTH is affected by VS (substrate bias effect, VS represents the source potential). When modeling without considering these effects, a single 1T1R cell has a large deviation, and it is foreseeable that larger errors will occur in larger arrays. Therefore, in order to make the modeling more accurate, it is very necessary for us to consider these influencing factors. For this reason, we propose a method of directly fitting several critical parameters of the equations by comparing to the SPICE-based circuit simulation results (in Cadence), and then embedding their corresponding equations. For NMOSFET operating in the subthreshold region (see (1)), the key to the subthreshold equation of NMOSFETs lies in the values of I0, VTH and ζ, which can be obtained in the corresponding process library file or calculated in Cadence by using specific process (in this work, Silterra 0.18μm process, transistor nhp). The value of VTH is 553mV (in subthreshold, the impact of VGS on VTH is ignored), and I0=29.4116μA is obtained when VGS=VTH. Non-ideal factor ζ is susceptible to changes in VGS. To establish the function ζ(VGS), one can perform parameter sweep analysis on VGS to obtain corresponding ID, and then get ζ through ζ=(VGS-VTH)/[VTln(I0/ID)]. The data of ζ and VGS calculated through this equation are fitted by using the nonlinear least square method as implemented in MATLAB, and Fig. 3 shows the data points and the fitted function curve. These two variables conform to an exponential relationship, and the ζ(VGS) function obtained by fitting is where a, b, c are fit parameters. Thus the subthreshold formula (1) can be rewritten as Later, we use (4) to represent the unselected cells of 1T1R array working in cutoff and subthreshold regions. For NMOSFETs operating in the triode region (see (2)), K=μnCox is affected by VDS and VGS, and VTH is affected by VGS. In Cadence, parameter scanning analysis is used for VGS and VDS, and it is found that the change of K with VGS is much smaller than that of VDS. In addition, in a typical 1T1R array, WLs (connect to gate of NMOSFETs) are used as control terminals, and their voltage fluctuation are small, while BLs (connect to the TE of eNVM devices) are used as read and write operation terminals, and their voltage fluctuation are relatively large. In other words, the fluctuation of VDS (also susceptible to changes in the resistance of the respective series-connected eNVM devices) is much greater than that of VGS. Therefore, in the fitting process, only the function relationship between K and VDS is established, and the slight influence of VGS is ignored. To establish the function K(VDS), one can perform parameter sweep analysis on VDS to obain corresponding ID, and then get K through where VTH, W and L are set as fixed values (VTH=553mV, W=4μm and L=0.18μm) to derive the relationship between K and VDS. The data of K and VDS calculated through this equation are fitted by using the linear least square method as implemented in MATLAB, and Fig. 4(a) shows the data points and the fitted function curve. These two variables conform to an linear relationship, and the K(VDS) function obtained by fitting is 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 where p and q are fit parameters. Similarly, to establish the function VTH(VS), one can perform parameter sweep analysis on VS, and then get the fit result (linear least square method) where z and v are fit parameters. Fig. 4(b) shows the data points and the fitted function curve. Through the previous fitting, the triode formula (2) can be rewritten as It is worth noting that the abscissa ranges in Figs. 3-4 are determined according to the working range of the corresponding variables in the 1T1R crossbar array. Table 1 lists the fit functions and corresponding parameter values for the transistors in 1T1R cells working in subthreshold and triode regions.  In order to evaluate the effect of the above fitting, we simulated the I-V characteristic curves of a single 1T1R cell in the subthreshold and triode regions, including unfit and fit results in MATLAB and circuit simulation results in Cadence, as shown in Fig. 5. Fig. 5(a) shows the ID-VGS simulation curves of the transistor working in subthreshold region, and Fig. 5(b) shows the ID-VDS simulation curves of the transistor working in triode region. However, the influence of source potential VS on threshold voltage VTH is difficult to be obtained by I-V characteristic curve simulation. Nonetheless, Fig. 5(a)(b) suffice to illustrate that by fitting several key fluctuation parameters of a single 1T1R cell (the fitted MATLAB simulation results are very close to the Cadence simulation results), the modeling accuracy of the 1T1R cell can be effectively improved.

Modeling for a typical 1T1R crossbar array
A typical connection mode of m×n 1T1R crossbar array is shown in Fig. 6, including m WLs and SLs, and n BLs. Different from 1D1R crossbar array, in 1T1R crossbar array, the WL is connected to the gate of the control transistor to achieve on and off control, the SL collects current, and the BL is used for inputting read or write voltages.
It is well known that the control devices of 1T1R crossbar arrays are MOSFETs, whose gates are insulated gate electrodes and have very high resistance (approximately 10 8~1 0 9 Ω). Compared with the gate resistance, the line resistances along the WL direction are very small and almost negligible. Under normal WL voltage input, almost no current flows through the entire WL line, and there is almost no voltage drop on the line resistances (verification by Cadence simulation, e.g. for transistor nhp of Silterra 0.18μm process, when the voltage applied to the WL of selected cell is 5V and if the line resistance between every two adjacent junctions along with WL are 5Ω, the total voltage drop across WL is only 3μV for 100×100 array scale with NMOSFET parameters W/L=4μm/0.18μm). So that we can get the schematic diagram of simplified 1T1R crossbar array that shown in Fig. 7, in which the voltages on each junction and the currents on the connecting line are marked. According to Kirchhoff's current law (KCL), two current continuity equations can be obtained as ( , ) ( , ) ( , 1).

Define input parameters
The steady-state electrical properties of a crossbar array can be totally described by a series of voltage variables at every crosspoint of the array. There are three voltages at each crosspoint in a crossbar array, {VWL(i, j)}, {VBL(i, j)}, and {VSL(i, j)}, with 1≤i≤m and 1≤j≤n, on the WL, BL, and SL plane, respectively. In addition to these 3mn voltage variables, some input parameters must be defined to calculate the entire array, including applied voltages (Vapp_WL, Vapp_BL, and Vapp_SL), terminal resistances (Rsens_WL, Rsens_BL, and Rsens_SL) and line resistances (RWL, RBL, and RSL). The corresponding array parameters are shown in Fig. 8, where the line resistance is defined as the resistance between adjacent two corsspoint along WLs, BLs, and SLs. Voltages can be applied from one side or both sides of WLs/BLs/SLs. Terminal resistance is defined here as the resistance between WLs/BLs/SLs and power or ground. In order to make the calculation converge to be able to simulate the floating situation, the terminal resistances between floating lines and voltage sources can be set to a ultra high value (e.g. 100MΩ). For reading operations, terminal resistance of the selected SL needs to be set to a sensing resistance (e.g. 100Ω). Fig. 8 Schematic diagram of the defined input parameters for 1T1R crossbar array.

Matrix algebraic method
According to the Kirchhoff's voltage law (KVL) and Ohm's law, 6 types of equations as shown in (10) can be obtained on the basis of (9), and the additional four types of equations relate to cells at the edge of the array where their voltages correspond to the applied voltages at both sides of the BLs (Vapp_BL1 and Vapp_BL2) and SLs (In 1T1R crossbar array, WLs are used for gating operation, and SLs are only connected to sensing resistors, which are equivalent to Vapp_SL1 and Vapp_SL2 grounding.). The voltage and resistance variables and line resistance in (10) correspond to those in Figs. 7 and 8. These 6 types of equations include 2mn equations after expansion, which can be used to solve 2mn unknown junction voltage variables. The resistance of each cell is given by {R(i, j), 1≤i≤m, 1≤j≤n}, where R(i, j)∈{Rlow, Rhigh}. A detailed matrix-based array solution method is given in the appendix, in which the equations have modified for the reason of the different array architecture between this work and [24]. These equations can be written in matrix form using MATLAB to calculate and solve.
The code in [25] enables the "parfor-loop" and uses "sparse" method for matrix storage to improve the calculation efficiency. This paper carries out the calculation and simulation of the 1T1R crossbar array based on the code (with many modifications, transfer modeling), and combines with SPICE-based circuit simulation to verify its effectiveness. Lambert-W function is used in 1D1R crossbar array modeling, piecewise function method is used in our 1T1R crossbar array modeling (Equations (4) and (8) for subthreshold and triode regions, respectively.).
In our code, the gate-source voltage and cell voltage of each cell are initially as VGS(i, j)=Vapp_WL1(i) -Vapp_SL1(i) and VBS(i, j)=Vapp_BL1(j) -Vapp_SL1(i), respectively. The equivalent resistances of the cells' resistor and transistor are calculated through the differential of (4) with respect to VGS or (8) with respect to VBS. The new array resistances are then feedbacked to the code, and a new set of array WL, BL, and SL voltages are calculated. The equivalent resistance of each cell is then calculated using the new set of voltages until a solution that satisfies an acceptable error is obtained (in this work, the error is set as 10 -5 ).

Simulation results of 1T1R crossbar array
It is easy to know that a target 1T1R cell can be uniquely determined and selected by cross-connected WL and BL. The apparent resistance (Rar) of selected cell is defined as which is calculated by dividing the voltage difference between the applied potential at the selected BL and voltage at the sensing resistor (Vsens_SL_select) by the total current collected through selected SL's end. In (11), iRar represents the row index that shares SL with selected cell. Table 2 shows 12 resistance patterns, which are used to calculate the corresponding apparent resistances. These 12 patterns can be divided into two categories, one of which measures the apparent resistance of the cell closest to the voltage source and current-collecting resistor (1,1), and the other measures the apparent resistances of the furthest cell (m, n). These cases are dedicated to measuring the apparent resistance when the selected cell (Rselect) is high or low and the unselected cells (Runselect) are all high, all low, or random pattern of the two resistance states. In this work, [0] and [1] are used to indicate HRS and LRS, respectively. With reference to the carefully selected input parameters in [25], the input parameters used in this work are shown in Table 3. Unless otherwise specified, most of the simulations in this paper are based on the single-side supply voltage and single-side grounding methods in this table. Based on the input parameters given in Table 3, the apparent resistances of the 12 cases for a 100×100 array are simulated, as shown in Fig. 9 (in MATLAB and Cadence). Generally, since the length and material of the connecting line between arbitrary two adjacent cells are the same, the line resistances of any two adjacent cells along WL (RWL), BL (RBL), and SL (RSL) can be represented by the same symbol Rl. If necessary, they can also be adjusted independently. In Fig. 9, the Rar values for the cases Rselect=Rhigh are plotted in red (namely, cases 1, 3, 5, 7, 9 and 11), whereas the Rar values for the cases Rselect=Rlow are plotted in blue (namely, cases 2, 4, 6, 8, 10 and 12). Table 3 The input parameters for measuring the apparent resistances of selected cell.

Input parameters Values Input parameters Values
Transistor temperature 300K Transistor channel size 4μm/0.18μm In order to verify the reliability of the matrix algebra solution of the 1T1R crossbar array, we set the same device and array parameters to simulate the identical array in Cadence, as shown in the right column of Fig. 9. The absence of the Rar values of cases 5, 6, 11, and 12 in the right column of the figure is due to the fact that the random resistance pattern is difficult to conveniently set in Cadence. It can be seen intuitively from the figure that the corresponding Rar values are basically the same in MATLAB and Cadence. The relative error inset proves that our matrix algebra solution for 1T1R crossbar array owns high accuracy (less than 3%). The relative error calculation formula we use is ar, MATLAB ar, Cadence ar ar, Cadence For crossbar arrays, sensing margin is a key performance factor, which determines the minimum sense voltage window for reading operations. If only the LRS and HRS of the eNVM devices are used to represent the 1bit memory/weight, the sensing margin in this work is calculated from the difference between the potential on the sensing resistor of Case 10 and Case 1 relative to the selected BL's applied voltage [7,30], which is   Unlike the 1D1R crossbar array, the apparent resistances of the 12 cases in 1T1R crossbar array are higher than their preset values whether Rselect=Rhigh or Rlow. This trend will significantly reduce the sensing margin, so it seems that the purpose of using MOSFET as the controller of the eNVM device-based crossbar array has not been achieved. The accuracy of read and write operations should be as high as possible, so as to achieve a single cell multi-level memory or computation. According to existing studies, the 1T1R crossbar array can effectively suppress the sneak-path current, but we estimate that the large reading error in the aforementioned simulation is not caused by the inherent sneak-path current issue of crossbar arrays (at least not the main reason) [31][32][33]. In general, the sneak-path current issue would at least make the apparent resistance lower than the preset value when Rselect=Rhigh.  Table 3 and the 1T1R fit functions in Table 1 are used to simulate the apparent resistances of the 12 cases in ...

MATLAB, Right column: Cadence). The input parameters in
Share WL/SL with target cell Selected cell The potential of the selected SL determined by I target and R sens_SL1_select , it is easy to raise the SL potential with the increase of the R sens_SL1_select . Analyzing the architecture of the 1T1R crossbar array based on the input parameters in Table  3, we found that three factors together lead to the appearance of the backflow current, so that the measured Rar in various cases are greater than the preset resistances. The three factors are (a) the source and drain of the MOSFET can be interchanged due to the symmetry of its structure, (b) the resistance value of the sensing resistor Rsens_SL1_select of the selected row, and (c) the resistance values of the terminal resistors Rsens_BL1_unselect of unselected columns. Figure 10 is a schematic diagram for the phenomenon of backflow current in 1T1R crossbar array, which occurs on the unselected BLs sharing with selected WL/SL. In this figure, the issue of common sneak-path current in crossbar arrays is not illustrated. If the Rsens_SL1_select value is larger, the current collected by it will produce a greater voltage drop across itself. Moreover, when the terminal resistances of the unselected BLs are too small (For example, Rsens_BL1=10Ω in previous simulation.), the source potential of the controller (NMOSFET) of the unselected cells that share with selected WL/SL are likely to be higher than the drain potential. As a result, backflow current occurs, which in turn reduces the current collected on Rsens_SL1_select, leading the apparent resistance of the selected cell is higher than the preset value. However, the 1D1R crossbar array does not exist this obvious backflow current issue due to the unipolar characteristics of the selector diode. Therefore, special attention should be paid to this issue when modeling and simulating 1T1R crossbar arrays.  In view of the three factors that cause the issue of backflow current, there are two solutions that can be tried, one is to reduce the Rsens value of the selected SL row at current collect terminal, the other is to use a floating method for the unselected BL (i.e. all Rsens_BL2 and Rsens_BL1_unselect are set as 100MΩ). When other parameters unchanged, the calculated apparent resistances of the 12 cases by changing Rsens_SL1_select values are shown in Fig. 11. The results show that when Rsens_SL1_select decreases, the backflow current phenomenon becomes weaker, and when Rsens_SL1_select increases, the backflow current phenomenon becomes more obvious. This accords with our prediction, but obviously adopting this scheme can't solve the backflow current issue well (and this will make the sensing margin smaller as described later). With the unselected BL1 and all BL2 floating, the calculated apparent resistances of the corresponding 12 cases, including MATLAB and Cadence simulations, are shown in Fig. 12, both floating and non-floating schemes are shown for comparison. After adopting the floating scheme, the issue of backflow current is well solved, so that the resistance of the selected cell can be read with high accuracy regardless of whether it is HRS or LRS. If the relative error here is defined as [(Rar−Rpreset)/Rpreset]×100%, the maximum relative errors at LRS and HRS are 10.308% and 0.3887% by using matrix algebra method, respectively. In addition, regardless of whether the floating or non-floating scheme is adopted, the simulation results of 1T1R crossbar array in MATLAB and Cadence are very consistent (the Cadence simulation does not include random cases), which confirms the validity of 1T1R crossbar array modeling.  After verifying the effectiveness of using matrix algebra method to model 1T1R crossbar array, we further verify its flexibility. Figures 13 and 14 respectively show the effects of changing the gating voltages Vapp_WL1_select or Vapp_BL1_select and changing the line resistances Rl on the array read operation. It is easy to see that the changes of the voltage Vapp_WL1_select or Vapp_BL1_select have little influence on the reading accuracy for the array cell resistance, while the change of the line resistance Rl has a greater influence on the reading accuracy. It can be seen from the inset in Fig. 14 that the effect of changing RSL on the array read operation is basically the same as changing the line resistance Rl (RWL, RBL, and RSL). Therefore, it is proved here that the applied voltage on the selected WL and BL and the line resistance along WL and BL have little influence on the read operation of 1T1R crossbar array, while the line resistance along SL has a greater influence on the array read operation. The former shows that the transistor in the 1T1R crossbar array has a powerful control action, while the latter is because the line resistances on the SL are directly in series with the current collecting terminal.
The aforementioned simulations are all based on the 1T1R crossbar array with 100×100 size, and here we change the array size to evaluate the read operation, as shown in Fig. 15. It can be seen that the larger the array under the non-floating scheme, the more prominent the backflow current issue and the lower the reading accuracy. For the 200×200 and 500×500 1T1R crossbar arrays, the floating scheme is used to simulate apparent resistances, and we find that the backflow current issue is drastically reduced, whereas the array read operation is mainly sensitive to array sizes. Selecting two parameters (Rsens_SL1_select and array sizes) that have a greater impact on the read operation as X and Y coordinates, and the corresponding sensing margin as Z coordinate to simulate the 1T1R crossbar array, as shown in Fig. 16, including non-floating and floating schemes (i.e. all Rsens_BL2 and Rsens_BL1_unselect are set as 100MΩ). In Fig. 16(a), we find that the smaller the array and the larger the sensing resistance, the greater the sensing margin. However, it can be seen from Fig. 16(b) that in the case of the floating scheme, the influence of the array sizes on the sensing margin is reduced, and it has higher sensing margin under the same array size and Rsens_SL1_select as in non-floating scheme.
Looking back at (13) and combining the simulation results in Fig. 16, it can be seen that the sensing margin and Rsens_SL1_select are strongly correlated. At the same time, it is easy to find that the measure and calculate of the sensing margin using this method has certain limitations, that is, the maximum sensing margin varies with Rsens_SL1_select. For example, in an ideal situation, if Rsens_SL1_select is 100Ω, the maximum sensing margin is 1% (That is, under extreme conditions, Vsens_SL_case10=(Vapp_BL_select/Rlow)×Rsens_SL1_select=(0.2V/10KΩ)×100Ω=2mV, and Vsens_SL_case1=0, so as per (13), the maximum sensing margin is (2mV/0.2V)×100%=1%.). By analogy, if Rsens_SL1_select is 1KΩ, the maximum sensing margin is 10%, and if Rsens_SL1_select is 2KΩ, the maximum sensing margin is 20%. These simulation results and related analysis remind us that the use of floating scheme and high Rsens_SL1_select is very necessary for the read operation of 1T1R crossbar arrays.
All the foregoing simulations are performed under the conditions of single-side applied voltages of the 1T1R crossbar arrays and the single-side connection of sensing resistors. In order to further improve the accuracy of the array read operation, the identical voltages can be applied to both sides of crossbar arrays. To this end, we simulate the current maps of each 1T1R cell of the 100×100 array with randomly distributed resistance states under four applied voltage schemes to evaluate the qualitative and quantitative information, as shown in Fig. 17. The four applied voltage schemes used in this paper are shown in Table 4, and other input parameters that are not stated refer to Table 3. It should be noted that when the WL voltages and the BL voltages are single-side, it means that they are applied from the left side and the upper side, respectively. A qualitative conclusion can be drawn from the change trend of the four current maps in Fig.  17, the cell closer to the applied voltage has a larger current, and the cell farther from the applied voltage has a lower current (The premise for drawing this conclusion is to compare the currents of HRS cells and LRS cells, respectively.). That is to say, applying voltage on both sides of the array can effectively improve the uniformity of the current flowing through each cell. The applied voltages are delivered to each cell to a greater extent, thereby improving the array reading accuracy. However, from the perspective of access or calculation based on 1T1R crossbar arrays, it will lead to an increase in peripheral circuits and a larger chip footprint, which is related to a trade-off between two aspects. From a quantitative point of view, when voltages are applied to both sides of WL and BL, the ratio between the maximum and minimum current flowing through a cell is reduced from 5.95 (voltages are applied to single-side WL and BL) to 2.24 for LRS cells, while from 1.05 to 1.01 for HRS cells. It is very convenient to use the matrix algebra method to model 1T1R crossbar arrays to simulate their cell current map figure, while the use of hardware circuit simulation requires a large amount of data to be measured, and it is a time-consuming task to convert these data into a map. 10

Conclusion and discussion
Based on the matrix algebra modeling method of 1D1R crossbar arrays, this paper realizes the matrix algebra modeling of 1T1R crossbar arrays (line resistance is taken into account). The process we implemented is called transfer modeling in this paper. In view of the main difficulties in the transfer modeling from 1D1R to 1T1R crossbar array mentioned in Section 1, we give the corresponding solutions in the maintext: (i) The parameters that susceptible to the terminal voltage in the MOSFET when working in the sub-threshold region (cutoff region is merged into this region) and triode region are respectively fitted to realize the V-dependent model of the 1T1R cell; (ii) According to the location of each cell of crossbar arrays, call the function of the corresponding work region to express. The calling method of (ii) is: the selected cell calls the fitted triode equation (4), and the unselected cells (including share WL and SL, BL, and no line with selected cell) call the fitted subthreshold equation (8).
In order to verify the accuracy of the 1T1R crossbar array modeling, the modeling based on the matrix algebra method and SPICE-based circuit simulation are compared and analyzed. After that, the flexibly parameter-adjustable tool is used to simulate and analyze the quantitative effects of line resistance, array size, and voltage application scheme on the array read operation under 12 resistance patterns. In addition, this tool has been used to analyze and study the backflow current issue that easily occur in 1T1R cross arrays. Finally, simulation and analysis of the current maps with random resistance pattern under single-side and double-side voltage application schemes is carried out. As a convenient and accurate numerical analysis method, the research outcomes in this paper can provide new ideas for groups that conduct memory or computing research based on the 1T1R crossbar array architecture.
In this paper, the parameter fitting method is used to fit the MOSFET parameters under a specific process to make the modeling based on the matrix algebra method more accurate. By referring to this method, the fitting of MOSFET parameters under other processes can be realized. The HRS and LRS of the eNVM device used in the simulation are the typical expected values of GeSbTe-type PCRAM [34,35], and the resistance values of other eNVM devices can be set when needed. It should be noted that the proposed method can only analyze 1T1R crossbar arrays under steady state, and dynamic analysis is beyond the scope of this paper. The 1T1R crossbar array can effectively suppress the sneak-path current, so that the eNVM device-based 1T1R crossbar array may achieve multi-bit memory or computing in each cell. The matrix algebra-based numerical analysis of the multi-bit or even continuous resistance for array operation will be carried out in our further work.
Appendix Due to the difference between 1T1R and 1D1R array architecture, the matrix calculation method of [24] is used for reference and modified below, and the related matrixes used for calculation are listed one by one.
The 2mn unknown voltage variables VBL(i, j) and VSL(i, j) in a m×n crossbar array can be written as a 2mn×1 column vector thus (10) Similarly, the vector E can also be split into smaller matrixes to facilitate understanding and programming. The corresponding splitting process and expansion equations are shown in (A8)~(A10). With these equations, the 2mn voltage variables and related electrical parameters can be solved in MATLAB. B1 S1  Schematic diagram of a model with a NMOSFET in series with a binary resistor.

Figure 3
Calculated ζ-VGS data points and its t curve, for VGS(0, 0.55V) with step 0.025V.  Schematic diagram of 1T1R crossbar array with m WLs and SLs (horizontal line) and n BLs (vertical lines). If the upper left cell is selected, the other cells can be divide into three groups: those share WL and SL with selected cell, those share BL with selected cell, and those share no line with selected cell.

Figure 7
Schematic diagram of Kirchhoff's law for electrical parameters at each junction in simpli ed 1T1R crossbar array with BL and SL line resistances.

Figure 8
Schematic diagram of the de ned input parameters for 1T1R crossbar array.

Figure 9
Caculated apparent resistances of the same crossbar array in MATLAB and Cadence (left column: MATLAB, Right column: Cadence). The input parameters in Table 3 and the 1T1R t functions in Table 1 are used to simulate the apparent resistances of the 12 cases in Table 2. The array size and line resistance are 100×100 and Rl=5Ω, respectively (red: high cases, blue:low cases). Inset: the relative error of apparent resistance of cases 1~4 and 7~10 in MATLAB and Cadence simulations.

Figure 14
Calculated apparent resistance values for a 100×100 1T1R crossbar array by changing line resistance Rl. Inset: Calculated apparent resistance values for a 100×100 1T1R crossbar array by independently changing the line resistance RSL.

Figure 15
Calculated apparent resistances for 1T1R crossbar array with different array sizes under non-oating and oating schemes. Figure 16 3D plots of sensing margin for 1T1R crossbar array under different Rsens_SL1_select and array sizes adopting (a) non-oating and (b) oating schemes.