Extraction of Preisach model parameters for fluorite-structure ferroelectrics and antiferroelectrics

Flourite-structure ferroelectrics (FEs) and antiferroelectrics (AFEs) such as HfO2 and its variants have gained copious attention from the semiconductor community, because they enable complementary metal-oxide-semiconductor (CMOS)-compatible platforms for high-density, high-performance non-volatile and volatile memory technologies. While many individual experiments have been conducted to characterize and understand fluorite-structure FEs and AFEs, there has been little effort to aggregate the information needed to benchmark and provide insights into their properties. We present a fast and robust modeling framework that automatically fits the Preisach model to the experimental polarization (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q_{FE}$$\end{document}QFE) versus electric field (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{FE}$$\end{document}EFE) hysteresis characterizations of fluorite-structure FEs. The modifications to the original Preisach model allow the double hysteresis loops in fluorite-structure antiferroelectrics to be captured as well. By fitting the measured data reported in the literature, we observe that ferroelectric polarization and dielectric constant decrease as the coercive field rises in general.

Recent rapid development in data-intensive computing necessitates novel memory technologies beyond traditional flash memory. Ferroelectric field-effect transistors (FeFETs), as emerging memory, find a niche in such applications due to their ultra-fast program/erase time, low operation voltage, and low power consumption [1][2][3] . Despite the fact that hafnium oxide 4,5 and its doped variants (Al-doped 6,7 , Gd-doped 8 , La-doped 9 , Si-doped 1,10-12 , Sr-doped 13 , Y-doped 14,15 , Zr-doped 4,10,16-26 ) have been extensively studied and characterized over the past few years, little has been done to aggregate those data into ferroelectric properties to provide the insight necessary to create a predictive model for ferroelectrics. Such a predictive model cannot be realized without the accurate determination of a multitude of ferroelectric parameters from various experimental hysteresis loops ( Q FE -E FE ). However, direct determination of ferroelectric parameters, such as polarization ( P s and P r ) and coercive field ( E c ), from ferroelectric hysteresis characteristics can be inaccurate. Conventionally, the remnant polarization P r is evaluated at E FE = 0 MV/cm, the coercive field E c is determined at P FE = 0 μC/cm 2 , and the saturated polarization P s is the maximum P FE . Due to the contribution from the linear dielectric component of ferroelectrics ( ε 0 ε FE E FE ), evaluating the coercive field E c at Q FE = 0 μC/cm 2 directly would underestimate the actual coercive field and the value of P s is not obvious on the Q FE -E FE hysteresis loop (Fig. 1a). In addition, non-ideal effects such as asymmetrical E c± and polarization offset P offset hinder the direct parameter identification (Fig. 1b). Furthermore, the actual hysteresis loop measured could be a non-saturated minor loop, thus extracting ferroelectric parameters from such a loop would lead to significant underestimation (Fig. 1c). In summary, direct determination of ferroelectric parameters from a single hysteresis loop, a technique ubiquitously used in the literature, is prone to undervaluing ferroelectric parameters and unable to accurately extract some ferroelectric parameters.
To overcome the problem, we present a fast and robust modeling framework that automatically fits experimental ferroelectric hysteresis loops based on a computationally efficient Preisach model that had been widely adopted to model hysteresis of ferroelectrics in the past 27 , and we further extend this model to describe antiferroelectrics. We benchmarked our modeling framework and extracted ferroelectric model parameters from experimental ferroelectric hysteresis reported in the literature, demonstrating great agreement between the model and the measurement. We further observed that ferroelectric polarization and dielectric constant decrease with increasing coercive field in general. Such an observation can provide insight for developing predictive models for ferroelectrics and optimizing ferroelectric memory.

