Analytical Procedure for the Extraction of Material Parameters in Antiferroelectric ZrO2

Here, we present an analytical procedure to extract the anisotropy constants of antiferroelectric (AFE) materials from a few key features of the experimental polarization versus field curves. Our approach is validated for two experimental datasets of ZrO2 capacitors, and the extracted parameters are consistent with the microscopically nonpolar nature of the zero-field state of the AFE ZrO2. The methodology has applications in AFE nonvolatile memories and memristors, as well as in electron devices exploiting the negative capacitance (NC) operation of ZrO2.

Preisach's model [14]. Recently, the AFE behavior has been also observed in hafnium-and zirconium-based materials [15], which exhibit also ferroelectricity and are of great interest due to their scalability and CMOS process compatibility. The microscopic picture behind antiferroelectricity in ZrO 2 is fundamentally different compared with PbZrO 3 and similar perovskites. In fact, ab initio calculations have revealed that the energy ground state of thin ZrO 2 films is tetragonal [16], which has been also confirmed by GIXRD measurements [5], [15], [17], so that at zero applied field, the material is microscopically nonpolar [16], [18]. By applying an electric field to the ZrO 2 , a phase transition is induced from the nonpolar tetragonal phase to a polar orthorombic phase, which is the phase also responsible for ferroelectricity in hafnium-zirconium oxides (HZOs).
While Kittel's model gives an adequate description for AFE materials having an anti-polar alignment of the domains, it may not be suitable to describe the physical picture governing the antiferroelectricity in ZrO 2 , which, as stated before, is quite different from the one observed in perovskites. In this article, we propose a procedure to extract the material parameters of the AFE ZrO 2 in the framework of the multidomain Landau, Ginzburg, Devonshire (LGD) model, that can be applied to AFE materials with microscopically nonpolar ground state, such as ZrO 2 . The calibrated LGD model can reproduce fairly well both the quasi-static polarization-field curves in [9] and [17], and the transient negative capacitance (NC) behavior reported in [17]. Moreover, the parameterization of ZrO 2 is consistent with its microscopically nonpolar state at zero applied field. This article is organized as follows. In Section II, we propose a methodology to extract the anisotropy constants of the LGD model for an AFE material with nonpolar ground state. In Section III, we provide a quick overview of the simulation framework used to validate the proposed extraction procedure. In Section IV, we show comparisons between simulations and experiments for different ZrO 2 thicknesses and operation regimes. In Section V, we offer a few concluding remarks.

