MPT Tool: A Parameter Extraction Tool for Accurate Modeling of Magnetic Tunnel Junction Devices

Until now, many models of magnetic tunnel junction (MTJ) devices have been reported in the literatures, which are very helpful for behavior simulations in spintronic circuit design and simulation. However, there are few reports on how to obtain the accurate modeling parameters, as well as the accurate simulation of the complete flow from device parameters to device modeling and hybrid MTJ/CMOS circuits. One of the difficulties in obtaining the accurate MTJ modeling results is the lack of experimental parameters. The parameters in the MTJ models should be precisely identified to obtain accurate simulation results. However, in the reported models, many parameters are theoretically set instead of actual extraction. In this paper, a parameter extraction tool, named MPT tool, is proposed for the accurate modeling of the MTJ devices. In the MPT tool, the measured results of the MTJ devices are adopted as the input. There are two main functions of the MPT tool. One is the parameter extraction based on the input experimental results. The other is the parameter optimization based on the specific algorithm. In the tool, the simulated annealing (SA) optimization algorithm is applied to optimize the extracted parameters. The parameters of the extraction procedure and the modeling results based on the parameters of the optimization procedure are compared with the experimental data, respectively. The proposed extraction procedure and the optimization procedure together with the adopted MTJ model is as an efficient approach to improve the accuracy of the simulation result.


D
MTJ device diameter (nm) t FL MTJ free layer thickness (nm) t ox MTJ oxide barrier thickness (nm) RA Resistance-area product ( µm 2 ) TMR Tunnel magnetoresistance ratio etc. [1]. As the key element in each MRAM cell, magnetic tunnel junction (MTJ) device plays important role in the memory. Currently, perpendicular magnetic anisotropy MTJ (P-MTJ) shows better performance than in-plane MTJ (I-MTJ). P-MTJ device is considered as the most promising MTJ device with high TMR ratio, high density, and low switching threshold current [2]. Many kinds of P-MTJ model have been proposed, as listed in Table 1. The models in Table 1 is divided into two categories based on the magnetic behavior description. One is the static magnetic model (SM), in which the electrical switching behaviors and the corresponding MTJ resistances, switching delay time are included. The other is the dynamic magnetic model (DM), in which the magnetization vector is keeping tracking on time [3].
With the developments and improvements of the MTJ models, various variables that affect the characteristics of the MTJ device are included, including voltage, temperature, size, structure, pulse width, stochastic initial magnetization angle, and stochastic thermal fluctuations, etc.. Based on Table 1, it can be seen that most MTJ models are based on Verilog-A implementation, by which the physical equations are expressed directly. Some models are based on SPICE circuits with faster simulation speed and low memory cost. In addition, several MTJ models are implemented in other ways, such as VHDL-AMS, MATLAB, finite element simulation and Fokker-Planck simulation framework, etc. However, concerning the problem that the parameters used in the models mainly are set theoretically, a method for determining model parameters should be studied. The accurate modeling parameters can be combined with the physical models to reflect the characteristics of the MTJ device, which would be very helpful for the designers in circuit and system simulation. It should be mentioned that the accurate modeling parameters are not the accurate physically parameters to some extent. There are still undiscovered physical phenomena in the device models, which lead to a gap between the model and device. To cover the gap, a tool for the extraction of the modeling parameters, need to be developed. In the tool, two steps are needed, including the parameter extraction and the parameter optimization. The extraction step is based on the experimental results of the device. Because of the undiscovered physical phenomena of the device, based on the device model with the current extracted parameters, there is obvious deviation between the experimental results and the simulation results. To solve this problem, unique optimization step is needed, by which the parameters are optimized to match the embedded model. In this way, the embedded device model combined with the optimized parameters can be used to simulate the device, showing less deviation from the experimental results. In this sense, the parameters obtained by the optimization are referred as "the accurate modeling parameters". Fig. 1 shows the flowchart of the hybrid MTJ/CMOS circuit simulation. Besides the standard CMOS simulation platform, such as Cadence or HSPICE, the simulation of the MTJ device should be applicable to different circuits  and integrated with the CMOS process. To conduct the MTJ-related simulation, both the MTJ device model and the modeling parameters are needed to guarantee a reliable simulation result.
For the MTJ device models, as shown in Table 1, many groups have reported the physical mechanisms related to the device behavior. On the other hand, the key parameters in the models should be extracted based on the experimental results. In this way, the combination of the physical model and the extracted parameters can be helpful for the MTJ device simulation. The MTJ parameter extraction tool, named MPT, is developed based on MATLAB in the paper. Firstly, the initial parameters are extracted based on the MTJ device experimental data, which are then applied in the developed MTJ model to get the initial simulation results. The initial simulated results are compared with the experimental data to check its deviations. The optimizations of the extraction parameters are then conducted based on specific algorithm to minimize the deviation. Finally, the optimized parameters are used to the simulation platform for the accurate simulation results. In this manner, with the help of the parameter extraction and optimization, the circuit simulation results would be more reliable than the initial simulation results. Fig. 2 shows the flowchart of the developed MPT tool. Firstly, the necessary modeling parameters including structural, electrical and magnetic parameters, are extracted initially by fitting the experimental MTJ device data. Then, these initial parameters are input into the MTJ model for the initial device simulation. Comparison between the initial simulation results and the experimental data is made to find the deviations. Based on the observed deviations, the objective function f obj is generated to evaluate the modeling accuracy. The simulated annealing (SA) algorithm is used to optimize the parameters to get more consistent modeling results. The optimized modeling results with the optimized parameters are then compared with the experimental MTJ data again to verify the developed extraction tool. In this case, a closed-loop for the parameter extraction is presented.
Some of the parameters of the MTJ device can be directly obtained by the measurements, including the resistances (R P and R AP ) and switching current, etc. However, the other parameters, such as the saturation magnetization M S , the energy barrier E b , the shape anisotropy, the interface anisotropy, etc., are also important for the simulation. But their corresponding values can't be obtained experimentally. Currently, the theorical values are always used in the simulation if the parameters can't be measured directly. This is the main reason that there exists obvious deviation between the SPICE simulation and the measurement results of the characteristics of the MTJ device.
To solve the problem, the parameter extraction approach is proposed in the paper, aiming to obtain the values of the parameters that can't be measured directly. To satisfy the requirement of the accuracy of the simulation, the parameter extraction approach is divided into two steps. Firstly, parameter extractions are conducted based on the experiment results. Then, the optimizations on the extracted parameters are followed. Finally, the optimized parameters could be used in the simulation and the simulation results coincide with the experiment results very well.
The rest of the paper is organized as follows. Section II introduces the MTJ model strategies, while Section III reports the accurate modeling process with the parameter extraction and the parameter optimization. The evaluation of the accurate modeling simulation is also presented in Section III to verify the developed accurate model by the proposed MPT tool. Finally, conclusions are made in Section IV.