Results
Preisach model. To accurately extract experimental ferroelectric hysteresis, the original computationally efficient Preisach model 27 is modified to account for asymmetries in hysteresis loops. In the static Preisach model, the saturated hysteresis loop P sat -E FE is described by Eq. (1), where P s is the saturation polarization, P r is the remnant polarization, s is the slope parameter of the P FE -E FE hysteresis loop, E FE is the electric field in the ferroelectrics, P offset accounts for the shift of hysteresis loop center along the polarization axis, and E c± are the forward ( �E FE > 0 ) and backward ( �E FE < 0 ) sweep coercive fields of the ferroelectric that dictate the shift of hysteresis loop center along the electric field axis, respectively. To determine the non-saturated minor hysteresis loops, a linearly scaled version of P sat is determined by Eq. (2) where m is the proportionality factor and b is the offset polarization that both depend on previous hysteresis loops. Given that the turning points determined by previous hysteresis loops are ( V + , P + ) and ( V − ,P − ), the coefficients m and b can be determined by Eq. (3) The turning points ( V + , P + ) and ( V − , P − ) are the smallest possible endpoints from previous hysteresis loops that contain current hysteresis loop. If the current hysteresis loop is larger than all previous hysteresis loops, the turning points are the endpoints of the saturated major loop ( V s , P s ) and ( −V s , −P s ), where V s is much larger than t FE E c . Then the ferroelectric Q FE -E FE relation and the saturated Q FE -E FE can be obtained by summing the polarization contribution and the linear dielectric response (Eq. 4), , where E c,ext directly determined at Q FE = 0 μC/cm 2 is smaller than the actually E c from the P FE -E FE loops. In addition, P s is not obvious from the Q FE -E FE curve. (b) A hysteresis loop with asymmetrical E c± and non-zero P offset , where 2 P r,ext is smaller than 2 P r , leading to underestimation. The green square indicates the endpoint of the 2 P r when the start point (top yellow square) of the 2 P r aligned with that (top yellow circle) of the 2 P r,ext . (c) A minor hysteresis loop (blue) and the corresponding saturated hysteresis loop (dashed red), where 2 P r,minor is smaller than 2 P r and E c,minor is less than E c,ext . Extracting ferroelectric parameters from a minor loop leads to underestimation. where ε 0 is the permittivity of free space and ε FE is the relative linear dielectric constant of the ferroelectrics. Simulated ferroelectric hysteresis for both saturated and non-saturated hysteresis loops are shown in Fig. 2. The parameter extraction procedure is shown in Fig. 3, where the interior-point algorithm 28 is used to minimize the total squared error between experimental data and simulated results. Since the interior point algorithm is a local solver, a set of good initial model parameters can significantly reduce the computation time which is estimated from the measured hysteresis loops. The initial set of model parameters are generated automatically by Eq. (5).

Stop Conditions
Achieved?

Yes
Optimal  Fig. 4a, indicating good agreement between the model and experimental data. If only minor hysteresis loops are utilized to determine ferroelectric parameters, the extracted parameters will vary somewhat depending on the number of minor loops used. In the cases shown in Fig. 4b,c, where three and two minor loops are used for parameter extraction, respectively, the variations for P s , P r , and E c ± are within 6%, the variation for P offset is within ±1.05 μC/cm 2 , and the variation for ε FE is within 13 % comparing to the case where five loops are used (Fig. 4a). Regardless of the number of minor loops used, a saturated hysteresis loop that is much larger than the minor loops is reconstructed (Fig. 4b,c). This can drastically reduce the potential underestimation of ferroelectric parameters from measured minor loops.
Antiferroelectrics model. The aforementioned Preisach model was modified to capture the double hysteresis in antiferroelectrics through shifting. Antiferroelectricity is an electric field induce phase transition between a non-polar and a polar phase while ferroelectricity is an electric field-driven rotation of the polarization vector without a change in the crystalline symmetry. The microscopic mechanism behind these two phenomena is different, albeit related. At a mesoscopic level, both polarization rotation in ferroelectrics and the field-induced phase transition in antiferroelectrics are mediated by domain nucleation and propagation. On the other hand, the Preisach model is agnostic of the microscopic mechansim and captures only the domain dynamics occurring at the mesoscopic scale. As such, we believe that the modified Preisach model, within its limitations, captures the behavior of antiferroeletric materials as well. The saturated double hysteresis loops P sat -E FE is described by Eq. (6), where s is the slope parameter of the P FE -E AFE double hysteresis loop, E AFE is the electric field in the antiferroelectric, Dir is dictated by the sweep direction of the antiferroelectric electric fields, and sgn is the signum    (7) where m is the proportionality factor that depends on the endpoint of the former hysteresis loop. Since antiferroelectrics are volatile, for simplicity, only one endpoint ( V end , P end ) is considered. The coefficient m can be calculated by Eq. (8).
Similarly, the corresponding antiferroelectric Q AFE -E AFE relationship is given by Eq. (9), where ε 0 is the permittivity of free space and ε AFE is the relative linear dielectric constant of the antiferroelectrics. The modified Preisach model was implemented to identify model parameters based on experimental data using the aforementioned optimization technique. Figure 5 shows excellent agreement between the experimental data from a 10 nm thick antiferroelectric ZrO 2 sample and simulated results for both saturated major and unsaturated minor double hysteresis loops.