II. EXTRACTION OF ANISOTROPY CONSTANTS
Here, let us consider a capacitor with metal electrodes and an FE or AFE dielectric (DE). For a simple homogeneous polarization picture, the Gibbs' free energy of the system consisting of the capacitor and the external battery can be  Table I, and for either a zero or a positive applied field E F . At zero field (blue curve), the point E F = P = 0 is a free energy minimum, and thus, it is a stable steady-state point for the system. The application of a positive E F (purple curve) shifts the energy minimum to a positive P. (b) Measured total polarization versus applied electric field in a TiN/ZrO 2 (9.5 nm)/TiN stack [9]. The meaning of points B and C is discussed in the text. written as follows [19]: where ε 0 is the vacuum permittivity, and P, ε F , and E F are the spontaneous polarization, background permittivity, and electric field of the FE or AFE material, while α, β, and γ are the anisotropy constants. The quasi-static P-E F trajectories are identified by the conditions (dG/d P) = 0 and (d 2 G/d P 2 ) > 0, [20], namely Quasi-static experiments in a metal-FE-metal or metal-AFE-metal (M-AFE-M) stack probe the overall charge in the system, usually denoted as total polarization P T ≈ Q = P + ε 0 ε F E F . Fig. 1(a) shows an example of the free energy landscape for an M-AFE-M system, and (2) prescribes that α be positive in order to have a microscopically nonpolar stable state at E F ≈ 0 and P ≈ 0. Fig. 1(b) displays the experimental P T versus E F curve recently reported for a ZrO 2 capacitor [9]. In Fig. 1(b), we denote by E B and E C the coercive fields corresponding, respectively, to the nonpolar to positive and positive to nonpolar transition in the P T -E F curve. In practice, the points (E B , P T,B ) and (E C , P T,C ) can be identified as the points where the P T versus E F curve exhibits a clear change in the slope. In order to define an analytical procedure for the extraction of material parameters in AFE ZrO 2 , we now assume that points B and C correspond, respectively, to a maximum and a minimum of the static E F − P T relation implied by the LGD polynomial. In the Appendix, we show that such maximum and minimum of the E F − P T relation coincide with those of the E F − P relation, which, in turn, are readily identified by the condition (∂ E F /∂ P) = 0 in (2b). Hence, the conditions ensuring that the quasi-static P T − E F trajectories include points B and C become From (3a) to (3c), we can readily express α, β, and γ as follows: Equation (4) provides the anisotropy constants in terms of E C , P C , and P B . However, the spontaneous polarizations P C and P B cannot be directly identified in the experimental curves of Fig. 1(b), but they must be calculated by using P = P T − ε F ε 0 E F . This implies that α, β, and γ in (4) are given in terms of E C , P T,C , and P T,B and of the remaining parameter ε F . In this latter respect, it has been theoretically argued that ε F should be considered an adjustable parameter rather a true material constant [21], and in practice, it is difficult to extract ε F independently of α, β, and γ. Therefore, we now substitute α, β, and γ from (4) into (3d) and rearrange it as follows: By can now be solved for ε F . Namely, ε F can be used as the fourth adjustable parameter determined by (3), so as to ensure that the quasi-static P T − E F trajectories include the points B and C in Fig. 1 As it can be seen, (5) implies also (5P 2 B − P 2 C ) > 0 (because 5P 2 C is by definition larger than P 2 B ), which, in turn, results in positive α and γ values and in a negative β value [see (4)]. As already mentioned, the positive α value is consistent with the microscopically nonpolar nature of thin ZrO 2 films at a zero applied field, and it is also consistent with previous  (6), and describing the domain wall energy contribution. d and w denote, respectively, the domain size and the width of the domain wall region [19].
literature for AFE ZrO 2 [22], [23]. Moreover, a positive α value is the only possible choice in order to obtain no remnant polarization.
In summary, the procedure to extract the anisotropy constants from experiments requires to first identify the points B and C in the measured P T -E F curves [see Fig. 1(b)]. Then, (5) can be solved numerically to determine ε F , and once ε F is known, (4) provides expressions for α, β, and γ.