II. P-MTJ DEVICE MODEL
There are two approaches to model the magnetic dynamic switching characteristics of the P-MTJ device, as listed in Table 1, including the SM and the DM. For the SM approach, the analytical expression of the switching threshold is derived from Landau-Lifshitz-Gilbert (LLG) equation, and the switching time is also included in the model. For the DM approach, the LLG equation is included in the model directly, which means the real-time magnetization simulation is conducted [24]. Although the LLG equation is unable to model the thermal activations, researchers have developed the dynamic model by adding the stochastic thermal fluctuation term to account for thermal effects [3], [14], [25].
The division of SM and DM in Table 1 is not based on the electrical dynamic properties. It is based on the magnetism transfer properties of the MTJ device. In the other words, SM models in Table 1 refer to the "magnetism static model", in which the real-time magnetism propagation process isn't included in the model. On the contrary, the DM models in Table 1 refer to "magnetism dynamic model", in which the magnetism propagation process is simulated by LLG equation. The SM models in Table 1 are commonly used in the electrical circuit simulation with the benefits of the fast simulation and high accuracy. It can be used to simulate both the electrical static properties and the electrical dynamic properties, including the electrical switching behaviors of the MTJ device. So the MPT tool can be used both on SM and DM models, providing an effective approach for the accurate modeling of the static and dynamic electrical simulations of the circuit.
The obtained modeling parameters in the tool are suitable for the adopted model to get accurate simulation results.
In fact, regardless of the accuracy of the model, the steps of parameter extraction and optimization are required for the accurate simulation of the fabricated MTJ devices. The theoretically set reasonable parameters can only obtain theoretically reasonable simulation results. The purpose of the tool is to obtain accurate circuit simulation results that can be adapted to the embedded model rather than to obtain parameters that accurately describe the physical information of the MTJ device.
The key equations of the embedded MTJ model are listed below: Equations (1) and (2) present the resistance model, in which R P and R AP represent the MTJ resistance in parallel (P) state and anti-parallel (AP) state, respectively. TMR is the tunneling magnetic ratio, which is dependent on the voltage across the MTJ device. Equations (3) and (4) can be used to calculate the critical current density J C0 and the switching current density J SW , in which is the thermal stability factor, [4], and E b = −δ N M 2 S t FL /2μ 0 + K b t FL + K i Area for perpendicular shape anisotropy (PSA) MTJ [26], where δ N is the shape anisotropy coefficient. g is the spin polarization efficiency factor, g = g SV ± g tunnel , where When the forward current flows through MTJ with J > |J P SW |, the P to AP switch occurs, with the resistance R = R AP . When the reverse current flows through MTJ with |J| > |J AP SW |, the AP to P switch occurs, with R = R P . μ 0 = 4p ×10 −7 H/m the permeability in free space, For spin orbit torque (SOT) MTJ device, by adding a heavy metal (HM) layer on the P-MTJ, the static model of SOT-MTJ device has a high degree of coincidence with (1)(2)(3)(4)(5). However, the switching current in the SOT-MTJ device is different from the STT-MTJ device. The current in SOT-MTJ device flows through the HM layer, while the current in STT-MTJ device is directly flows through the MTJ device. The relationship between the driving current density J c in the HM layer and the induced spin current density J s in the P-MTJ is as follows [27]: where θ SH the spin hall angle. The measurement method of θ SH is given in [28]. The total resistance of the SOT-MTJ device includes the resistances of the HM layer (R HM ) and the MTJ device (R P or R AP ). The resistance of the HM layer can be expressed as: where ρ the resistivity of HM, and L HM , W HM , d HM are the length, width, height of the HM layer, respectively. It can be seen that the main difference between SOT and STT devices of the static model is the existence of θ SH and R HM in the SOT device. The methodology of the proposed MPT tool can be also applicable to the SOT-MTJ device. The simulation on the SOT-MTJ device is not included in the paper because of the shortage of the SOT-MTJ experimental data. Fig. 3 shows the flow of the parameter processing in the MPT tool. The initial extracted parameters are obtained based on the uploaded experimental data of the MTJ device. It should be noted that the initial structural parameters, including D, t ox , t FL , are directly input in MTP tool.

