Modeling of Static Negative Bias Temperature Stressing in p-channel VDMOSFETs using Least Square Method

Negative bias temperature instability (NBTI) is a phenomenon commonly observed in p-channel metal-oxide semiconductor (MOS) devices simultaneously exposed to elevated temperature and negative gate voltage. This paper studies threshold voltage shift under static stress associated with the NBT stress induced buildup of both interface traps and oxide trapped charge in the commercial p-channel power VDMOSFETs IRF9520, with the goal to design an electrical model. Change of threshold voltage follow power law tn, where parameter n is different depending on the stressing phase and stressing conditions. Two modeling circuits are proposed and modeling circuit elements values are analyzed. Values of modeling circuits elements are calculated using least square method approximation conducted on obtained experimental results. Modeling results of both circuits are compared with the measured results and then further discussed.


Introduction
As the device dimensions in CMOS technologies have been continuously scaled down, a phenomenon called Negative Bias Temperature Instability (NBTI) has gained in importance as one of the most important degradation mechanisms. Degradation of transistor parameter values due to NBTI has emerged as a major reliability concern in current and future technology generations, especially in p-channel MOSFETs [1][2][3]. During the stress period, the transistor parameters slowly deviate from the nominal value. Longer the stress period, higher is the impact of NBTI on transistor parameters. NBTI effects are manifested as the changes in device threshold voltage (V T ), transconductance (g m ) and drain current (I D ), and have been observed mostly in p-channel MOSFETs operated under negative gate oxide fields in the range 2 -6 MV/cm at temperatures around 100 °C or higher [1][2][3][4][5]. Change in these parameters is dependent on the stress parameters (time, temperature, gate voltage). Considering the effects of NBTI related degradation on device electrical parameters, NBT stressinduced threshold voltage shift (∆V T ) seems to be the most critical one [6,7].
Despite notable scale down in the device dimensions, ultra-thick gate oxide reliability studies are still of significance because of broad use of MOS technology. Taking this into consideration, our earlier papers were dealing with NBTI in p-channel power vertical double diffused MOS (VDMOS) transistors [6,[8][9][10][11][12][13]. Their reliability has been investigated under various stress conditions, such as irradiation, high electric field, NBTI, NBTI under low magnetic field, and NBTI and irradiation [13][14][15][16]. NBTI is critical for normal operation of power MOS-FETs since they are routinely operated at high current and voltage levels, which lead to both increased gate oxide fields and self-heating, thus favor NBTI.
First thorough explanation of processes on negative bias temperature instabilities was made by Jeppson and Svensson [17]. Even though more than 40 years have passed since then, many mechanisms of NBTI are not very well understood yet. In the last decade, many different working groups are addressing NBTI effects, with accent on both description and modeling of voltage threshold shifts [19,20]. Swami presented a model for nano MOSFET for FinFET technologies [20], while Aleksandrov reported a model that is based on a reaction-diffusion principle [21]. Still, an appropriate electrical model to describe instabilities corresponding to different stressing conditions is lacking [22]. A model by Danković is RC based model and describes both static and pulsed NBT stressing [23]. However, further mathematical improvements are needed in order to overcome time-bounded flaws of the model, and to make use of the dependencies between temperatures and bias values. Ma presented a model that attempts to explore these dependencies [24], although they aren't analyzed on the accelerated NBT stressing conditions. Maricau used similar approach of RC based model [25], but it analyses short stressing periods, while this paper focuses on modeling of the effects caused by longer stressing periods.
The use of this type modeling is to create a model that can be mathematically calculated in order to incorporate the model into SPICE. That will enable the design-ers to consider these instabilities of the circuit during the design phase.
The primary part of this study is to propose an equivalent electrical circuit for modeling of ∆V T based on experimental data using least square method. In the section 2, experimental setup and results will be briefly described and discussed. Section 3 deals with the modeling approach, and evaluates the model and gives comparison between measured and modeled results.

Experimental setup
Used devices for this research are commercial p-channel power VDMOSFETs IRF9520 with initial threshold voltage of V T0 = -3.6 V [26]. Devices are built in standard silicon-gate technology with 100 nm thick gate oxide and packed in plastic TO-220 packages.
To properly investigate NBTI in p-channel power VD-MOS transistors in which gate oxides thickness are 100 nm, special stress voltages are needed. These voltages need to be several times larger than typical operating voltage of these transistors (more than -40 V). To ensure these specific signals, it is needed to develop an additional separate circuit that provides appropriate higher stress voltages for stressing. In years during the research, we have developed a system that automates both NBT stress and measurement on p-channel power VDMOS transistors [12].
System consists of high voltage stress circuit and of low voltage measurement circuit that are controlled with software-controlled switches. Stressing circuit includes an external amplifier between the stress voltage source unit and the device under test (DUT). The transfer I-V characteristics were measured at the drain voltage value of 100 mV, so the device was kept in the linear region of operation. Gate voltage was swept from −2 to −4.75 V, with −50 mV step. Designed system allows full range measurement of transfer I-V characteristics, that is used to extract threshold voltage using second derivative method [27]. This measurement method has previously been used by several research groups [9-13, 15, 16].
Several sets of p-channel devices were tested under different conditions. Devices were subjected to stress for 24 hours during which 36 interim measurements were performed. Two different negative voltages were applied to gate (-45 V and -50 V) while source and drain terminals were grounded. Experiments were done at two different temperatures (150 °C and 175 °C) and values of threshold voltage shifts are obtained. Our ear-lier experiments included testing devices under many different conditions in terms of both gate voltage and temperature [9][10][11][12][13]. Experimental results are given in Figure 1. For the comprehensiveness of the proposed model, only four combinations of stress conditions are shown. NBT stress under static conditions results in notable threshold voltages shifts. These changes become more pronounced at higher stress voltages and temperatures, which is in line with other investigations [18,22,28,29]. A lot of results, including ours indicate that ∆V T saturates with increase in stress time [4,6,11].
NBT stress causes threshold voltage shifts with widely different rate in different time periods [28,29]. Evolution of ∆V T through time is presented in Figure 1 Inserted graphic, and given with power law (t n ). During this evolution, distinct phases can be distinguished [4,13,30]. For every phase, the value of parameter n is different. Through each of the phases, parameter n is as close as constant (not exactly constant, but in the very limited range). During our earlier researches that involved extensive NBT stress [4,6,8,[9][10][11][12][13], in the case of long-term experiments, three different phases are detected in evolution of n [6]. Starting phase, where n = 0.4. In this phase, n is highly dependent on the stressing temperature and bias [4,6]. Second phase, where n drops to n = 0.25, and is almost independent on stressing conditions. Third, final phase, where n in again dependent on the stressing conditions and steadily declines to n = 0.14, continuing to saturation. Similar results are obtained for all stressing temperatures in the experiments. All of the results are suggesting similar development of the parameter n [6,13,31]. So, it can be concluded that parameter n is overall higher at the start phase of the stressing, but tends to saturate in later phases of stressing.
Described progress of the parameter n is caused by originated oxide trapped charge and interface traps, which directly influence the changes in the threshold voltage during NBT stress [4]. These occurrences are product of numerous electrochemical processes and reactions concerning oxide and interface defects, holes and other and species related to hydrogen. Depending on the stressing conditions and the number of defects, these reactions can occur in either forward or reverse direction. Since reversed reactions are characteristic for recovery of the degradation that occurs during pulsed stressing, in the case of static stress, forward reactions are dominant [4]. Starting phase of stressing, which explain creation of interface traps is explained with the reaction: The released hydrogen atoms ( • ) H are highly reactive, and they also can dissociate the SiH bonds at the interface or in the oxide near the interface. These reactions lead to creation of additional interface traps or positively charged oxide defects.
Released unstable hydrogen atoms react with the holes to form ions.
However, hydrogen ions also dissociate SiH bonds in the oxide near the interface leading to creation of positively charged oxide defects.
Oxide trapped charge and interface traps buildup described in the given reactions is notably enhanced in the early phase while concentration of SiH trap precursors is still high. As the stress time increases, the number of both positively charged defects and interface traps is getting higher. However, probability for reverse reaction (passivation processes) occurring rises as well.
The H 2 molecules released in reactions (2), (3) and (5) diffuse deeper into oxide and can be cracked at positively charged oxide traps: As a product of reaction (6), H • is released. It can take part in either forward reaction or reverse reaction, thereby rounding the chain of reaction. Increased amount of reverse reaction occurences lead to change in the slope of a function interpreting parameter n. Key step for appropriate modeling is to tackle the change of parameter n in phases, in order to follow the evolution of ∆V T .
Since the change of ∆V T is given with the power law (t n ), a capacitor C charged through resistor R is chosen for the central element of the modeling circuit. Capacitor charging equation is given with: Capacitor is chosen because the capacitor voltage is given in exponential form, which is a type of the power law, needed for ∆V T modeling. So, the capacitor voltage, V C models ∆V T . In given modeling circuit, rise of the capacitor voltage should correspond to the rise of the ∆V T , so that in any moment t 1 , V C should be as much as accurate as ∆V T . To accomplish that, specific controlled charging rate of the capacitor must be achieved. Value of time constant, τ, must be calculated first, and then values of capacitance of capacitor C and resistance of the resistor R must be fitted. Value of V S is acquired using stretch exponential (SE) equations [3,6] and listed in Table 1. The SE fit predicts that value of ∆V T will saturate after extended stressing and it estimates the saturation value. Increased stressing time leads to better estimation of the saturation value, as can be seen for parameter n. Stretch exponential equation is given with: In the equation (8), β, τ 0 and ∆V Tmax are fitting parameters. Parameter β is defined as a distribution width, and τ o represents a characteristic time constant of the distribution. Parameter ∆V Tmax is a value of ∆V T saturation [6,42]. Value of V S is actually value of ∆V Tmax given in equation 8. Through this value, dependence of the modeling results on bias and temperature values is given. Therefore, for different stressing conditions, value of V S is different as well. Although at first this looks as a serious limitation, based on our earlier studies, it is possible to estimate the interdependence of time, voltage, temperature and ∆V T of investigated VDMOSFETs using the results obtained by accelerated NBT stressing [39,40].
To do modeling based on experimental data, it is needed to find a function that describes experimen-tal data sets, and then to calculate modeling circuit element values based on a fitting function parameter. Most appropriate method for this fitting is Least Square Method (LSM) [41]. LSM is one of the most widely used method to find or estimate the numerical values of the parameters. It can also be used to fit a function to a set of data and to characterize the statistical properties of the estimates [42][43][44]. LSM is a mathematical method for finding the best-fitting curve to a given set of points by minimizing the sum of the squares of the offsets of the points from the curve. Basic principle of the LSM implies that a set of I pairs of data points given with (x 1 , y 1 ), (x 2 , y 2 ), … (x i , y i ) is used to find a function that describes dependence of y from independent variables x. A curve that is created from the experimentally measured data points can be presented in the generic power form given as: In equation (9), y* represent modeled function, while A and B are free fitting parameters. This form of the modeling function is the most suitable one since the parameter t n that should be relevant with modeled data also has power evolution. Parameters A and B specify the slope of the modeling function and determine the regression line. LSM defines the estimate of these parameters as the value which minimizes the sum of squares between the measurement (y) and the model (y*) which leads to expression: In the equation (10), ε stands for error which is a value to be minimized. The most suitable way to calculate parameters A and B is to introduce matrix notation in the following way:  After acquiring fitting parameters, A and B, it is needed to equalize equation (9) with the equation that describes capacitor charging (7). Even though it is mathematically simpler to use exponential form with least square approximation, negative exponent in the capacitor charging formula reverses the convexity of the function. Reversed convexity can introduce a lot of mathematical problems, delivering inadequate poor modeling results. After solving equations, variable τ is calculated. However, even after that calculation is concluded, problem of calculating precise resistance R and capacitance C that compose τ still remains. For mod-eling purpose, value of capacitance is set to 1 mF. Basic modeling circuit is given in Figure 2. Since the variety of the time constants can be found in the progress of the ∆V T , modeling with only one specific RC connection could lead to considerable dissent in the particular parts of the curve. Parameters A and B are calculated for the full period of stressing time.
The concept of improved modeling is to present the experimental curve as the product of multiple parts. These parts are to be determined by the phases of parameter n evolution, in favor of increasing of the accuracy of the model. To conduct this improvement, additions to the modeling circuit must be made. To follow evolution of ∆V T , charging rate of the capacitor is to be different in different time intervals. A method for enabling charging rate diversity is to increase the number of RC connections. With the range of different time constants, capacitor voltage can adapt more precisely to ∆V T . Adding switches and enabling only specific RC connections in determined time periods leads to capacitor voltage being additionally determined. Greater number of RC connections can be achieved in multiple ways, either by increasing number of resistors that are charging one capacitor, either by increasing number of capacitors that are charged through one resistor or either by increasing number of both resistors and capacitors. Since the principle of this modeling is that capacitor voltage models ∆V T , number of capacitors is set to one. This way, complicate procedure of summing of capacitor voltages that corresponds to ∆V T in different time periods is avoided, while concept is sustained. With only one capacitor, only way to increase number of RC connections is with increasing number of resistors.
With the goal to link phases of degradation with number of RC connections, modeling circuit is expanded with two additional resistors and switches. Expanded circuit is given in Figure 3.
In the starting moment, both switches, S 1 and S 2 are closed. Capacitor C is charged through all of three resistors. Since the resistors are connected in parallel, equivalent resistance is lower than the resistance of the resistor with lowest resistance. Growth of ∆V T is largest in starting phase of stressing, so the capacitor C is charged through lowest resistance, which is in line with model expectations. After starting phase, switch S 1 opens, so that capacitor C continues to charge through parallel resistance of resistors R 2 and R 3 only.
To provide suitable results and follow properly ∆V T , it is important that R 3 > R 2 > R 1 . That way, equivalent resistance of resistors R 2 and R 3 is greater than the equivalent resistance of all three resistors. After opening switch S 1 , capacitor is charged through higher resistance than before opening the switch, leading to slower capacitor charging and lower slope of the charging curve.
Higher charging resistance leads to even slower charging, which is, once again, in line with the model expectations, and describes development of ∆V T during stressing. To find appropriate values of resistance LSM is used again. This time, every phase is approximated separately, LSM is applied 12 more times, leading to parameters A 1 , A 2 , A 3 and B 1 , B 2 and B 3 . According to these parameters and in consideration with equivalent resistance equations, resistance of resistors R 1 , R 2 and R 3 are calculated. Results of different stressing conditions are given in the Table 3. Error of the estimated parameter is calculated using R-squared method and is between 0.95043 and 0.98268 for modeling with basic circuit and between 0.98764 and 0.99262 for modeling with expanded circuit. Modeled results using expanded modeling circuit are given in Figure 4 and 5. Results given in the Figure 4 and Figure 5 show that modeling error is considerably reduced if the modeling is done with taking in mind phase division of the pa-  rameter n evolution. This is notable especially during the second phase of ∆V T development. Since the third phase occurs after approximately 10 hours of stressing [6,38], during the experiment, samples are subdued to this phase the longest. When approximating full range of experimental data, LSM adapts parameters better for the longest lasting phase, and thereby, creating a slightly greater mismatch for the rest of the data. With this type of phase division, modeling error is more than twice decreased, as can be seen in Figure 6. Results given in Figure 6 show that modeling error has very similar form to the evolution of the parameter n, given in the Figure 1 (inserted graphic), which is fundamental signature of NBTI. With additional resistors and switches, used in the expanded modeling circuit, slope of the curve that shows modeling error is decreased (peak of the error is reduced by half, and error is reduced even more in the other parts of the curve). This result is in line with the modeling approach. Further expansion of the circuit (adding more RC connections, and splitting degradation phases into more phases that are shorter, would impact in even more decreased slope).
However, even with further increase in number of resistors, in the specific shorter time interval, time constant will be uniform. To design even more accurate model, number of RC connections is to depend not only on stressing time or phase, but on the actual numbers of individual defects in the circuit itself [29]. With increasing the number of defects, number of RC connections rises, thereby increases the overall sum of the voltages that comprise ∆V T .

Conclusions
Impacts of static negative bias temperature stressing in p-channel power VDMOSFETs IRF9520 have been reported. An equivalent electrical circuit is designed in order to model the behavior of ∆V T . Mathematical method of least square approximation is described and conducted to determine and to acquire parameters for values of modeling circuit elements. Phase division of parameter n evolution during modeling leads to better overall results in terms of accuracy and precision, regardless of stressing conditions. Future steps of work include improving the model with adapting it to the modeling of pulsed stress as well, since p-channel power VDMOSFETs are widely used in switching circuits because of their good switching characteristics.

Acknowledgments
This work has been supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia and in part by the Serbian Academy of Science and Arts.

Conflict of Interest
The author declares no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results