III. FRAMEWORK FOR NUMERICAL MODELING
In Section IV, we will illustrate several comparisons between simulations and experiments aimed at a validation of the extraction procedure for the anisotropy constants of ZrO 2 . All simulations were carried out by using the solver for the multidomain LGD equations that has been already discussed in [19], [24], and [25]. In this section, we recall only a few aspects of the simulation framework, which are relevant for the cases at study in this article. For an AFE or an FE material consisting of n D domains, as shown in Fig. 2, the dynamics of the polarization P i in each domain is described by the following equations [19], [24], [25]: where k is the domain wall coupling coefficient, ρ is the switching resistivity, while 1/C (dep) i, j = 1/2(1/C j,i + 1/C i, j ), and the terms C i, j are the capacitive couplings between domains. Given the similarity between the crystal chemistry of ZrO 2 and HfO 2 [26], in simulations, we used k ≈ 0, as suggested by recent first principle calculations for HfO 2 [27]. For each domain, the α i , β i , and γ i values were calculated by using a Gaussian distribution of the coercive fields with the mean E C and E B values used to extract the parameters in Table I and with a ratio σ EC = σ EB between the standard deviation and mean value; ε F is the same in all domains. All simulations were performed using n D = 400 domains i, j , which, as stated before, describe the capacitive coupling between the ith and jth domains. However, such a capacitive coupling decreases quite steeply with the distance between domains, so that simulations become insensitive to n D for large enough n D values. Moreover, for an M-AFE-M stack (without the DE layer in Fig. 2), the terms 1/C (dep) i, j tend to zero, because there is no electrostatic coupling between the domains through the DE layer, which further reduces the sensitivity to n D of the simulations results. The experimental P T versus E F curves for AFE ZrO 2 sometimes exhibit a non-negligible polarization at zero field, that is ascribed to the presence of FE domains. Therefore, in our simulations, we accounted for a small fraction of FE domains, which can be included in our model by setting appropriate values of the anisotropy constants for such domains. More precisely, for the LGD parameterization of FE domains in ZrO 2 , we used educated guesses for the remnant polarization P r ≃ 25 µC/cm 2 and coercive field E C,FE ≃ 1.2 MV/cm consistent with [28] and [29], resulting in the following LGD parameterization: α FE = −5.94 · 10 8 m 3 /F, β FE = 4.28 · 10 9 m 5 /(FC 2 ), and γ FE = 1.16 · 10 9 m 9 /(FC 4 ). Even for FE domains, we introduced a Gaussian distribution of the coercive field, with the same σ EC value used for AFE domains. The domain-dependent anisotropy constants have a spatially random distribution across the domain grid, and we have verified that their spatial distribution does not practically affect the simulation results. This is not unexpected especially for M-AFE-M stacks, where there is no electrostatic coupling between the domains. In our simulations, it is also possible to include a small built-in electric field in the FE material, possibly arising from a slight work-function difference at the two electrodes or from fixed charges in the DE stack. Table I reports the material parameters extracted with the methodology of this work from two experimental datasets, namely, the P T -E F curves recently reported in [9] and [17]. Quite interestingly, from the parameters in Table I, one can calculate the zero-field permittivity of ZrO 2 , which is defined as 1/ε 0 (∂ P T /∂ E F ) at E F = P T = 0. By recalling P T = P + ε 0 ε F E F and using (2a) for (∂ P/∂ E F ), the zero-field permittivity is readily expressed as (ε F +1/(2αε 0 )). As already mentioned in Section II, while ε F is related to the zerofield permittivity, which is the quantity actually measured in experiments, it is not equivalent to it and can be thought as a fitting parameter [21]. The zero-field permittivity obtained from the parameters in Table I ranges between 30 and 40, which is in good agreement with experimental values in [30] and [31]. The anisotropy constants in Table I provide the mean values of the domain-dependent α i , β i , and γ i parameters used in the numerical simulations.