III. ACCURATE MODELING BY PARAMETER EXTRACTION AND PARAMETER OPTIMIZATION
With the optimization process, new parameters are obtained. The SA optimization is completed in two steps based on f obj , where the extracted parameters are as the input initial parameters. The first SA optimization is conducted with symmetrical operation for the parameters affecting the resistance. The second SA optimization is conducted for the parameters affecting the switching current density, where the updated parameters obtained in the first SA and part of the extracted parameters are as the input initial parameters.

A. PARAMETER EXTRACTION
In order to present the practical application of the developed MPT tool, it is applied with the experimental data of SB-MTJ device in [29] and the experimental data of PSA-MTJ device in [26]. With the experimental data uploaded, extraction process is performed in the MPT tool.
The parameter resistance area product RA is extracted by linear fitting according to the linear function of D TEM and G 0.5 described as following: where D TEM is the physical diameter of MgO barrier layer observed by transmission electron microscopy (TEM), G is the measured conductance, D 0 is the length of the dead zone. Fig. 4 shows the D TEM − G 0.5 data [26] and its linear fit, D 0 is obtained from the x-intercept, and RA is extracted from the slope of the fitted line. The parameters TMR 0 , V h , P can be extracted based on J-R experimental data, as shown in Fig. 5. Fig. 5 (a) shows the related results of SB-MTJ device, while Fig. 5 (b) shows the results of PSA-MTJ device. Due to the difference in V h , the J-R curves behave differently. The larger the V h , the smaller the influence of the voltage V on R AP . The blue dots in V-TMR coordinate is extracted by TMR = (R AP − R P )/R P , where R P and R AP are obtained by symmetrical processing. The red line in V-TMR coordinate is fitted by the TMR model [4], [6]. As a result, TMR 0 and V h can be obtained by solving Eqn. (1). P is extracted by [31].
The parameter saturation magnetization M S is extracted from the in-plane magnetization data and the out-of-plane magnetization data. Fig. 6 (a) and Fig. 6 (b) show the inplane and out-of-plane magnetization data from [26] with various t FL , respectively. It can be seen there are four segments of saturation region for each t FL . The saturation value (m S ) is extracted by taking the average of the absolute values of the magnetization in the four regions. The extracted t FLm S and its linear fitting are shown in Fig. 6 (c). M S is equal to the slope of the fitted line.   The parameter interface anisotropy energy density K i and the parameter bulk anisotropy energy density K b can be extracted by linear fitting based on the t FL -Kt FL data, where K is the anisotropy constant. Fig. 7(a) shows the t FL -Kt FL data and its linear fitting. Due to K = K b − M 2 S /2μ 0 + K i /t FL [29], the values of K b − M 2 S /2μ 0 and K i can be extracted by the slope and the intercept to the vertical axis of the blue line in Fig. 7(a), respectively. Based on the ferromagnetic resonance (FMR) measurements, the effective anisotropy field H Keff and the damping factor α can be extracted. Fig. 7 (b) shows the θ − μ 0 H R curve, where θ the magnetic field angles, H R the resonance field. H Keff is extracted by solving a binary linear equation [32]: where f the frequency. Considering (9) as a binary system of linear equations about f and H Keff , H Keff can be obtained by solving the equations. Fig. 7(c) shows the θ − FWHM (full-width at half-maximum) curve, the magnetic damping constant α can be extracted [26]. Based on the above extraction approaches, the MPT tool is implemented to carry out the function of the parameter extraction.