Discussion
The proposed modeling framework was used to extract ferroelectric parameters for Al-doped 6,7 , Gd-doped 8 , La-doped 9 , Si-doped 1,10-12 , Sr-doped 13 , Y-doped 14,15 , Zr-doped 4,10,16-26 , and undoped 4,5 HfO 2 thin films reported in the literature as shown in Fig. 6. By mapping the ferroelectric model parameters in parameter space, we ruled out unrealistic and far-fetched ferroelectric parameter combinations. We observed that as the coercive field E c increases, polarization ( P s and P r ) and the dielectric constant ε FE decrease (Fig. 6a-c). The coercive field E c and dielectric constant ε FE show lower bounds at ∼ 1 MV/cm, and 13, respectively. Saturated polarization has an upper bound around 30 μC/cm 2 . Such systematic relation between coercive field E c , polarization ( P s and P r ), and the dielectric constant ε FE is interesting and suggests that there are natural limits to the material, and understanding the microscopic origin of such interdependence will require further research. To achieve large memory windows and low operation voltage in ferroelectric memory, ferroelectric materials with low dielectric constant, high polarization, and small polarization difference ( P s − P r ) are preferred. The coercive field typically has an optimal value between 1 MV/cm and 2 MV/cm depending on other parameters because high E c hinders ferroelectric switching. Such a parametric study is crucial for ferroelectric memory design to achieve large memory windows and low program voltage under material restrictions as well as develop predictive models for ferroelectrics. Figure 7 demonstrated some fitted hysteresis loops that are used in Fig. 6. The model achieved an excellent agreement with experimental hysteresis and was able to reconstruct the saturated loops based on identified ferroelectric parameters if a non-saturated loop was measured.

Conclusion
We developed a fast and robust modeling framework for the automated determination of ferroelectric and antiferroelectric parameters from experimental data by modifying the Preisach model. We demonstrated excellent consistency between modeled results and measured data reported from the literature. By plotting the extracted ferroelectric model parameters in parameter space, we observed that ferroelectric polarization and dielectric constant tends to decrease with the increasing coercive field. This parametric study suggests that it is important to maximize the ferroelectric polarization and minimize the dielectric constant and polarization difference ( P s -P r ) to achieve a large memory window and low voltage of operation of ferroelectric memory.

Methods
Sample preparation. Both the HZO (metal-ferroelectric-metal) and ZrO 2 (metal-antiferroelectric-metal) capacitor structures were fabricated on p+ Si (100) substrates with a native SiO 2 layer. For the HZO sample, the bottom TiN electrode (12 nm), 10 nm Hf 0.5 Zr 0.5 O 2 film, and the top TiN electrode (12 nm) were deposited subsequently using plasma-enhanced atomic layer deposition (PEALD) technique in a Veeco Fiji G2 Plasma ALD system. The deposition was carried out at 250 • C using Tetrakis(dimethylamido) hafnium, zirconium, and titanium precursors with water as the oxygen source. A post-TiN metallization annealing was done on the sample at 450 • C for 30 seconds in nitrogen atmosphere to crystallize the Hf 0.5 Zr 0.5 O 2 layer. For the ZrO 2 sample, the bottom TiN layers (10 nm) and 10 nm ZrO 2 film were deposited in a 300 mm TEL Trias TM clean-room tool at a temperatures of ∼430 • C and 350 • C, respectively. The temperature is high enough to crystallize ZrO 2 during growth such that no additional annealing was required to achieve antiferroelectricity. Afterward, the top TiN layer of this sample was deposited using a Plasma Enhanced ALD (PEALD) system at Georgia Tech. After ALD and crystallization of both the samples, an Al metal layer (100 nm) patterned into rectangular shapes were evaporated to define the capacitor with areas of (200 μm) 2 , (100 μm) 2 , (50 μm) 2