IV. COMPARISON WITH EXPERIMENTAL RESULTS
In Fig. 3, we show a comparison between simulations and experiments for the P T -E F curves of the M-AFE-M stacks reported in [9] and for the materials parameters in Table I. In Fig. 3(a), we considered 3% of the overall domains to be FE with the parameters discussed in Section III, while for Fig. 3(b), we did not include FE domains, as it can be seen that the hysteresis of the P T -E F curve is completely closed at zero applied electric field.
In our model, the timescale for the polarization dynamics is t ρ = ρ/(2⟨α⟩) [24]. For ρ ≈ 400 m (see Table I), we have t ρ ≈ 70 ns, which is consistent with the literature for large-area devices [32] and ensures that simulations in Fig. 3 are quasi-static. While it could be argued that each stack could have its own ρ value, there is no direct measurement to extract it, rather it is usually inferred from polarization switching measurements [33]. Given the lack of a direct information about the value of ρ, we kept its value fixed for all stacks. The agreement between simulations and experiments is fairly good for both t ZrO 2 values because of a good symmetry of the experiments along both the P T and E F axes. Fig. 4 reports a similar comparison for the experimental dataset in [17], where we considered 4% of the domains to be FE. The agreement between simulations and experiments is still fairly good, but we also observe a discrepancy in the negative E F hysteresis. This is mainly due to an asymmetry in the measured P T values for positive and negative E F at large |E F |, possibly due to a non-negligible influence of leakage. In fact, while an asymmetry along the E F axis can be included in our model through a built-in field E B I (see Table I), the LGD model is instead inherently symmetric in the P T values.
Hoffmann et al. [17] also reported transient NC experiments, that we here analyze by using the LGD model, as previously reported for the NC behavior in FE devices [24], [34], [35], [36]. In the TiN/ZrO 2 /Al 2 O 3 /HfO 2 /TiN stack, the undoped HfO 2 layer is paraelectric, and the thicknesses are as follows: t ZrO 2 = 10 nm, t Al 2 0 3 ≈ 1 nm, and t HfO 2 = 8 nm. The timescale of the voltage pulses in these experiments (now comparable to the t ρ ), and the relatively thick DE were on purposely chosen to minimize the role of charge injection and trapping [17]. Therefore, our simulations neglect trapping, which has been shown to be instead important in quasi-static measurements for FE-oxide stacks having a thin DE layer [25], [37]. Fig. 5(a) and (c) compares the simulated and experimental P T -V max curves for a pulse width of, respectively, 275 ns and 1 µs, where V max is the amplitude of the voltage pulse, and the simulated P T values have been extracted following the definition in [17]. Fig. 5(b) and (d) displays the corresponding plots for the P T versus an effective electric field, E EFF , across the ZrO 2 layer. In experiments, the E EFF cannot be directly probed; hence, it was estimated as E EFF ≈ (V max − P T /C D )/t ZrO 2 , where C D ≈ 1.78 µF/cm 2 is the effective Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply. capacitance of the Al 2 O 3 -HfO 2 series [17]. In the simulations of Fig. 5(b) and (d), the E EFF was calculated according to the same definition given in [17]. Fig. 5 shows that the same ZrO 2 parameters already employed in Fig. 4, both LGD mean values and their statistical distribution, can provide a fairly good agreement also for transient NC experiments, with a matching between simulations and experiments that is similarly good for the two different pulse widths. The results in Fig. 5 reinforce our confidence in the extraction procedure for the ZrO 2 parameters and in the overall simulation framework.

V. DISCUSSION AND CONCLUSION
We have proposed a procedure to extract the material parameters for the LGD model of AFE ZrO 2 films, which is consistent with the microscopically nonpolar nature of the zero-field state in AFE ZrO 2 [17]. The points (E B , P T,B ) and (E C , P T,C ) necessary to extract the anisotropy constants can be reliably identified by a distinct change in P T versus E F slope of the P T -E F curves [see Fig. 1(b)], provided that the curves are not significantly distorted by leakage and that they display a full hysteresis loop, as opposed to minor loops. Our methodology was successfully validated by considering quasi-static P T -E F curves in M-AFE-M stacks, where a small fraction of FE domains was also included in the model to explain and reproduce the residual polarization at zero field observed in some AFE ZrO 2 films [9]. Moreover, we analyzed very recent experiments reporting both P T -E F curves in M-AFE-M stacks and transient NC experiments in a TiN/ZrO 2 /Al 2 O 3 /HfO 2 /TiN stack. With a single set of ZrO 2 parameters extracted from the P T -E F curves, our simulations could reproduce fairly well also the transient NC experiments and for different pulse widths.
Our analytical extraction procedure has a clear physical background, and it is easy to implement, although it may lead to fitting results that are not as accurate as those obtained with more phenomenological approaches [13], [14], where the quite many parameters of the models are typically extracted by using numerical procedures targeting a minimization of the mean-squared error between simulations and experiments.
We believe that the methodology proposed in this article to extract the anisotropy constants of antiferroelectic ZrO 2 layers will have useful applications in FE nonvolatile memories and memristors, as well as in possible devices exploiting the ZrO 2 NC behavior.

APPENDICES
In order to show that the maximum and minimum of the static E F -P T curve coincide with those of the E F -P curve, we can substitute P = P T − ε 0 ε F E F in (2a) and obtain (7) Then, we derive both sides of (7) with respect to P T , and we have which can be refactored in ∂ E F ∂ P T 1 + ε 0 ε F (2α + 12β P 2 + 30γ P 4 ) = 2α + 12β P 2 5 + 30γ P 4 . (9) Equation (9) clearly shows that the condition ∂ E F ∂ P = 2α + 12β P 2 + 30γ P 4 = 0 implies also (∂ E F /∂ P T ) = 0.