B. PARAMETER OPTIMIZATION
Although the MTJ device model is gradually matured, the reliable simulation results cannot be obtained without the input of accurate modeling parameters. Based on the abovementioned parameter extraction approach, the extracted parameters are not accurate enough for the circuit simulation. There is obvious deviation of the simulation results based on the extracted parameters from the experimental results. Two main reasons accounting for the deviation are discussed. The first one is the deviation from the extraction process. The second one is the deviation of the embedded model, there are still many undiscovered phenomena in the nanoscale spintronic device. The existence of the gap between the device model and the actual device should be solved during the simulation. The optimization of the MTJ parameters is adopted in the paper to minimize the gap. The deviation between the modeling results and the original experimental data is calculated by the objective function [33]: where N the number of the data, g s (R i ) is the simulated data and g e (R i ) is the experimental data.
The objective function f obj is used to evaluate the goodness of the modeling result. In other words, it is used to evaluate the deviation between the modeling results and the experimental results. Therefore, searching the minimum value of f obj under nonlinear constraints is the principle of the optimization. For complex problem involving large space, nonlinearity, global optimization, and combinatorial optimization, etc., there are many available intelligent optimization algorithms, such as genetic algorithms (GA), particle swarm optimization (PSO), simulated annealing (SA), taboo search (TS) and artificial neural network (ANN), etc [34]. In this paper, the SA algorithm is adopted since it can provide effective approximate solution for non-deterministic polynomial complexity problems [35].
The SB-MTJ device is used as the example to show the process of the parameter optimization. The optimization process of the PSA-MTJ device is similar. In the MTJ model, the parameters related to the MTJ resistance include D, t ox , TMR 0 and V h . And the parameters affecting the switching current J SW include t FL , α, M S , H Keff and P. The SA algorithm is performed in two steps with f obj , firstly the parameters affecting the resistance is optimized, and then the parameters affecting the J SW is optimized. In this way, the resistance and the J SW are adjusted. The concrete steps are as follows. Firstly, the searching algorithm is performed for the values of D, t ox , TMR 0 , V h with the extracted parameters input as the initial parameters. Afterwards, the second SA is performed for the values of t FL , α, M S , H Keff and P, where the searched D, t ox , TMR 0 , V h and the extracted t FL , α, M S , H Keff , P input as the initial parameters. Fig. 8 shows the process of applying SA algorithm to search for the minimum value of f obj , in which the initial extracted parameters act as the input parameters for the SA process. The new parameters are generated based on the initial parameters, and the new f obj is also calculated for comparison, synchronously. Then a judgment is made as to whether to accept the new parameters as a new solution. If accepted, keep the new parameters firstly, and then judge whether the number of updating parameters reaches the set number of iterations. If the number of iterations is also satisfied, the current parameters are taken as the final parameters after optimization and saved in the tool. If not accepted, then the parameters are updated again, and the trial-anderror process of "update parameters -calculate f obj -retain or discard the solution" is repeated until the iteration ends. The update of parameters is based on the random selection of points in the search interval, and the Metropolis sampling rule is used in SA algorithm to make the random selection gradually converge to the local optimal solution. Fig. 9 shows the relationship between the f obj value and the iteration times of the second SA searching process. As shown in Fig. 9(a), due to the uncertainty caused by nonlinearity, the current f obj value fluctuates largely even with 4000 times of iteration operation. Fig. 9(b) shows the searching process of the minimum value of f obj .  Whichever the algorithm is selected for the optimization, its speed is directly affected by the defined parameter variation ranges. In the SA algorithm, the initial values of the parameters show much more influence on the optimization speed than the variation ranges of the parameters. Fig. 10 shows the optimal searching processes by using the GA and SA algorithm, separately, with the same variation ranges of the parameters. Compared with the SA algorithm, the application of the GA algorithm doesn't need the inputs of the initial values of the parameters, and the iteration can be terminated by its own judgment. As shown in Fig. 10(a), the optimization process of the GA algorithm is automatically terminated at the 54th iteration step. The iteration totally takes 12.9 minutes in Fig. 10(a), in which the first 20 iterations cost 5.1 minutes. The SA algorithm takes about 100 iteration steps to find the optimal solution, as shown in Fig. 10(b). To verify the stability of the process, 300 iteration steps are used for the optimization. The iteration time of the SA is much faster than that of the GA algorithm. The total time of the 300 iterations only takes 1.9 minutes.
Based on the above comparison, the GA algorithm shows the advantage of the self-termination, with the shortage of the slow iteration speed. However, the self-termination depends on the variation ranges of the parameters closely. If the experimental values are scattered in a relative wide range, which are usually observed in lab, the minimum value of f obj can't be obtained easily by the GA algorithm. So, the SA algorithm is adopted in the optimization, with the benefits of its fast speed and high tolerance to the wide scatter ranges of the experimental data.  Fig. 11 (a) shows the results before the optimization, while Fig. 11 (b) shows the results after the optimization. The obvious deviation between the experimental results and the simulated results is presented in Fig. 11 (a), where the experimental data is from [29]. The parameters in Fig. 11(a) are the initial extracted parameters without optimization. Fig. 11 (b) shows the simulation results with the optimized parameters. After the optimization, the value of f obj is decreased from 0.3231 to 0.0022. The embedded model in Fig. 11 (a) and Fig. 11 (b) are the same. The difference mainly comes from the parameters.

C. MODELING RESULT AND EVALUATION
For the PSA-MTJ device, the results are shown in Fig. 12, in which the experimental data is referenced from [26]. The structure of the device as shown in [26] is: thermally oxidized Si substrate/Ta ( [26], as D = 10.4 nm, t FL =15 nm, t ox1 = 0.93 nm, t ox2 = 0.90 nm, as shown in Fig. 12(a). After the optimization, the value of f obj is decreased from 0.4569 to 0.0596. In addition, we also compared our result with others in following to further verify the improved accuracy with the approach in the paper. For P to AP switch, the experimental data of J SW is ∼3.2× 10 11 (A/m 2 ) for the 10.4 nm PSA-MTJ device in [26], and the theoretical calculated J SW by the model in [26] is 4.4×10 11 (A/m 2 ). The J SW in [11] is ∼2.0×10 11 (A/m 2 ), in [23] is 2.68×10 11 A/m 2 , and it is 3.3049×10 11 (A/m 2 ) in this paper. Fig. 11 and Fig. 12 confirm the usefulness of the proposed MPT tool. The precise identification of the modeling parameters allows adopting it to design and simulation of MRAM circuits. One thing should be mentioned that the variation range of parameters should conform to the physical meaning. The obvious inconsistency between extracted parameters and optimized parameters indicates that there is some undiscovered physics for the PSA-MTJ device operation and should be investigated in future. Fig. 13 shows the simulation results of magnetism dynamic model, in which the propagation magnetic vectors are simulated based on LLG equation. m denotes the unit magnetic vector of free layer (FL), while m x , m y , m z the x-component, y-component, z-component of m, respectively. It should be noted the unit magnetic vector of pinned layer (PL) m p = (0, 0, −1) in the embedded dynamic MTJ model. m z witches from near −1 to 1 denotes MTJ switches from P to AP. The switching time is longer after parameter optimization, which is consistent with the obviously increased J SW in Fig. 12(b).

D. HYBRID MTJ/CMOS CIRCUIT SIMULATION
Finally, in order to verify the validity of the developed MPT tool, a hybrid MTJ/CMOS circuit is simulated based on TSMC 65-nm CMOS technology on HSPICE platform. As shown in Fig. 14, the state of the MTJ device is switched by changing the voltage of BL and SL on the two ends of the MTJ device. State "0" refers to low resistance state, while state "1" refers to high resistance state. It can be seen that the state is switched in correct logic. For the SB-MTJ  device, Fig. 14(a) shows the comparison of circuit simulation results with modeling parameters before and after optimization. For the PSA-MTJ device, Fig. 14(b) shows the circuit simulation results with parameters before and after optimization, respectively. There is obvious difference in the transient properties before and after the parameter optimization, both for SB-MTJ and PSA-MTJ devices. After the parameter optimization, the switching current as well as the switching time are different. These differences are critical for the evaluation of the MTJ device.
It should be noted that the simulation results of the circuit, as shown in Fig. 14, can't be used to demonstrate the improved accuracy by using the MPT tool on the parameter extraction. It can only indicate that there exists obvious difference before and after the application of the MPT tool.
Here the authors take the chance to be open for any possible support of the relevant experimental data from the Foundry factories, which can be very helpful for the updating of the MPT tool.
The comparison of the circuit simulation results before and after parameter optimization in Fig. 14 is consistent with the comparison of the MPT tool in Fig. 11 and Fig. 12. For example, the modeling results of SB-MTJ after parameter optimization show larger resistance than that before optimization, as shown in Fig. 11(a) and (b). Correspondingly, the MTJ current I mtj after parameter optimization is smaller than that before optimization, as shown in Fig. 14(a). The modeling results of PSA-MTJ after parameter optimization show longer switching time than that before optimization, as shown in Fig. 13. The results in Fig. 14(b) correspond to the prolonged switching time. The modeling of the MTJ device based on the MPT tool can be helpful for the spintronic circuit design and simulation.

IV. CONCLUSION
With the development of the MTJ device, various efforts are contributed to the device model, which should be accurate for the device modeling. However, the input parameters of the model are paid less attention until now. Accurate parameters are also the critical part of the model. Besides the resistance and the switching behavior of the MTJ device are concerned in the application, the influences of the modeling parameters on MTJ characteristics are more concerned by researchers. The accurate parameters are especially helpful to develop new MTJ devices and investigate the undiscovered device physics.
The MPT tool for the parameter extraction of the MTJ device is developed in the paper, which can simulate the switching characteristics of the MTJ devices with variable parameters. Two types of the MTJ devices are studied to verify the MPT tool, including SB-MTJ and PSA-MTJ. The parameter extraction and validation are conducted by comparing the experimental and modeling results. The good agreement is observed after optimization, which confirms the usefulness of the proposed extraction procedure together with adopted MTJ model as an efficient tool to help design of MRAM systems. Besides the above two kinds of STT-MTJ devices, the parameter extraction approach is also suitable for the other MTJ devices, such as SOT-MTJ device. It has potential application for the development of the MTJ device and the simulation of the hybrid MTJ/CMOS circuit applications. The MPT tool could be helpful for the Foundry factory to provide accurate modeling parameters for the design house.
The operation of the MPT tool is based on the principle of the minimum deviation between the simulation and the experimental results. Sometimes there exists obvious difference between the optimized parameters and the physical parameters of the device. This hints some undiscovered physics related to the inconsistency and should be investigated. Currently, the MPT tool is independent of the circuit simulation platform. In future, it can be applied directly on MRAM technology in Foundry factories for accurate parameter extraction for device simulation.