An Extended Kolmogorov–Avrami–Ishibashi (EKAI) Model to Simulate Dynamic Characteristics of Polycrystalline-Ferroelectric-Gate Field-Effect Transistors

A physics-based model on polarization switching in ferroelectric polycrystalline films is proposed. The calculation results by the model agree well with experimental results regarding dynamic operations of ferroelectric-gate field-effect transistors (FeFETs). In the model, an angle θ for each grain in the ferroelectric polycrystal is defined, where θ is the angle between the spontaneous polarization and the film normal direction. Under a constant electric field for a single-crystal film with θ = 0, phenomena regarding polarization domain nucleation and wall propagation are well described by the Kolmogorov–Avrami–Ishibashi theory. Since the electric fields are time-dependent in FeFET operations and the θ values are distributed in the polycrystalline film, the model in this paper forms an extended Kolmogorov–Avrami–Ishibashi (EKAI) model. Under a low electric field, the nucleation and domain propagation proceed according to thermally activated processes, meaning that switching the time scale of a grain with the angle θ is proportional to an exponential form as exp(const./Ezcosθ) [Ez: the film-normal electric field]. Wide θ distribution makes the time response quite broad even on the logarithmic scale, which relates well with the broad switching time experimentally shown by FeFETs. The EKAI model is physics based and need not assume non-physical distribution functions in it.


Introduction 1.Necessity of a Physics-Based Model for Ferroelectric Device Dynamics
Ferroelectric-gate field-effect transistors (FeFETs) have attracted attention not only because of their potential functionality [1,2] such as being compact, non-volatile, and nondestructive-read memory cells, but also due to recent rich accumulations of experimental data.FeFETs have been characterized by various measurements.Quasi-static characteristics of drain current (I d ) vs. gate voltage (V g ) are measured with slow V g sweeping by a semiconductor parameter analyzer.Also, Q m vs. V g are measured by a quasi-static procedure with the slow V g sweeping by a ferroelectric tester.The Q m is the metal-gate charge density induced by ferroelectric layer polarization.Dynamic time-dependent properties of FeFETs are characterized by a pulse-write-and-V th -read (PWVR; see Figure A3) or pulse-write-and-I d -read (PWIR) method.The pulse-write (PW) is done by a pulse generator.In a FeFET, while a positive or negative V g pulse is given with the V g height of V h (>0) or V l (<0), ferroelectric polarization responds dynamically as a function of pulse height (V h or −V l ) and time width (t w ).After ceasing PW, considerable amounts of ferroelectric polarization are retained due to ferroelectricity.This retained polarization is a function of V h , V l , and t w .A read operation is done after PW, where V g of the transistor is swept in a narrow voltage range slowly compared to the PW time scale.
A threshold voltage (V th ) is read for PWVR, and the drain current (I d ) is read at a fixed V g for PWIR.A threshold voltage difference (∆V th ) for PWVR or a current ratio between the highand low-I d -current states for PWIR can be derived for a pair of the write voltages; V h and V l .∆V th by PWVR or the I d ratio by PWIR is a performance indicator of an FeFET as a nonvolatile memory transistor.
A lot of experimental results by PWVR and PWIR are available to investigate dynamic polarization switching in metal ferroelectric insulator-semiconductor (MFIS) type FeFETs, whose ferroelectric is SrBi 2 Ta 2 O 9 or (Sr 1−x Ca x )Bi 2 Ta 2 O 9 [3][4][5].The former and the latter are abbreviated as SBT and CSBT, respectively.The SBT and CSBT belong to a family of the Bilayered perovskite oxides with similar electrical properties.The coercive field of CSBT for x = 0.1-0.2 is 10% larger than that of SBT.Hereafter, FeFETs comprising SBT or CSBT in the ferroelectric layers are called SBT-FeFETs.The electric properties of SBT-FeFETs are suitable for model verification because they are reproducible with high switching endurance due to negligibly small charge injection and trapping in the FeFET gate stacks.
Despite the many experimental works of PWVR and PWIR reported in detail, there have been no correct theories based on physics for simulating them.As one of the existing models, a phenomenological Landau-Ginzburg-Devonshire (LGD) theory in ferroelectrics constructs a Gibbs free energy for the ferroelectrics where the primary energy term is expressed as a polynomial expansion form of polarization as an order parameter [6][7][8].By minimizing the free energy, equilibrium or quasi-static properties are derived.The LGD theory described the transition between the paraelectric and ferroelectric phases with temperature variation across the Curie temperature [7].It also described phase transitions between the tetragonal and orthorhombic and between the orthorhombic and rhombohedral in BaTiO 3 , a cubic-based ferroelectric material [6].In the Gibbs free energy, the spatial differentiation term brought non-homogeneity, by which domain patterns were described [9,10].The Gibbs energy for ferroelectrics can include the elastic energy term and the coupling terms of elasticity and polarization since the polarization accompanies lattice displacement of crystals.In the case that the ferroelectric layers are divided into LGD segments, the LGD theories are called phase field models [11,12].The models showed coexistence of 180 • -and non-180 • domains under constraints of substrates [11] and polycrystalline grain growth [13].
The LGD theory has its own representation of time dependence.In the case that the electric properties without elastic terms are matters of concern, the description of time dependence is sometimes called the Laudau-Khalatnikov (LK) equation [14][15][16].In the phasefield models, it is called the time-dependent Ginzburg-Landau (TDGL) equation [11,17].LK and TDGL are essentially the same.A Lagrangian function, which includes a simple viscosity term of classical mechanics [18], may help us comprehensively understand the origin of the time-dependent polarization switching.The Gibbs free energy for the ferroelectrics replaces the potential term in the Lagrangian.In deriving equations of motion, the inertia term is usually neglected for polarization reversal problems, and the LK and TDGL equations are obtained [19].Overdamped (i.e., strong viscosity) cases in the Lagrangian may allow us to omit the inertia term in the equation of motion.As is obvious from the derivation of these time-dependent equations, however, the LK and TDGL equations are deterministic and do not have thermal activation forms in them.As indicated by previous papers [20][21][22], nucleation and domain-wall propagation under an electric field occur via thermally activated steps.Thus, the LK and TDGL equations are not suitable (at least quantitatively).Additions of Gaussian white-noise terms to these equations are necessary for describing the thermal activation steps, and in fact the effect of the noise term was demonstrated for a small segment number case [23].However, since the LDG segments should be finer than the domain-wall width and the wall width size is one crystal lattice or a few times of it [24,25], the number of the LDG segments for general objects including polycrystalline films becomes so huge that we cannot make calculations within a realistic time scale.This paper proposes a physical model for describing polarization variation with time in a ferroelectric film and for calculating electronic device operations of FeFETs, metal-ferroelectric-metal (MFM) capacitors, and metal-ferroelectric-insulator-metal (MFIM) capacitors.In the model, the ferroelectric film can be a polycrystalline one, and to each grain, an angle θ is assigned where θ is the angle between the film normal and the direction along the spontaneous polarization.The θ distribution in the polycrystalline film is given by experiments.In fact, for SBT-based FeFETs, the angle distribution was derived by an electron backscattering diffraction (EBSD) patterns technique [26].
In each grain having θ, it is also assumed that, under an electric field, a seed for polarization reversal grows along the spontaneous direction and forms a narrow columnar 180 • reversed domain, which is expanded along the sidewise direction, because 180 • switching and ±180 • domain formation are subjected to occur in SBT and CSBT ferroelectrics [27][28][29].Experimental observations of the 180 • switching and 180 • domain-wall moving for θ ̸ = 0 cases are found in [30,31].
The Kolmogorov-Avrami-Ishibashi (KAI) model [32][33][34][35][36] fits the physical picture mentioned above (i.e., 180 • domain nucleation and domain wall propagation).The KAI model provides mathematical functions of the θ = 0 case describing domain expansion with time under a constant electric field.As for epitaxial thin films, many papers indicated good agreement of the KAI model with experiments [37,38], but for polycrystalline films, the KAI model is asserted to be unsuitable [39].This complicated issue concerning polycrystal films will be overcome in our present paper by introducing θ terms in grains.In order to explain wide-range log (t) distributions of polarization switching times that Pt/Pb(Zr,Ti)O 3 /Pt capacitors showed, the inhomogeneous field mechanism (IFM) model [40] and the nucleation-limited switching (NSL) model [41] are proposed.Although switching reversal experiments in polycrystalline films were explained [42,43] using the IFM, the model assumed inhomogeneous electric-field distribution of a Gaussian type, which was not physics based.The NLS is a non-KAI model.A (111)-oriented Pb(Zr,Ti)O 3 was used.The waiting times of elementary regions are stochastic.By giving an exponentially broad distribution of the waiting time, the switching phenomena across a wide log (t) scale were described.Since the distribution function is assumed so that the numerical results reproduce the experimental ones, understanding the physics of the distribution function is not easy.If we infer a possible inclusion of grains with other crystal orientation in the film, physics in the wide log (t) characteristics might be explained by an effect of the orientation distribution.A HfO 2 -based FeFET is reported to be consistent with the NLS [44].The paper indicates that the grain size is ≈20 nm and the field for creating a nucleus is ≈1 MV/cm.This suggests that the nucleation size is already the same as the grain size, and thus there is no space left for the created nucleus to induce a wall expansion supposed in the KAI model.
Herein we propose an extended KAI (EKAI) model, which is applicable to represent characteristics of FeFETs with polycrystalline ferroelectrics under time-dependent electric fields.In the EKAI model, KAI-like pictures are adopted only inside individual grains.The EKAI describes 180 • switching and ±180 • domain wall propagations in every grain separately.No wall motions are assumed to propagate across adjacent grains.There are two fundamental premises of the EKAI.One is that only the electric field component along the spontaneous polarization works for the polarization switching.The other is that we pick up only the film-normal component of the switched polarization because the film-normal component is the important response to the externally applied potential.The in-plane components of the switching polarization are expected to be randomly distributed among grains.The effect of the in-plane components is thought to be weak because of the random distribution, and thus the boundary conditions between neighboring grains are ignored in the present EKAI model.The polarization of the KAI model is an averaged quantity where the polarization is represented by the average of the volume fraction of ±180 • domain regions.The EKAI has a distinct perspective on ferroelectrics from the phase-field model in which the ferroelectrics are divided into small segments connected by strict boundary conditions according to the idea of the LGD.The EKAI has a much shorter computation time than the phase-field model because it ignores the grain boundary conditions and averages the ferroelectric polarizations in grains.The PWVR and PWIR data show a wide time range of polarization switching characteristics from 50 ns to 10 ms [3][4][5].We suppose that these rather slow and wide-ranged time responses are attributed to ferroelectric polycrystals consisting of grains that have broad distributions in the crystal orientations.The response times of the paraelectric components in the ferroelectric grains and the dielectrics in the insulator are supposed to be much shorter than those of the polarizationswitching components.The potential formation time in the semiconductor is also supposed to be shorter.Therefore, the EKAI model in this paper assumes that the parameters except for the polarization switching in the ferroelectric grains vary instantaneously.Correctness of the EKAI model is supported by good agreement with experimental results of quasi-static I d -V g , Q m -V g , and PWVR for SBT-based FeFETs as discussed later in Section 4. The EKAI model in this paper assumes no free charges existing in the gate insulator and ferroelectric; thus the EKAI is not directly applicable to HfO 2 -based FeFETs in which the ferroelectric polarization switching is always accompanied by charge injection currents [45][46][47][48][49].

Presentation of the KAI with the Characteristic Times for the EKAI Model
Let us write some equations of the original KAI model, because they are necessary in the succeeding section for the EKAI derivation.In the original KAI, polarization switching nucleation is instantaneous in comparison with domain growth motion.A constant electric field and the domain wall velocity (v wall ) depending on this constant field are assumed.As shown in Figure 1, the ferroelectric volume is constituted by domain regions.The switching polarizations in the downward and upward regions are P s and −P s , respectively.P s is the spontaneous polarization.The polarization direction is parallel to the z-axis.The downward (or upward) domain regions expand with time under a positive (or negative) field, after application of step function with a constant positive field, E z+ .Consider a case that at the initial (t < 0) the volume is occupied fully by the upward regions, and a step function with a constant positive field, E z+ is applied at t = 0.Then, the volume fraction of the downward domain (R dn∥z ) is varied with time for t ≥ 0 under a constant E z+ as with t o is the characteristic time for the polarization switching and is a function of the constant field E z+ .The power exponent, n, is a parameter relating to the domain growth dimension.The volume fraction of the upward domain (R up∥z ) is The switching polarization (P z ) averaged over the ferroelectric volume is Similarly, in the case that at the initial (t < 0) the volume is occupied fully by the downward regions, and step function with a constant negative field, E z− is applied at t = 0, the volume fraction of the upward domain (R up∥z ) is for t ≥ 0 under a constant E z− : R up∥z = c(t). ( t o in Equation ( 2) is a function of the constant field E z− .The switching polarization (P z ) averaged over the ferroelectric volume is with is the characteristic time for the polarization switching and is a function of the constant field  .The power exponent,  , is a parameter relating to the domain growth dimension.The volume fraction of the upward domain ( ∥ ) is According to Ishibashi and Takagi [36] and Ishibashi [50], in category 1, domain nucleation occurs with a fixed probability, and n = 3 when the wall shape is two-dimensional (i.e., circular) and n = 2 when the wall shape is one-dimensional (i.e., straight line).In category 2, there exist latent nuclei, and no new nucleation appears.In category 2, n = 2 when the shape is two-dimensional, and n = 1 when it is one-dimensional.In real materials, domain nucleation may occur, and latent nuclei may also exist.Some parts are twodimension and others are one-dimensional-like.Hence the value n is not an integer of the range 1 ≤ n ≤ 3. When t o is small in Equation ( 2), the transient time of the polarization reversal is short, and when t o is large, it is long.That is, t o is a characteristic time that gives a time scale of the polarization variation.t o has a relationship with v wall .
For category 1, For category 2, In Equation (7), γ 1 and γ 2 are constants.The KAI model did not show explicit mathematical forms of t o .The wall velocity v wall as a function of the applied film-normal electric field (E z ) has been investigated mainly by switching current measurement of single crystals [51,52] and by piezoresponse force microscopy (PFM) of epitaxial thin films including random disorder by defects [21,53,54].A consensus view in the case of films including the disorder is that v wall under high electric fields quickly increases in linear equation whereas it creeps up in reciprocal exponential at low fields [52,55].In the case of ferroelectric films including defects, v wall under high field is expressed as v wall ∝ E z for E z ≫ E dpin , where E dpin is a critical field over which depinning of wall motions occurs at 0 K [55,56].In an intermediate E z region, v wall is expressed using a power exponent τ ′ as and in a low-E z creep region, v wall is expressed at finite temperature as where U ′ is a scale of energy barrier and σ is a power exponent originated from random disorder defects in ferroelectric films [55,56] (k B : the Boltzmann constant and T: the absolute temperature).
From Equations ( 7) and ( 8), the v wall expressions are changed to t o expressions.Using a renormalized constant τ, Equation (8a) is, irrespective of the category 1 or 2, Similarly, irrespective of the categories, Equation (8b) is converted to t o , using a renormalized constant U, as where t in f is a constant.
In the present work, experimental results of SBT-based FeFETs are used in order to demonstrate the EKAI model's credibility.High endurance of the SBT FeFETs can be realized on small or at least moderate write-voltage conditions where the charge injection and trapping in the gate stack are suppressed [2,57].Since a small or moderate write voltage brings a low electric field in the ferroelectric, we adopt Equation (9b) rather than Equation (9a) hereafter in this paper; i.e., Equation (9b) will be used in calculations shown later.Further, since the separate determination of U and E dpin is difficult, we have the following equation for t o in this paper: with an activation field constant (10).

EKAI Model
Let us consider a ferroelectric polycrystal film as shown in Figure 2a.The thickness of the film is d f .The film is divided into plural grains labeled l from one to l max .All grains have the common thickness d f .Grain boundaries are along the film normal, i.e., the z-axis.The area of grain l is A l .The direction of the spontaneous polarization of grain l is parallel to the u l axis.The angle between the u l axis and the z-axis is θ l .The range of θ l is defined as 0 ≤ θ l ≤ 90 • .All grains have a same spontaneous polarization, P s or −P s .The polarization of the direction of the u l axis is defined as P s .
In Section 4, the EKAI model will be compared to the experimental results of MFIS FeFETs, where the ferroelectric layers consist of ferroelectric polycrystals.Experimental FeFETs with 135 nm thick ferroelectric SBT layers are available for discussion in the present EKAI model.Averaged in-plane diameters of the SBT grains are about 200 nm, which is larger than the ferroelectric layer thickness in the FeFETs.Therefore, we shall assume a single grain occupation along the z-direction or the film normal in the SBT FeFETs.Each grain stands as a pillar with a constant cross-section from the film top to the bottom.
Figure 2b is an expanded schematic picture of grain l.The grain consists of upward domain regions and downward domain regions.In the downward (or upward) domain regions, the spontaneous polarization P s (or −P s ) is parallel (or anti-parallel) to the u l axis.Figure 2c is a schematic picture focusing on a cylindrical-shape downward domain existing in the grain.The EKAI model in this section describes polarization variation in one grain.In the model, there are two fundamental assumptions (described in Sections 2.1 and 2.2) prior to the significant description (Section 2.3) of the polarization dynamics under a varying electric field with time.Section 3 provides calculation schemes for the following specific devices including ferroelectric polycrystal: (1) MFM capacitors, (2) MFIM capacitors, and (3) MFIS FeFETs.
Let us consider a ferroelectric polycrystal film as shown in Figure 2a.The thickness of the film is  .The film is divided into plural grains labeled  from one to  .All grains have the common thickness  .Grain boundaries are along the film normal, i.e., the z-axis.The area of grain  is  .The direction of the spontaneous polarization of grain  is parallel to the  axis.The angle between the  axis and the z-axis is  .The range of  is defined as 0  90°.All grains have a same spontaneous polarization,  or − .The polarization of the direction of the  axis is defined as  .
(c)  l has an orientation angle θ l and area A l .θ l is the angle between the z-axis and the u l axis.The u l axis is parallel to the spontaneous polarization P s .(b) Grain l consists of downward (P s ) and upward (−P s ) domain regions.The P s direction in the upward (downward) domains is parallel (anti-parallel) to the u l axis.(c) A scheme of the EKAI model.A polarization nucleus is formed along the u-direction and the domain wall spreads laterally to the u-direction.

Polarization Variation under a Constant Electric Field
In grain l, the polarization direction is tilted from the z-axis by θ l .Nevertheless, it is assumed that, under a constant field, E z , the wall-motion equations are formally the same as Equations ( 1)-( 6) of the KAI model.
Consider a case that a positive constant field E z+ as a step function is applied at t = 0.At t < 0, the volume is occupied fully by the upward regions.Then, the volume fraction of the downward domain (R dn (l)) is varied with time for t ≥ 0 as The volume fraction of the upward domain ( R up (l) is The switching polarization (P z (l)) averaged over the volume of grain l is Consider the opposite case that a negative constant field E z− as a step function is applied at t = 0.At t < 0, the volume is occupied fully by the downward regions.Then, the volume fraction of the upward domain (R up (l)) is varied with time for t ≥ 0 as The switching polarization (P z (l)) averaged over the volume of grain l is

Effect of the Spontaneous Polarization Direction Different from the Z-Axis
The function c(t) in Equations ( 11)-( 15) is 1 − exp −(t/t o ) n , which is the same as Equation (2).However, the u l axis is tilted from the z-axis by θ l .In the EKAI model, the characteristic time for the polarization switching t o is assumed to be i.e., |E z |cos(θ l ) replaces |E z | in Equation (10).

Switching Polarization under Time-Dependent Electric Fields
Let us consider switching polarization evolution in grain l in the case that the electric field varies with time.Using the information of (R dn (l)) now , R up (l) now , and (E z (l)) now , (R dn (l)) next and R up (l) next are derived by the method of this subsection, where t next = t now + ∆t (∆t: a small time-increment), and (• • • ) now or (• • • ) next means the amounts of (• • • ) at t = t now or t = t next , respectively.
The EKAI model assumes that the polarization varies during a short period ∆t as if the polarization varies under a constant field (E z (l)) now .Note, however, that we must consider two cases, i.e., the cases of (E z (l)) now > 0 and (E z (l)) now < 0.
Figure 3 shows an explanation for the case (E z (l)) now > 0. Let us consider a case in grain l that, at t = t now , R dn (l) is (R dn (l)) now under E z (l) = (E z (l)) now .This status is point A in the graph.Draw a curve of the EKAI function (Equation ( 11)) at a constant field, E z (l) = (E z (l)) now , as the blue solid line c(t).Note that c(t) always starts from the origin of the graph (i.e., R dn (l) = 0 at Time = 0.) The line c(t) has (R dn (l)) now at point B. The time at point B is (t Kdn ) now , which is obtained by solving the c(t) (Equation (2)) as We assume that the growth of (R dn (l)) now from t = t now to t = t now + ∆t = t next under E z (l) = (E z (l)) now at point A is the same as the growth of (R dn (l)) now at point B under the same constant field (E z (l)) now .The distance between points A and B is t now − (t Kdn ) now .The R dn (l) growth during ∆t at the point A can thus be calculated by a parallel-shifted function, c(t − t now + (t Kdn ) now ), and thus (R dn (l)) next in the case of (E z (l)) now > 0 is obtained as The volume fraction of the upward domain and the z-axis component of the switching polarization averaged over grain l are and Quite similarly, we derive the following equations in the case of (E z (l)) now < 0. The volume fractions of the upward domains and the downward domains are and (R dn (l) The z-axis component of the switching polarization averaged over grain l is The electric field (E z (l)) next can be obtained using the obtained P z (l) next (Equation (21)  or Equation (25)) and the electrostatic equation, which depends on the device structure considered.See the next section.Once (E z (l)) next is obtained, we know now (R dn (l)) next , R up (l) next , and (E z (l)) next .Then, the quantities of (• • • ) next are regarded as the quantities at t = t now of (• • • ) now ; we can repeat calculation.
When ∆t is infinitesimally small, Equations ( 19) and ( 22) can be expressed as a differential style of R dn (l) and R up (l), that is, for (E z (l)) now > 0, and for (E z (l)) now < 0, The physical meaning of Equations ( 26) and ( 27) is that the differential of c(t − t now + (t Kdn ) now ) at point A in Figure 3 equals the differential of c(t) at point B.

Switching Polarization under Time-Dependent Electric Fields
Let us consider switching polarization evolution in grain  in the case that the electric field varies with time.Using the information of ( ()) ,  () , and ( ()) , ( ()) and  () are derived by the method of this subsection, where  =  + ∆ (∆: a small time-increment), and (⋯ ) or (⋯ ) means the amounts of (⋯ ) at  =  or  =  , respectively.The EKAI model assumes that the polarization varies during a short period ∆ as if the polarization varies under a constant field ( ()) .Note, however, that we must consider two cases, i.e., the cases of ( ()) > 0 and ( ()) < 0. .This status is point A in the graph.Draw a curve of the EKAI function (Equation ( 11)) at a constant field,  () = ( ) , as the blue solid line () ((Equations ( 2) and ( 11)).The point on the curve () at  () =  () is B, and the time at point B is ( ) .We assume that the growth of  () from  =  to  =  + ∆ =  with  () =  () at point A is the same as the growth of  () at point B under the same constant field  () . The  () growth during ∆ at point A can thus be calculated by a parallelshifted function, ( −  + ( ) ), and  () is obtained as Equation (19).At  =  + ∆ the same procedure is repeated under a varied new- (), and we obtain  () at  =  + ∆.The same is performed for  () in the case of  () < 0, and  () is derived as Equation (22).
We assume that the growth of  () from  =  to  =  + ∆ =  under  () =  () at point A is the same as the growth of  () at point B under the same constant field  () . The distance between points A and B is  − ( ) . The  () growth during ∆ at the point A can thus be calculated by a In the case of E z varying with time, let us consider a case in grain l that, at t = t now , R dn (l) is (R dn (l)) now under E z (l) = (E z (l)) now .This status is point A in the graph.Draw a curve of the EKAI function (Equation ( 11)) at a constant field, E z (l) = (E z ) now , as the blue solid line c(t) ((Equations ( 2) and ( 11)).The point on the curve c(t) at R dn (l) = (R dn (l)) now is B, and the time at point B is (t Kdn ) now .We assume that the growth of (R dn (l)) now from t = t now to t = t now + ∆t = t next with E z (l) = (E z (l)) now at point A is the same as the growth of (R dn (l)) now at point B under the same constant field (E z (l)) now .The R dn (l) growth during ∆t at point A can thus be calculated by a parallel-shifted function, c(t − t now + (t Kdn ) now ), and (R dn (l)) next is obtained as Equation ( 19).At t = t next + ∆t the same procedure is repeated under a varied new-E z (l), and we obtain R dn (l) at t = t next + ∆t.The same is performed for R up (l) in the case of E z (l) < 0, and R up (l) next is derived as Equation ( 22).

Total Calculation Scheme for Describing Time-Varying Switching Polarizations of Specific Devices
Figure 4 shows a schematic drawing of (a) an MFM capacitor, (b) an MFIM capacitor, and (c) an MFIS FeFET.As stated before, the ferroelectric film consists of poly-crystal grains indexed by l with l = 1, 2, • • • , l max .We apply a time-dependent voltage, V g , to the top metal electrode of the MFM and MFIM capacitors and to the gate metal electrode of the MFIS FeFET.We provide the ground voltage, 0 V, to the bottom metal electrode of the MFM and MFIM capacitors and to the semiconductor substrate terminal.
Following the EKAI model in the previous section, we derived, at a moment, the z-axis component of the switching polarization P z (l) averaged in the volume of grain l.
We assume that the paraelectric component of the ferroelectric is isotropic.The electric displacement D z (l) averaged in the volume of grain l is represented as where ε o is the permittivity in vacuum and ε f di is the dielectric constant of the paraelectric component of the ferroelectric.
Materials 2024, 17, x FOR PEER REVIEW 10 of 32 Following the EKAI model in the previous section, we derived, at a moment, the zaxis component of the switching polarization  () averaged in the volume of grain .
We assume that the paraelectric component of the ferroelectric is isotropic.The electric displacement  () averaged in the volume of grain  is represented as where  is the permittivity in vacuum and  is the dielectric constant of the paraelectric component of the ferroelectric.and is derived by Equation (38).Then  () at  =  is obtained by Equation (37).For an MFIS FeFET case (c),  at the top surface of the insulator is averaged over all the grains due to a transitional thin layer between the ferroelectric and the insulator (see Figure 5).The surface potential  of the semiconductor is assumed to be uniform because of insufficient acceptor density  .The uniform  and  lead to a grain-independent field  [ () =  for all ]. at  =  is derived by Equation (42).Then,  is obtained from Equation ( 44), and  is obtained from Equation (40).(c) an MFIS FeFET.The EKAI model under a varying field E z with time provides P z (l), R dn (l), and R up (l) at t = t next using R dn (l), R up (l), and E z (l) at t = t now .To repeat the calculation, it is necessary to obtain E z (l) at t = t next .Depending on the structure of the specific devices, the method of deriving E z (l) is different.For an MFM capacitor case (a), E z (l) does not depend on l, and E z at t = t next is obtained by Equation (33).For an MFIM capacitor case (b), the induced charges of the bottom electrode of grain l, −Q(l) is a good approximation.E z (l) depends on l.Q(l) at t = t next and is derived by Equation (38).Then E z (l) at t = t next is obtained by Equation (37).For an MFIS FeFET case (c), D z at the top surface of the insulator is averaged over all the grains due to a transitional thin layer between the ferroelectric and the insulator (see Figure 5).The surface potential ψ s of the semiconductor is assumed to be uniform because of insufficient acceptor density N A .The uniform D z and ψ s lead to a grain-independent field E z [E z (l) = E z for all l].ψ s at t = t next is derived by Equation (42).Then, Q m is obtained from Equation (44), and E z is obtained from Equation (40).
The metal electrodes of all the three devices are assumed to have large free carrier densities.At the metal electrode surfaces facing the ferroelectric or insulator, the electricfield penetrations are negligibly small, and thus no potential variations inside the metal can be assumed.The potential in the metal is uniformly V g .Following Equation ( 28) and Gauss's law, the induced charge per unit area facing the grain l is Because of the averaged value of D z (l), Q(l) is also averaged in the volume of grain l.The total induced charges of the top metal electrode of the potential, V g , is the summation of the induced charge Q(l) over all the grains.The total induced charges are defined as the charge per unit area, Q m , which is where A l is the area that grain l faces the electrode and A tot = ∑ l max l=1 A l .This Q m is the quantities obtained experimentally in Q m − V g measurements, which are explained in (2) in Appendix A. The switching polarization (i.e., not including the paraelectric polarization) averaged over the entire ferroelectric, P z , is as Using the EKAI and electrostatic equations across MFM, MFIM, and MFIS stacks, calculation steps are expressed as follows.The metal electrodes of all the three devices are assumed to have large free carrier densities.At the metal electrode surfaces facing the ferroelectric or insulator, the electricfield penetrations are negligibly small, and thus no potential variations inside the metal can be assumed.The potential in the metal is uniformly  .Following Equation ( 28) and Gauss's law, the induced charge per unit area facing the grain  is Because of the averaged value of  (),  () is also averaged in the volume of grain  .The total induced charges of the top metal electrode of the potential,  , is the summation of the induced charge  () over all the grains.The total induced charges are defined as the charge per unit area,  , which is where  is the area that grain  faces the electrode and  = ∑  .This  is the quantities obtained experimentally in  −  measurements, which are explained in (2) in Appendix A. The switching polarization (i.e., not including the paraelectric polarization) averaged over the entire ferroelectric,  , is as Using the EKAI and electrostatic equations across MFM, MFIM, and MFIS stacks, calculation steps are expressed as follows.

MFM Capacitor
The electrostatic equation of the MFM capacitor is The approximate structure assumed in the EKAI.In the experiment, the insulator is the bilayer of a 4 nm thick HfO 2 layer and a 2.6 nm thick IL (SiO 2 ) layer.The ferroelectric CSBT is 135 nm thick.A transitional layer exists at the interface between the CSBT and HfO 2 layers.A photo of a cross-section TEM confirmed this layer presence [58].The transitional layer (≈5 nm thick) is constituted by fine grains.The TEM photo contrast indicates that the main atom elements are of the CSBT.Since a CSBT-originated high permittivity can be assumed, D z (l) is expected to be averaged of this transitional layer.Therefore, we can approximate the structure as (b) in this figure.D z is averaged over all the grains at the top surface of the insulator.

MFM Capacitor
The electrostatic equation of the MFM capacitor is where V f is the voltage across the ferroelectric layer.The flat band voltage, V f b , is the work function difference between the metal top electrode and the metal bottom electrode.Since V f does not depend on the grains, the electric field is also independent of l as We have Equation (30) as The second term of Equation ( 34) is the switching polarization averaged over the entire ferroelectric, P z (Equation ( 31)).
Equations ( 21) and ( 25) of the EKAI model in the proceeding section can give us the P z (l) (for l = 1 • • • l max ) at t = t next .We also know V g at t = t next , because V g is externally applied.Then, E z and Q m at t = t next are obtained by Equation (33) and Equation (34), respectively.We know the behavior of the MFM capacitor at t = t next .
R dn (l), R up (l), and E z (l) at t = t next are regarded as those at t = t now .The EKAI model provides R dn (l) and R dn (l), and P z (l) at t = t next .Equation (33) gives us E z (l) at t = t next .We repeat these procedures and obtain fully the time-dependent solution of the MFM capacitor.

MFIM Capacitor
Charges appear at the top surface of the bottom electrode.We shall extend the z-axisorientated boundaries between the adjacent grains into the insulator.We have the potential and E z relationships for each grain; i.e., we have the following for l = 1, 2, • • • , l max : and Here, V f (l) is the voltage across the grain l.V i (l) is the voltage across the region in the insulator belonging to grain l.At the top surface of the bottom electrode corresponding to grain l, the induced charge is −Q(l), and we have for l = 1, 2, • • • , l max : with the insulator capacitance, C i : and We eliminate E z (l) from Equations ( 35)- (37).Then, we derive with C f di is the paraelectric component capacitance of the ferroelectric as The EKAI model provides R dn (l) and R dn (l), and P z (l) at t = t next .We also know V g at t = t next , because V g is externally applied.Equation (38) leads to Q(l) at t = t next for all l, and Equation (37) gives us E z (l) at t = t next , using the obtained Q(l) for all l.All physical quantities are derived at t = t next .Then, similarly to the case of the MFM capacitor, we can return the EKAI model calculation, and consequently we derive the time-dependent numerical solution for the MFIM capacitor.

MFIS FET
Figure 6 summarizes the total calculation scheme for MFIS FeFETs.Regarding the semiconductor, in the MFIS stack, the acceptor concentration (N A ) and carrier density in the semiconductor are much smaller than the carrier density in the metal layer.Thus, the electric field penetration into the semiconductor inevitably occurs so that the semiconductor surface potential ψ s is different from the substrate potential (V sub = 0).Imagine further a surface case where a region of an inversion state adjoins a region of an accumulation state.The surface potential cannot change abruptly on the surface.A good measure of the potential variation on the surface is the maximum depletion width, where ε s is the relative permittivity of the semiconductor, n i the intrinsic carrier concen-tration, ζ = e o /(k B T), and e o the elementary charge [56].For example, W m ∼ = 300 nm when N A = 10 16 cm −3 T = 300 K using ε s = 11.9 and n i = 1.45 × 10 10 cm −3 for Si.If the size of ferroelectric grains is smaller than W m , a uniform surface potential ψ s over all the grains is a good approximation (Figure 4c). Figure 5 shows a schematic drawing of a ferroelectric-insulator-semiconductor (FIS) as part of the experimental MFIS FeFETs.As shown in Appendix A, the gate stack of the experimental FeFETs is Ir/CSBT/HfO 2 /Si.Via the crystallization annealing of the CSBT, a 2.6 nm thick SiO 2 interfacial layer (IL) was formed between HfO 2 and Si.The HfO 2 and CSBT layers are 4 nm thick and 135 nm thick, respectively.The bilayer of the HfO 2 and IL forms the insulator in the MFIS.Transmission electron microscope photos confirmed a thin (about 5 nm thick) transitional layer as shown in Figure 5a [58].The transitional layer is constituted by fine grains (≈5 nm) whose main elements are originated from the ferroelectric CSBT.Due to the fine sizes, these grains may be non-ferroelectric but work as a high permittivity material.The dielectric constants of this transitional layer, the HfO 2 layer, and the SiO 2 IL are typically 180, 25, and 3.9, respectively.The different values of P z (l) (Equation ( 28)) at the bottom of the ferroelectric among the grains are averaged via the in-plane pass in the transitional layer (Figure 5b).Therefore, the switching polarization at the interface between the CSBT and the insulator consisting of the bilayer of HfO 2 and IL can be reasonably assumed to equal the z-axis component of the polarization P z (l) averaged over all the grains, P z , that is already defined in Equation (31).Since P z at the top surface of the insulator and ψ s at the bottom of it are uniform laterally, the electric field E i along the z-axis in the insulator and the potential V i across the insulator are also laterally uniform, which leads to a grain-independent field E z [E z (l) = E z for all l] in the ferroelectric and a grain-independent potential V f = E z d f across the ferroelectric.The z-axis component of the electric displacement D z averaged over all the grains as well as Q m from Equations ( 28)-( 30) is FeFETs have the gate, source, drain, and substrate terminals.Let us consider the case where the gate voltage V g is generally varied with time and other terminals are grounded (V s = V d = V sub = 0).Here V s , V d , and V sub are the voltages applied on the source, drain, and substrate, respectively.The electrostatic equation across the MFIS stack is where V f b is the flat-band voltage, i.e., the difference between the metal work function and semiconductor fermi level.By eliminating E z from Equations ( 40) and ( 41) and using Gauss's law while noticing the capacitance definition of Equations (36b) and (39), we have If there is no interface state density between the insulator and semiconductor, the induced charge density, Q s , in the semiconductor surface region has the same magnitude as Q m with the opposite polarity (Q m = −Q s ).Q s is a function of semiconductor surface potential ψ s as follows [59]: with a negative sign for ψ s > 0. This equation is for n-channel FeFETs formed in p-type substrates.The equation for p-channel FeFETs can be rewritten appropriately.The Debye /2 and The donor, acceptor, and intrinsic carrier densities in the semiconductor are N D , N A , and n i , respectively, and p po and n po are the equilibrium densities of holes and electrons, respectively.
where  is the flat-band voltage, i.e., the difference between the metal work function and semiconductor fermi level.By eliminating  from Equations ( 40) and ( 41) and using Gauss's law while noticing the capacitance definition of Equations (36b) and ( 39), we have If there is no interface state density between the insulator and semiconductor, the induced charge density,  , in the semiconductor surface region has the same magnitude as  with the opposite polarity ( = − ). is a function of semiconductor surface potential  as follows [59]: with a negative sign for  > 0. This equation is for n-channel FeFETs formed in p-type substrates.The equation for p-channel FeFETs can be rewritten appropriately.The Debye length is  =      ⁄ with  =  −  + (( −  ) + 4 ) ⁄ /2 and  =   ⁄ .The donor, acceptor, and intrinsic carrier densities in the semiconductor are  ,  , and  , respectively, and  and  are the equilibrium densities of holes and electrons, respectively.
In the case that the interface states between the semiconductor and insulator are considered, the equation  = − is modified as where  is the trapped charge at the interface per area that is expressed as Here, ℰ is an electron-energy variable,  (ℰ) is the area density of interface-states per electron energy, and  is the Fermi-Dirac distribution function for acceptors [59].In the case that the interface states between the semiconductor and insulator are considered, the equation where Q it is the trapped charge at the interface per area that is expressed as Here, E is an electron-energy variable, N it (E ) is the area density of interface-states per electron energy, and F SA is the Fermi-Dirac distribution function for acceptors [59].
In the case that N it (E ) is approximated as a constant D it [V −1 cm −2 ] with respect to the energy, Q it can be approximated as For convenience of numerical root-finding calculations, we emphasize that Q s (Equation ( 43)) monotonically decreases with the increase of ψ s , and Q it (Equation (45) or Equation ( 46)) also has the same monotonical property.Equation (44) shows that Q m is a monotonically increasing function of ψ s .
See Equation (42).The right-hand side quantity is a monotonically increasing function of ψ s .At t = t next in the previous section, we derived P z (l) for all l by either Equation (21)  or Equation (25).V g at t = t next is given by an external voltage source.Thus, the left-hand side of Equation ( 42) is a known constant at t = t next .The right-hand side of Equation ( 42) consists of only one variable ψ s and is a monotonically increasing function of ψ s .Hence, ψ s can be uniquely determined at t = t next .Once ψ s is determined, we derive Q m from Equations ( 43)- (46), and E z from Equation (40) at t = t next .
We now have E z at t = t next , and, in the previous section, we have R dn (l) and R up (l) at t = t next for all l.The set of E z , R dn (l), and R up (l) at t = t next replaces a set of E z , R dn (l), and R up (l) at t = t now .If E z > 0, R dn (l) and R up (l) at t = t next are obtained by Equations ( 19) and (20).If E z < 0, R up (l) and R dn (l) at t = t next are obtained by Equations ( 22) and (23).The procedure of this section gives us E z at t = t next .Repeating these procedures provides the dynamics of the MFIS FeFETs.
The status of the polarization in ferroelectric or the electric displacement can be monitored by the drain current, I d , of the sub-threshold region, which is represented as functions of the drain voltage V d and the surface potential ψ s by [59]: In this paper, all numerical results of FeFETs are obtained on the premise of V s = V d = V sub = 0, whereas the practical I d measurements of n-channel FeFETs are usually on the condition of V d = 0.1 V and V s = V sub = 0.Such a small difference in V d condition does not affect the validity of comparing results from the calculations and the measurements.
The EKAI, the formulae from Equations ( 11)-( 27), and the total calculation scheme for MFIS FeFETs, those from Equations ( 28)-( 31) and from Equations ( 40)-( 47), describe the general transient response of MFIS FeFETs under time-dependent V g conditions.The formulae also cover slowly changing phenomena.There is a so-called data retention mode, in which all the quantities such as Q m , E z , R dn (l), and R up (l) vary very slowly, and the ∆t values in Equations ( 19) and ( 22) can be chosen flexibly with minimal changes of those quantities; thus, retention results of the period days and years are obtained in a practical computation time.

Parameters Being Able to Be Assigned
The EKAI model and the total calculation scheme of this paper are verified using experimental data of FeFETs consisting of MFIS gate stacks of Ir/CSBT/HfO 2 /IL/Si.The experimental details are reviewed in Appendix A. An insulator in the modeled MFIS FeFET corresponds to the bilayer of IL and HfO 2 in a real CSBT FeFET.Capacitance of the insulator (C i ) was evaluated as C i = 0.99 µF/cm 2 .Regarding the ferroelectric, the dielectric constant ε f di of the paraelectric component is determined by the Q m -V g curves at V g < 0 or in the third quadrant of Figure A2 because the semiconductor depletion layer does not affect at V g < 0. In the third quadrant, by taking the gradient of the curve of a V g sweep amplitude, the combined capacitance of C i and C f di (Equations (36b) and ( 39)) was 0.54 µF/cm 2 .Since C i = 0.99 µF/cm 2 and d f = 135 nm, ε f di = 180 at room temperature was obtained and used for calculation in the EKAI.The P s direction of SBT and CSBT ferroelectrics is the a-axis direction of the crystal unit cell of SBT and CSBT [60,61].
With regard to the angle θ l and the area A l , the technique of EBSD patterns could characterize the crystal orientation and grain size.The bar graph plots in Figure 7 show the distribution function of grains with the orientation angle θ.The quantity of the vertical axis is the area of grains whose θ is in the range from θ to θ + ∆θ (∆θ: the bar width).The scanning area of EBSD, ≈ 56 µm 2 , is not enough for statistical treatment, meaning that ∆V th versus log(t w ) curves simulated by the EKAI model and the calculation scheme using the bare bar graph plots were not smooth.As shown in Figure A4, the experimental ∆V th − log(t w ) curves are very smooth.Therefore, we used a fitted smooth curve in Figure 7 instead of the bare data.Although the smoothed curve is used in the actual calculation, the curve is digitized at every 3 • to save calculation time.
When we make quasi-static calculation of I d -V g and Q m -V g , a sinusoidal V g function is provided as V g = V f b + V amp sin(2π f s t).The frequency f s and the calculation time length are typically set as 10 Hz and 0.2 s, respectively.The time step ∆t = 1 × 10 −9 s is commonly used in the calculation of I d -V g , Q m -V g , and PWVR.Validity of the value ∆t and accuracy of the calculated results were verified by confirming that calculations with ∆t = 0.5 × 10 −9 s gave the same results as those with ∆t = 1 × 10 −9 s.
Experimental I d -V g curves do not suggest any electron-energy dependence of N it (E ) in Equation (45), and thus we use D it and Equation (46) regarding the interface states between the insulator and semiconductor.Several series of I d -V g simulations with various V f b and D it values were examined and compared to the experimentally obtained I d -V g .We found that D it = 4 × 10 12 V −1 cm −2 is a good value to fit the experiment.Regarding V f b , the suitable range was between −0.8 V and −1.0 V, and thus we used V f b = −0.8V in this work.
A4, the experimental ∆ − log( ) curves are very smooth.Therefore, we used a fitted smooth curve in Figure 7 instead of the bare data.Although the smoothed curve is used in the actual calculation, the curve is digitized at every 3° to save calculation time.
When we make quasi-static calculation of  - and  - , a sinusoidal  function is provided as  =  +  sin(2 ).The frequency  and the calculation time length are typically set as 10 Hz and 0.2 s, respectively.The time step ∆ = 1 × 10 −9 s is commonly used in the calculation of  - ,  - , and PWVR.Validity of the value ∆ and accuracy of the calculated results were verified by confirming that calculations with ∆ = 0.5 × 10 −9 s gave the same results as those with ∆ = 1 × 10 −9 s.Experimental  - curves do not suggest any electron-energy dependence of  (ℰ) in Equation ( 45), and thus we use  and Equation (46) regarding the interface states between the insulator and semiconductor.Several series of  - simulations with various  and  values were examined and compared to the experimentally obtained  - .We found that  = 4 × 10 12 V −1 cm −2 is a good value to fit the experiment.Regarding  , the suitable range was between −0.8 V and −1.0 V, and thus we used  = −0.8V in this work.
According to the PFM experiments and the analyses by phenomenological theories including random disorder potentials originated from imperfections in epitaxial films [53], the exponent σ of  in Equation (8b) and that of  in Equations (9b) and ( 10) are generally σ ≠ 1.The magnitude of σ depends on materials and material preparation methods [62].Some experiments indicated σ ≈ 0.5 of BaTiO3 [63], 0.5-0.6 in PbZrxTi1−xO3 (PZT) of x = 0.2 [64], and 0.20-0.28 of ferroelectric organic polymer [65].In the case that the defect densities were intentionally increased, σ was decreased more [53].However, other experiments showed σ ≈ 0.9 and ≈1.0 for epitaxial PZT films [21,55].Furthermore, the experiments switching current to single crystals of BaTiO3 and triglycine sulfate According to the PFM experiments and the analyses by phenomenological theories including random disorder potentials originated from imperfections in epitaxial films [53], the exponent σ of v wall in Equation (8b) and that of t o in Equations (9b) and ( 10) are generally σ ̸ = 1.The magnitude of σ depends on materials and material preparation methods [62].Some experiments indicated σ ≈ 0.5 of BaTiO 3 [63], 0.5-0.6 in PbZr x Ti 1−x O 3 (PZT) of x = 0.2 [64], and 0.20-0.28 of ferroelectric organic polymer [65].In the case that the defect densities were intentionally increased, σ was decreased more [53].However, other experiments showed σ ≈ 0.9 and ≈1.0 for epitaxial PZT films [21,55].Furthermore, the experiments switching current to single crystals of BaTiO 3 and triglycine sulfate showed an exponential relationship without the exponent σ [51,52].The theory of wall propagation in a defect-free crystal by Miller and Weinreich [20] led to an exponential form like Equation (8b) but did not include σ.Molecular dynamics simulation for a defect-free domain interface also demonstrated v wall as a form of Equation (8b) but without σ [22].Zhao et al. [66] showed that polarization switching times of ferroelectric organic polymer films obeyed simply the Merz exponential law [51].We did not use a method for intentionally introducing defects during annealing processes for ferroelectric layer crystallization of MFIS FeFETs [3][4][5].Like these, there is no reason to choose σ ̸ = 1 in the following calculations for the SBT-based FeFETs.σ = 1 is assumed in simulations shown later in this paper, but σ as the mathematical expression is maintained.The σ appears in Section 5.2.later.We adopted n = 1.3, which is a value derived experimentally for less than 200 kV/cm for an epitaxial PZT film [37] and which is also close to the average value, 1.25, of the n range, 1.0-1.5, reported for polymer ferroelectric films [66].Note also that the calculated results were insensitive to the n variation from 1 to 3. The acceptor density was also not sensitive to the results, and N a = 1 × 10 16 cm −2 is used for calculation.

Method for Determining Significant Three Parameters
The remaining parameters, E act , t in f , and P s can be determined by curve fitting of numerical results to the experimental about the ∆V th vs. t w in PWVR.The experimental results of ∆V th vs. t w are found in Figure A4 in Appendix A. Regarding the numerical results, PWVR simulations with varying E act , t in f , and P s are introduced in Figure 8a, Figure 8b and Figure 8c, respectively.Every marker corresponds to a calculated point of ∆V th vs. t w .Using the cases of Figure 8, we show how uniquely three parameters, E act , t in f , and P s , are determined.We present a reference curve (the blue solid line with filled square markers in Figure 8a, Figure 8b and Figure 8c, respectively), which was a numerical solution for well simulating an experimentally obtained ∆V th vs. t w .The reference curve was drawn using a set of parameters E act = E o , t in f = t 1 , and P s = P o with V h = −V l = 4 V.The E o , t 1 , and P o are constants for the sake of explanation.The E act and t in f are deeply involved in ferroelectric polarization dynamics.The P s is an inherent parameter of the ferroelectric.As shown in Figure 8a, if E act is as large as 1.7 E o , the increase of ∆V th with t w is very slow, and the curve is far from the reference curve in a realistic t w range of the experimental.If E act is as small as 0.43 E o , rapidly increasing ∆V th approaches a saturated value.This is far from the log-linear styles.An appropriate ∆V th can draw a log-linear curve like the reference curve at E act = E o .As shown in Figure 8b, t in f determines a quantity of the ∆V th vs. log(t w ) curve shift in parallel along the log(t w ) axis.The polarization growth in Equation (21) or Equation ( 25) is proportional to P s , meaning that if P s is larger, the separation between a ∆V th vs. log(t w ) curve written by pulses of V h (= −V l ) and the neighboring curve written by V h ± 1V(= −V l ∓ 1V) is wider.See Figure 8c, where the curves are drawn in the case of P s = P o and 1.7P o , respectively.The role of P s is to adjust this separation distance to fit the experimental results.By experiencing these processes, E act , t in f , and P s are uniquely determined under a fixed σ.Here, "uniquely determined" means that no solutions having a new parameter set exist far from the derived parameter set.Remember that  has a meaning of an activation-or a threshold-field for domain wall motions as Equations (9b) and ( 10) indicate.The domain wall energy is affected by elastic and electric-dipole contributions in atomic scales.Theoretical works [20,24,25,53,67] indicated that the domain wall energy included a power exponent of  , indicating that the  may also be a function of  .Since the treated temperature is only room temperature, three parameters,  ,  , and  , can be searched independently to fit the experimental data.In the case that temperature is varied,  is also changed with temperature.We must consider that  is a function of  via the domain wall energy in addition to the thermally activation term ( ) (Equation (9b)).Remember that E act has a meaning of an activation-or a threshold-field for domain wall motions as Equations (9b) and ( 10) indicate.The domain wall energy is affected by elastic and electric-dipole contributions in atomic scales.Theoretical works [20,24,25,53,67] indicated that the domain wall energy included a power exponent of P s , indicating that the E act may also be a function of P s .Since the treated temperature is only room temperature, three parameters, P s , E act , and t in f , can be searched independently to fit the experimental data.In the case that temperature is varied, P s is also changed with temperature.We must consider that E act is a function of P s via the domain wall energy in addition to the thermally activation term (k B T) −1 (Equation (9b)).

Calculation Using Optimized Parameters and Comparison with the Experimental Data
Figure 9 shows a fitting result of PWVR where the calculated curves compared to the experimental data at V h = 3 V, 4 V, 5 V, and 6 V. Markers represent the calculated points of ∆V th vs. t w using an optimum parameter set of E act = 828 kV/cm, t in f = 8.30 × 10 −12 s, and P s = 3.0 µC/cm 2 .Solid lines mean the experimental results of Figure A4 introduced in Appendix A. The calculation well simulated the experimental ∆V th vs. log(t w ) throughout the wide t w range from 50 ns to 0.5 ms.Other parameters for the calculation are summarized as follows: σ = 1, n = 1.3, V f b = −0.8V, D it = 4 × 10 12 V −1 cm −2 , N A = 1 × 10 16 /cm 3 , d f = 135 nm, ε f di = 180.d i = 3.5 nm, and ε i = 3.9.
Using the same values of parameters as those solved for fitting PWVR (Figure 9), other correlations were simulated that were quasi-static characteristics of I d -V g (Figure 10) and Q m -V g (Figure 11) with various V g sweep amplitude.The calculated results are drawn with thick and red-colored lines.The experimental results are expressed by thin lines colored in black.In the calculations, V th was defined as V g at which the semiconductor surface potential ψ s was equal to 85% of 2ψ B , the surface strong inversion condition (i.e., ψ s = 2ψ B × 0.85 where ψ B = (1/ζ)ln(N A /n i ) [59]).Figures 10 and 11 show moderate agreements of the calculated results with the experimental results.However, high curvature of Q m -V g around V g = 0 V seems more emphasized in the calculated than in the experimental results as shown in Figure 11.As an effort of the Q m -V g matching, for example, D it may be raised for decreasing the nonlinearity of Q m -V g .But the attempt enhances another mismatch in I d -V g as shown in Figure 10.The reason for the inconsistency is not clear now.Despite having some numerical mismatch remaining in the curve fitting, the EKAI model and the calculation scheme qualitatively and comprehensively well simulate FeFET characteristics such as dynamic PWVR and quasi-static I d -V g and Q m -V g .
dependences of Δ vs.  curves, respectively.As shown in (a), if  is not optimized, the loglinear relation of Δ vs.  cannot be derived across the  wide range from 50 ns to 1 ms.(b) indicates that  variation brings a parallel shift of Δ vs. log( ) curves.(c) indicates that the separation of Δ vs.  lines at  = − = 3 V, 4 V, and 5 V is very sensitive to  .The line separation by the dashed lines is larger than the solid lines, where the former  is 1.7 times larger than the latter. in (a),  in (b), and  in (c) are example constants for explanation convenience.
Remember that  has a meaning of an activation-or a threshold-field for domain wall motions as Equations (9b) and (10) indicate.The domain wall energy is affected by elastic and electric-dipole contributions in atomic scales.Theoretical works [20,24,25,53,67] indicated that the domain wall energy included a power exponent of  , indicating that the  may also be a function of  .Since the treated temperature is only room temperature, three parameters,  ,  , and  , can be searched independently to fit the experimental data.In the case that temperature is varied,  is also changed with temperature.We must consider that  is a function of  via the domain wall energy in addition to the thermally activation term ( ) (Equation (9b)).Using the same values of parameters as those solved for fitting PWVR (Figure 9), other correlations were simulated that were quasi-static characteristics of  - (Figure 10) and  - (Figure 11) with various  sweep amplitude.The calculated results are drawn with thick and red-colored lines.The experimental results are expressed by thin lines colored in black.In the calculations,  was defined as  at which the semiconductor surface potential  was equal to 85% of 2 , the surface strong inversion condition (i.e.,  = 2 × 0.85 where  = (1  ⁄ )(  ⁄ ) [59]).Figures 10   and 11 show moderate agreements of the calculated results with the experimental results.However, high curvature of  - around  = 0 V seems more emphasized in the calculated than in the experimental results as shown in Figure 11.As an effort of the  - matching, for example,  may be raised for decreasing the nonlinearity of  - .But the attempt enhances another mismatch in  - as shown in Figure 10.The reason for the inconsistency is not clear now.Despite having some numerical mismatch remaining in the curve fitting, the EKAI model and the calculation scheme qualitatively and comprehensively well simulate FeFET characteristics such as dynamic PWVR and quasi-static  - and  - . measurements for various  sweep amplitude.Red lines and dark grey lines are the calculation and experiment results, respectively.We used the same parameters as those for PWVR calculation in Figure 9.

Calculation Using Optimized Parameters and Comparison with the Experimental Data
Figure 10.Model calculation with the experiment (Figure A1) of quasi-static I d vs. V g measurements for various V g sweep amplitude.Red lines and dark grey lines are the calculation and experiment results, respectively.We used the same parameters as those for PWVR calculation in Figure 9.  measurements for various  sweep amplitude.Red lines and dark grey lines are the calculation and experiment results, respectively.We used the same parameters as those for PWVR calculation in Figure 9.In the EKAI model, a Q m vs. E z correlation is calculated along a hysteresis loop as shown in Figure 12a.The drawing sequence is explained by corresponding V g variations with checkpoints as shown in Figure 12b, which are, a', b, c, d, g, g' d', e, f, a, h, h' and back to a', repeated cyclically in this order.A positive pulse writing (PPW) draws a trajectory connecting a', b, c, and d.A negative pulse writing (NPW) draws a trajectory connecting d', e, f, and a.As stated in the introduction, we assume in this paper that the parameters except for the polarization switching varies instantaneously.On the Q m -E z curve, points instantaneously move from a' to b, from c to d, from d' to e, and from f to a.These curves can be written as where P zconst is a constant that each straight line has and P zconst equals P zb , P zc , P ze , and P z f for the line a'-b, c-d, d'-e, f-a, respectively.

Details of PWVR Operations Using 𝑄 vs. 𝐸 Domains
In the EKAI model, a  vs.  correlation is calculated along a hysteresis loop as shown in Figure 12a.The drawing sequence is explained by corresponding  variations with checkpoints as shown in Figure 12b, which are, a', b, c, d, g, g' d', e, f, a, h, h' and back to a', repeated cyclically in this order.A positive pulse writing (PPW) draws a trajectory connecting a', b, c, and d.A negative pulse writing (NPW) draws a trajectory connecting d', e, f, and a.As stated in the introduction, we assume in this paper that the parameters except for the polarization switching varies instantaneously.On the  - curve, points instantaneously move from a' to b, from c to d, from d' to e, and from f to a.These curves can be written as where  is a constant that each straight line has and  equals  ,  ,  , and  for the line a' -b, c-d, d'-e, f-a, respectively.
As discussed in Section 3.3, the electrostatic potential equation across the MFIS stack (Equation ( 41)) is valid at any moment where  =   and  =   ⁄ .The  is a function of  by a discussion at Equations ( 43)-( 46).Consequently, at any time,  - satisfies the following equation: Curves described by Equation ( 49) are called load lines in this paper.Four load lines I, II, III, and IV appear in Figure 12a.The load lines I, II, III, and IV are the lines when  in Equation ( 49) equals  , 0,  , and  _ , respectively.Equation (49) indicates that a solution point ( ,  ) locates at any moment ( ) on a load line having the  value at  =  .The EKAI model decides which point on the load line is really the solution point.
The point b, d, e, and a are decided by that Equations ( 48) and ( 49) are simultaneously satisfied.From b to c in PPW, the EKAI model describes a  increase from  to  during the period  of  application (Equation ( 19)).Similarly, from e to f in NPW, the model describes a  decrease from  to  during the period  of  application (Equation ( 22)).
grows on b-c according to Equation (19), and R up (l) grows on e-f according to Equation ( 22) by the EKAI.The line d-g-g'-d' or a-h-h'-a' corresponds to data retention (or holding) and V th read stage.At the read end, Q m rises to g' or h'.At the retention stage the trajectory also exists on the load line Due to the presence of an electric field, the polarization may be decreased, i.e., d (or a) may move to an inside point d' (or a') on the load line at V g = 0.The decreasing extent depends on E act , E z , cos(θ l ), and the time length between d and d' (or a and a').Q m rising to h' for the read operation may also decrease the polarization.
As discussed in Section 3.3, the electrostatic potential equation across the MFIS stack (Equation ( 41)) is valid at any moment where 43)- (46).Consequently, at any time, Q m -E z satisfies the following equation: Curves described by Equation ( 49) are called load lines in this paper.Four load lines I, II, III, and IV appear in Figure 12a.The load lines I, II, III, and IV are the lines when V g in Equation ( 49) equals V h , 0, V l , and V swp_end , respectively.Equation (49) indicates that a solution point (E z , Q m ) locates at any moment (t a ) on a load line having the V g value at t = t a .The EKAI model decides which point on the load line is really the solution point.
The point b, d, e, and a are decided by that Equations ( 48) and ( 49) are simultaneously satisfied.From b to c in PPW, the EKAI model describes a P z increase from P zb to P zc during the period t w1 of V h application (Equation ( 19)).Similarly, from e to f in NPW, the model describes a P z decrease from P ze to P z f during the period t w2 of V l application (Equation ( 22)).
Initially, idling write cycles are executed that consist of PPW and NPW without a V th reading (VR).The idling cycles have the role of making the Q m -E z trajectory converge into the steady loop shown in Figure 12a.The simulation can be started from the coordinate origin at t = 0, where E z = 0, Q m = 0, and R dn (l) = R up (l) = 1/2 for all l.The Q m -E z trajectory changes during the idling write cycles and is converged into a steady state loop after experiencing plural cycles.Then the PWVR operation starts.After one PPW is executed, a VR draws a trajectory connecting d, g, g', and d'.After one NPW is executed, another VR draws a trajectory connecting a, h, h', and a'.At read, V g is swept to V swp_end (i.e., V g at point g' and h').When V g = V swp_end , g' and h' are on the load line IV.V th values are decided at a reference level of Q m , i.e., at points g and h.Q m and ψ s have a single-valued function relationship with each other via Equations ( 43)- (45).The reference level of ψ s (ψ th ) can replace the reference of Q m .The ψ th is chosen in a sub-threshold region of the FeFETs so that ψ B ≤ ψ th ≤ 2ψ B .The EKAI model simulation indicates that, during a V g sweeping after NPW, there is a tendency that P z increases (i.e., the stored negative P z f decreases).If this is a visible case, point a' separates from a and shifts to the E z smaller direction on load line II.The EKAI model works during all the period d-g-g'-d' and a-h-h'-a', which means that this period is also in a data-retention stage.If the depolarization electric field is not small, Q m shifts to the |Q m | decreasing direction on load line II.In this sense, d' and a' differ from d and a, respectively.The decrease amounts depend on the relationship among E act , E z , and cos(θ l ), as shown in Equation (16).

Strategy to Reduce Charge Injections into FeFET Gate Stacks
Note that, in drawing the loop in Figure 12a, it passes through two points: the maximum induced charge (Q mmax ) and the minimum one (Q mmin ) (i.e., the negative maximum one).The Q mmax and −Q mmin decide the amount of an undesirable current of the direct tunneling type or the electric-field-assisted (Fowler-Northeim, FN) tunneling type.At is defined where the V imax is the maximum voltage drop across the insulator.If Q mmax = 2.0 µC/cm 2 and the insulator is 1.6-nm-thick SiO 2 , then V imax = 0.93 V and the corresponding field maximum of the insulator, E imax = 5.8 MV/cm.Let us estimate gate leakage currents thorough the insulator.Using a high permittivity insulator or a combination of such insulators weakens the field a little, but the field is still high enough to induce the leakage currents.Investigations of the tunneling current of polysilicon/SiO 2 /Si [68,69] are good references to know the impact of the charge injection for the MFIS stacks where the semiconductor is Si.The charge injection is mainly caused by the tunneling current through the IL (i.e., SiO 2 ).We propose a significant guideline in investigating ferroelectric FETs.Imagine what happens at Q m = Q mmax = 2.0 µC/cm 2 .At this moment, the charge injection is the highest.According to Ref. [68], in the case that SiO 2 was about 1.6 nm thick, the tunneling current was ≈ 10 A/cm 2 at the 5.8 MV/cm whereas it was ≈10 −5 A/cm 2 in the case of about 3.2 nm thick SiO 2 , at the same field.The voltage drops across the SiO 2 layer, Q m /C SiO2 , [C SiO2 : the SiO 2 layer capacitance] at Q m = 2.0 µC/cm 2 are 0.93 V and 1.85 V for the 1.6 nm thick and 3.2 nm thick SiO 2 layers, respectively.By choosing this twice thick insulator SiO 2 , the tunneling current decreases ≈10 −6 times, and the V g for writing increases only about 0.9 V.As this quantitative consideration indicates, there seems no way to avoid the tunneling current except for increasing the thickness for SiO 2 .To preserve the nonvolatile device reliabilities, the IL SiO 2 should be moderately thick enough to avoid charge injection caused by the tunneling current.A strategy of thinning SiO 2 is logically failed.
If charge injections are not negligible, a scenario is as follows: In PPW, electrons are injected from the silicon.The electrons may mostly arrive at the metal electrode and be absorbed.However, some of them are trapped in the ferroelectric layer, the insulator, and the interface between the ferroelectric layer and insulator.The trapped ones near the silicon side may return to the semiconductor by tunneling back after PPW [45], but other trapped electrons remain trapped, leading to the increase of V th of n-channel FeFETs, while NPW holes are injected from the silicon.Similarly, some holes are stably trapped, leading to the decrease of V th of n-channel FeFETs.Since the number of trapped electrons after PPW and that of trapped holes after NPW are not the same, unintended V th shifts appear with increasing the cycle of PPW and NPW in endurance tests.

Short Consideration in Negative-Capacitance-Transistors
On load lines I, II, III, and IV mentioned above, P z -and Q m -increase accompanies E z decrease, meaning that a capacitance dQ m /dE z is negative.Note that, as described in the EKAI model and the total calculation scheme, the time response of the linear dielectric part is instantaneous, but P z variation needs time.Thus, while V g is swept back and forth between negative and positive voltages, the derived Q m -V g and I d -V g curves inevitably draw hysteresis loops.The EKAI model does not realize the idea discussed in Ref. [70] regarding steep-slope non-hysteresis transistors during on-off operation for back-and-forth voltage sweeping.

Coercive Field in the EKAI Model
Although the coercive fields can be read in the Q m -V g hysteresis curves derived by the calculation of this model, the EKAI model of this paper does not contain the coercive field (E c * ) as an explicit parameter.Let us find the relationship between the model and E c * for a case of a single grain ferroelectric indexed by θ.To take Q m vs. E z curves, a triangular shape E z as a function of time, like Figure 13b, whose linear slope is K, is supplied to an MFM capacitor.Figure 13a shows a P z vs. E z curve.During E z increasing with time, R dn varies from 0 to 1, like Figure 13c.Since E z = K time + constant, R dn vs. time is converted to R dn vs. E z .At a narrow region of E z , R dn increases rapidly, as shown in Figure 13c.The R dn rapid increase corresponds to P z increase of the right-side branch of the hysteretic loop in Figure 13a, indicating that an E z at which the R dn rapid increase occurs is regarded as a coercive field.Since the R dn increase is rapid within a narrow range of E z , a defined coercive field well approximates E c * which is defined as the field at Q m = 0 (Figure 13a).The E z representing this narrow range field can be decided by the E z (E cm ) at which dP dn /dE z takes a maximum (Figure 13d).E cm is the coercive field derived analytically in the EKAI model.Exactly speaking, E cm approximates E c * , but E cm is not equal to E c * .
By starting from Equation ( 26) with the use of dE z = Kdt and by calculating the second derivative, d 2 R dn /dE z 2 = 0, we have, without any approximation, In Equation (50), K is the sweep slope, and E act , t in f and σ are the model simulation parameters.In g(n), n is the model parameters.Although g(n) contains R dn , g(n) weakly depends on R dn variation.This means g(n) can be regarded as a constant.(In the case of n = 1, g(n) equals 1, and Equation (50) does not contain the dimension parameter, n.) Hence, Equation (50) can be solved for E z , and the root of it is E cm .
Equation ( 50) is a finding-root problem of x e = E act /E z cosθ on the condition of g(n) = 1.Since E act is included in x e , E cm is changed rapidly with E act .Since K appears only in the coefficient in Equation ( 50), the coercive field E cm varies slowly with K.The solid curve in the inset of Figure 13e is the solution of Equation ( 50), where K is varied with the constants of σ = 1, E act = 828 kV/cm, and t in f = 8.30 × 10 −12 s.The filled-circle markers are the results of the full simulation of the EKAI model with MFM capacitor structures.Good agreement is confirmed between the solid line and the filled circles.The obtained E cm values that vary with K are about 50 kV/cm.This value agrees well with the experimentally known value of SBT.In fact, the Q m -V g curve of an MFM with a (100)/(010)-oriented SBT thin film [71] showed E c * = 48 kV/cm when the swing amplitude and frequency were 225 kV/cm and 20 Hz that corresponds to K = 1.8 × 10 4 (kV/cm)/s.This point is added in the inset graph of Figure 13e as the filled square, indicating good agreement with the solid curve by Equation (50).Slower sweeping cases are also solved by Equation (50) as shown in the solid curve of Figure 13e.E cm decreases very slowly with the decrease of log (K).The figure shows E cm = 48 kV/cm and 20 kV/cm at K = 1.8 × 10 4 (kV/cm)/s and 5.5 × 10 −8 (kV/cm)/s, respectively.This means that E cm is 20 kV/cm when the field is swept with a cycle period of 1.6 × 10 10 s (≈500 years) with sweeping ±225 kV/cm.That is to say, although the model does not have a coercive field as an explicit parameter, the derived results assure nonvolatile performance.26 values that vary with  are about 50 kV/cm.This value agrees well with the experimentally known value of SBT.In fact, the  - curve of an MFM with a (100)/(010)-oriented SBT thin film [71] showed  * = 48 kV/cm when the swing amplitude and frequency were 225 kV/cm and 20 Hz that corresponds to  = 1.8 × 10 4 (kV/cm)/s.This point is added in the inset graph of Figure 13e as the filled square, indicating good agreement with the solid curve by Equation (50).
Slower sweeping cases are also solved by Equation (50) as shown in the solid curve of Figure 13e. decreases very slowly with the decrease of log ().The figure shows  = 48 kV/cm and 20 kV/cm at  = 1.8 × 10 4 (kV/cm)/s and 5.5 × 10 −8 (kV/cm)/s, respectively.This means that  is 20 kV/cm when the field is swept with a cycle period of 1.6 × 10 10 s (≈500 years) with sweeping ±225 kV/cm.That is to say, although the model does not have a coercive field as an explicit parameter, the derived results assure nonvolatile performance.
(e) Figure 13.Semi-analytical discussion of a quantity equivalent to a coercive field in the EKAI model.For an MFM with a single grain ferroelectric, schematic curves of (a) P z /P s vs. E z , (b) triangular waveform E z (t), (c) R dn vs. E z , and (d) dR dn /dE z vs. E z .The slope, dE z /dt, is ±K.Since the field (E cm ) at which dR dn /dE z takes a maximum is very close to the coercive field (E c * ) defined as the field at Q m = 0, E cm can be regarded as a coercive field.An analytical equation, Equation (50), is derived, the root of which is E cm .The solid lines in (e) and the inset of (e) are the root curve of Equation (50) as a function of K.The inset in (e) is the graph in the range of 10 3 < K < 10 7 .For usual Q m -E z measurement, K is within this range.The filled circles are the E c * values obtained by the full numerical calculation of the EKAI and scheme (Figure 6), and the filled square mark in the inset is the experimental result of SBT [71].The solid line in (e) indicates that E cm is maintained as about 20 kV/cm under an ultraslow rate K = 5.5 × 10 −8 (kV/cm)/s, which corresponds to the case that E z of ±225 kV/cm sweeping is executed by spending 8 × 10 9 s.

Conclusions
An extended KAI (EKAI) model for describing the electrical properties of ferroelectricgate field effect transistors was proposed.The model is physics based and was validated via comparison with rich experimental data of metal-ferroelectric-insulator-semiconductor type FeFETs where the ferroelectric was of SBT or CSBT in the Bi-layered-perovskite oxides.The model features and the results of comparison with the experimental data are summarized as follows.

The EKAI Model and the Calculation Scheme for FeFETs
The orientation angle θ of each grain in the ferroelectric film and its size is assigned, where θ is the angle between the film normal and the spontaneous polarization direction.In each grain, polarization reversed domain nucleation occurs along the spontaneous polarization direction, and the domain wall expands to the lateral direction of the P s direction.In the case of large θ, the time scale of transient phenomena under an electric field is much larger than that in the θ = 0 case.This is the essential cause that the ∆V th vs. t w of PWVR indicates log-linear relationships.
Since the ferroelectric thickness of experimental FeFETs compared to the model is 135 nm, which is comparable with the average size of grains, we assumed a single grain occupation along the z-direction in the SBT-based FeFETs.Each grain is supposed to have a pillar shape with a constant area from the film top to the bottom.
The electrostatic condition of the MFIS stacked structure renders a time-varying electric field in each grain in the ferroelectric film.The KAI equation about the time evolution of polarization is presented on a condition of a fixed electric field.The EKAI model represents polarization variation under the time-varying field as follows: At t = t now , let grain l have the volume fraction of the downward polarization domain (R dn (l)) now at a positive field (E z ) now .The polarization changes from t now to t next = t now + ∆t as if the (R dn (l)) now would change under a constant (E z ) now .In the case of negative (E z ) now , the volume fraction variation of the downward domain R up (l) now is calculated similarly.At t = t next , using obtained (R dn (l)) next or (R dn (l)) next and the external gate voltage, the electrostatic equations of MFIS stack derives the electric field (E z ) next .By repeating this procedure from t now to t next , the time-dependent behavior of FeFETs can be derived.
The characteristic time t o ∝ exp(const./Ez cosθ) in the EKAI model is a measure of switching time of respective grains.Wide distribution of θ makes the time FeFET response quite broad on the log(t) scale.Regarding the connection of the ferroelectric to the insulator and semiconductor, the polarization is averaged over the area at the bottom of the ferroelectric film.This average procedure can be accepted because semiconductors have a much smaller ability to shield polarization than metals and the transitional layer at the bottom of the ferroelectric works for averaging the polarization variation among the grains.

Comparison with Experimental Data and Discussion
The parameters, P s , E act, and t in f can be uniquely determined in comparison with experimental data as P s = 3 µC/cm 2 , E act = 828 kV/cm, t in f = 8.30 × 10 −12 s.Using these parameters, I d vs. V g , Q m vs. V g , and PWVR were calculated.The calculation results were explained consistently and entirely by the corresponding experimental data.
Transient behavior can well be understood using Q m vs. E z planes.The paraelectric component in Q m (=D z ) of Equation (40) instantaneously follows the electric field variation Q m = ε o ε f di E z +P zconst (Equation ( 48)).Q m always moves on the load line On the PWVR measurement, V g is set to equal V h or V l , at the pulse write (PW) stage and is swept in a small voltage range to find the V th at the V th read (VR) stage.At the PW stages, P z and Q m grow via the EKAI model during the pulse width t w and finally reach the maximum or the negative maximum of Q m .
In order to suppress the effect of charge trapping in the FeFET operations, there should be restrictions on the values of Q mmax and −Q mmin .The guideline was proposed for avoiding the charge injection.For preserving the nonvolatile reliabilities, the IL SiO 2 should be moderately thick enough to avoid charge injection caused by the tunneling current.
On the load line mentioned above, the increase in P z and Q m accompanies the decrease in E Z , meaning that a capacitance dP z /dE z is negative.In our model, linear dielectric response is instantaneous, but P z variation needs time, and thus, during back-and-forth sweeping of V g , the derived curves of Q m -V g and I d -V g always draw hysteresis loops.
The EKAI model and the calculation scheme do not explicitly have a coercive electric field (E c * ).When a triangular waveform field is given across the ferroelectric film for a metal-ferroelectric-metal capacitor, a simple equation (Equation ( 50)) is derived, containing the field increasing rate, K, as well as E act , and t in f .The root of Equation ( 50) gives an electric field E cm that approximates E c * .E cm decreases with the decrease of K, but even if a very slow K corresponding to a time scale of more than 100 years is chosen, a sufficient E cm remains.This means that the EKAI model assures a non-volatile memory function that the ferroelectrics hold.
the insulator layer is characterized by d i = 3.5 nm, ε i = 3.9, and C i = 0.99 µF/cm 2 in the calculation scheme.Gate metal length L m was 10 µm.Channel length L was 8 µm as a distance between the source and drain edges, implying that the overlap lengths were 1 µm for both gate-and-drain and gate-and-source.The gate widths (W) were in the range from 10 µm to 200 µm.
To achieve thin ILs, FeFETs were annealed for the CSBT crystallization in an N 2 main gas mixed with a small amount of O 2 .We found that this mixing ratio of the ambient gas during the annealing changed the CSBT crystal-orientation distribution.The EBSD characterization showed the crystal-orientation distribution as a function of θ.The distribution was mostly flat in O 2 annealed FeFETs [26], but not so in N 2 dominant annealed FeFETs.The FeFETs adopted in the present work are the ones annealed in the N 2 dominant, where the orientation is unevenly distributed.Major grains have θ ≥ 65 • [79] as shown in Figure 7.
Three kinds of experimental data are compared to numerical results by the EKAI model calculations, i.e., (1) drain current (I d ) versus gate voltage (V g ), (2) Induced charge Q m in the metal gate versus V g , and (3) V g -pulse-write-and-V th read, i.e., PWVR.
(1) Experimental I d -V g I d -V g curves were obtained using a semiconductor parameter analyzer (4156C, Keysight Technologies, Santa Rosa, CA, USA), in which the V g sweeping was quasi-statically slow.The experimental curves of an n-channel FeFET are shown in Figure A1.The source voltage V s and the substrate voltage V sub were 0 V, and the drain voltage V d = 0.1 V.The I d was normalized by W/L, and V th was defined as the V g at I d = 10 −8 A. The memory window V w was defined by V thr − V thl , where the threshold voltage V thr (V thl ) is on the right-(left-) side on the I d -V g hysteresis curve.With increasing the V g sweep amplitude, V w was increased (the inset of Figure A1) due to increased polarization switching.By using Ir [80] or Pt [3] as a gate metal, SBT or CSBT FeFETs could have long retention and high endurance [3,80].Since the work functions of Pt and Ir [81] are larger than the fermi level of p-type Si substrates, SBT or CSBT FeFETs originally have a larger V thl than 1 V [3,82].Practically, FeFETs with nearly 0 V for (V thr + V thl )/2 were preferred, and thus the surface potentials of the channel were controlled by n-type shallow doping in the Si channel [58,83].By adjusting the dose and energy of the dopants, V thl and V thr can be shifted to negative by more than 1 V at the cost of increase in off-state I d .In the I d -V g curves as shown in Figure A1, a high off-state current of 10 −12 A order is due to the channel surface doping.In the calculation scheme of this paper, this surface doping effect is treated as an adjustment of the flat-band voltage V f b with negative values.Presence of a large surface state density, D it , in the channel is suggested by large SS (subthreshold swing) values of 140 mV/decade at around I d = 10 −9 A. As shown in Figure A1, I d increases with decreasing V g below about −1 V.This is the gate-induced drain leakage (GIDL) current [3,84].Due to the overlapped large area of the gate and drain, a p-type inversion layer is formed at the surface of the drain beneath the gate metal, when a negative V g is applied.This p-inversion forms a reverse-biased p − n junction in the n-drain region.Since the impurity concentration is heavy in the drain, the reverse-biased junction easily accompanies a tunnel current.The GIDL is closed in the substrate and related neither to the gate ferroelectric defects nor the gate leakage current between the gate and silicon.If one chooses a structure with negligible gate-and-drain overlapping, the GIDL current can be diminished.When we used a self-aligned gate structure, this current was indeed decreased [2,5].(2) Experimental  -  vs.  curves were obtained using a ferroelectric test system (RT6000S, Radiant Technologies, Albuquerque, NM, USA).The gate is connected to one terminal of the system, and the substrate, source, and drain are connected to the other terminal.The measurement time between the neighboring two points is 1.16 ms, meaning that if 200 measurement points are assigned for one loop, it takes 0.23 s to measure one loop.The experimentally obtained  vs.  curves with several  sweep amplitudes are shown in Figure A2.Because of the MFIS structure if we fix  and assume no charge injection,  is the summation of the voltages across the ferroelectric and across the insulator, and the surface potential of the silicon channel.The curves in the first quadrant of the  - coordinates include the depletion state of the semiconductor surface side, and those in the third quadrant do not include it.Hence, the curves are asymmetric for the origin of the coordinate.Despite the semiconductor surface depletion, the curve variation with changing  in Figure A2 is smoother than expected.This smoother phenomenon can be explained by charge occupation and release in a rather high interface state density presented at the silicon channel surface.The gradient of the curves in the third quadrant can be utilized to decide the permittivity of the linear paraelectric component,  , of the ferroelectric because the depletion state is not included in this side.(3) Experimental PWVR (2) Experimental Q m -V g Q m vs. V g curves were obtained using a ferroelectric test system (RT6000S, Radiant Technologies, Albuquerque, NM, USA).The gate is connected to one terminal of the system, and the substrate, source, and drain are connected to the other terminal.The measurement time between the neighboring two points is 1.16 ms, meaning that if 200 measurement points are assigned for one loop, it takes 0.23 s to measure one loop.The experimentally obtained Q m vs. V g curves with several V g sweep amplitudes are shown in Figure A2.Because of the MFIS structure if we fix Q m and assume no charge injection, V g is the summation of the voltages across the ferroelectric and across the insulator, and the surface potential of the silicon channel.The curves in the first quadrant of the Q m -V g coordinates include the depletion state of the semiconductor surface side, and those in the third quadrant do not include it.Hence, the curves are asymmetric for the origin of the coordinate.Despite the semiconductor surface depletion, the curve variation with changing V g in Figure A2 is smoother than expected.This smoother phenomenon can be explained by charge occupation and release in a rather high interface state density presented at the silicon channel surface.The gradient of the curves in the third quadrant can be utilized to decide the permittivity of the linear paraelectric component, ε f di , of the ferroelectric because the depletion state is not included in this side.(2) Experimental  -  vs.  curves were obtained using a ferroelectric test system (RT6000S, Radiant Technologies, Albuquerque, NM, USA).The gate is connected to one terminal of the system, and the substrate, source, and drain are connected to the other terminal.The measurement time between the neighboring two points is 1.16 ms, meaning that if 200 measurement points are assigned for one loop, it takes 0.23 s to measure one loop.The experimentally obtained  vs.  curves with several  sweep amplitudes are shown in Figure A2.Because of the MFIS structure if we fix  and assume no charge injection,  is the summation of the voltages across the ferroelectric and across the insulator, and the surface potential of the silicon channel.The curves in the first quadrant of the  - coordinates include the depletion state of the semiconductor surface side, and those in the third quadrant do not include it.Hence, the curves are asymmetric for the origin of the coordinate.Despite the semiconductor surface depletion, the curve variation with changing  in Figure A2 is smoother than expected.This smoother phenomenon can be explained by charge occupation and release in a rather high interface state density presented at the silicon channel surface.The gradient of the curves in the third quadrant can be utilized to decide the permittivity of the linear paraelectric component,  , of the ferroelectric because the depletion state is not included in this side.(3) Experimental PWVR (3) Experimental PWVR ∆V th versus t w with various V h and V l of (PWVR) were obtained experimentally.According to the schematic chart (Figure A3), the measurement was performed.In the case of FeFETs we usually cannot use a fully saturated state like that all grains are in their saturated polarization states.If we intend to make the fully saturated state, we must apply large V g , that leads to charge injection and traps as discussed in Section 5.1.2.This means that there are no well-defined reference voltages as the threshold voltage.Therefore, regarding ∆V th , we do not take the difference from a well-defined reference voltage level.Instead, in order to derive ∆V th , we take the write voltages of a pair, (V h , t w ) and (−V h , t w ).V l = −V h .That is, the pair has the same volage height but the opposite polarity with the same pulse width.Idling pulse cycles of this pair were given plural times before performing the read operation.The role of this idling-cycles is, regardless of the history of previous operations, to obtain the threshold voltages of writing by the pair.In this experiment the idling cycles were repeated twice.During the write and the idling, V d = V s = V sub = 0.After the write by a positive pulse and a negative pulse, respectively, V g was swept in narrow range from V swp_start to V swp_end with V d = 0.1 V, and V s = V sub = 0. Time length of the read by sweeping V g was of the order of 1 s.V g was applied by a pulse generator (81110A, Keysight Technologies, Santa Rosa, CA, USA), and I d was measured by the semiconductor parameter analyzer (4156C).Our homemade program codes [5] controlled these machines.Figure A4 shows the experimental results of PWVR.The figure shows that ∆V th varies linearly with log(t w ).These log-linear properties have been commonly obtained for SBT-FeFETs [85] and CSBT-FeFETs [2,5] reported previously.In writing, t w was changed from 100 ns to 1 ms, with V h = 3 V, 4 V, 5 V, 6 V.In reading, V swp_start = 0 V, V swp_end = 1.4 V, and V d = 0.1 V. Write pulses with combinations of long t w and high V h = −V l = 5 V and 6 V were not applied to avoid charge injection and trap in the FeFET.∆ versus  with various  and  of (PWVR) were obtained experimentally.According to the schematic chart (Figure A3), the measurement was performed.In the case of FeFETs we usually cannot use a fully saturated state like that all grains are in their saturated polarization states.If we intend to make the fully saturated state, we must apply large  , that leads to charge injection and traps as discussed in Section 5.1.2.This means that there are no well-defined reference voltages as the threshold voltage.Therefore, regarding ∆ , we do not take the difference from a well-defined reference voltage level.Instead, in order to derive ∆ , we take the write voltages of a pair, ( ,  ) and (− ,  ). = − .That is, the pair has the same volage height but the opposite polarity with the same pulse width.Idling pulse cycles of this pair were given plural times before performing the read operation.The role of this idling-cycles is, regardless of the history of previous operations, to obtain the threshold voltages of writing by the pair.In this experiment the idling cycles were repeated twice.During the write and the idling,  =  =  = 0.After the write by a positive pulse and a negative pulse, respectively,  was swept in narrow range from  _ to  _ with  = 0.1 V, and  =  = 0. Time length of the read by sweeping  was of the order of 1 s. was applied by a pulse generator (81110A, Keysight Technologies, Santa Rosa, CA, USA), and  was measured by the semiconductor parameter analyzer (4156C).Our homemade program codes [5] controlled these machines.Figure A4 shows the experimental results of PWVR.The figure shows that ∆ varies linearly with log( ).These log-linear properties have been commonly obtained for SBT-FeFETs [85] and CSBT-FeFETs [2,5] reported previously.In writing,  was changed from 100 ns to 1 ms, with  = 3 V, 4 V, 5 V, 6 V.In reading,  _ = 0 V,  _ = 1.4 V, and  = 0.1 V. Write pulses with combinations of long  and high  = − = 5 V and 6 V were not applied to avoid charge injection and trap in the FeFET.For a set of V h , V l , and t w , after plural times repeating the idling cycle by the (V l , t w ) and (V h , t w ) pulses, the procedure of V l (<0) pulse write, V th read, V h pulse write, and V th read was performed.
∆ versus  with various  and  of (PWVR) were obtained experimentally.According to the schematic chart (Figure A3), the measurement was performed.In the case of FeFETs we usually cannot use a fully saturated state like that all grains are in their saturated polarization states.If we intend to make the fully saturated state, we must apply large  , that leads to charge injection and traps as discussed in Section 5.1.2.This means that there are no well-defined reference voltages as the threshold voltage.Therefore, regarding ∆ , we do not take the difference from a well-defined reference voltage level.Instead, in order to derive ∆ , we take the write voltages of a pair, ( ,  ) and (− ,  ). = − .That is, the pair has the same volage height but the opposite polarity with the same pulse width.Idling pulse cycles of this pair were given plural times before performing the read operation.The role of this idling-cycles is, regardless of the history of previous operations, to obtain the threshold voltages of writing by the pair.In this experiment the idling cycles were repeated twice.During the write and the idling,  =  =  = 0.After the write by a positive pulse and a negative pulse, respectively,  was swept in narrow range from  _ to  _ with  = 0.1 V, and  =  = 0. Time length of the read by sweeping  was of the order of 1 s. was applied by a pulse generator (81110A, Keysight Technologies, Santa Rosa, CA, USA), and  was measured by the semiconductor parameter analyzer (4156C).Our homemade program codes [5] controlled these machines.Figure A4 shows the experimental results of PWVR.The figure shows that ∆ varies linearly with log( ).These log-linear properties have been commonly obtained for SBT-FeFETs [85] and CSBT-FeFETs [2,5] reported previously.In writing,  was changed from 100 ns to 1 ms, with  = 3 V, 4 V, 5 V, 6 V.In reading,  _ = 0 V,  _ = 1.4 V, and  = 0.1 V. Write pulses with combinations of long  and high  = − = 5 V and 6 V were not applied to avoid charge injection and trap in the FeFET.

Figure 1 .
Figure 1.KAI schematic drawing of domain nucleation and growth in a film under an electric field parallel to the z-axis (E z > 0).The uppers and lowers are top and side views, respectively.The lefts and rights represent two-and one-dimensional nucleation and growth manners, respectively.The shaded areas are on the way to expansion.An analytical function of the KAI model (Equation (2)) describes the volume fraction increase of the downward (upward) domain with time under a constant positive (negative) E z .

Figure 2 .
Figure 2. (a) EKAI schematic picture of a ferroelectric film constituted by polycrystal grains.Grain  has an orientation angle  and area  . is the angle between the z-axis and the  axis.The  axis is parallel to the spontaneous polarization  .(b) Grain  consists of downward ( ) and upward (− ) domain regions.The  direction in the upward (downward) domains is parallel

Figure 2 .
Figure 2. (a) EKAI schematic picture of a ferroelectric film constituted by polycrystal grains.Grain l has an orientation angle θ l and area A l .θ l is the angle between the z-axis and the u l axis.The u l axis is parallel to the spontaneous polarization P s .(b) Grain l consists of downward (P s ) and upward (−P s ) domain regions.The P s direction in the upward (downward) domains is parallel (anti-parallel) to the u l axis.(c) A scheme of the EKAI model.A polarization nucleus is formed along the u-direction and the domain wall spreads laterally to the u-direction.

Figure 3 .
Figure 3. Concept of the EKAI model when the electric field changes with time.In the case of E z varying with time, let us consider a case in grain l that, at t = t now , R dn (l) is (R dn (l)) now under E z (l) = (E z (l)) now .This status is point A in the graph.Draw a curve of the EKAI function (Equation (11)) at a constant field, E z (l) = (E z ) now , as the blue solid line c(t) ((Equations (2) and (11)).The point on the curve c(t) at R dn (l) = (R dn (l)) now is B, and the time at point B is (t Kdn ) now .We assume that the growth of (R dn (l)) now from t = t now to t = t now + ∆t = t next with E z (l) = (E z (l)) now at point A is the same as the growth of (R dn (l)) now at point B under the same constant field (E z (l)) now .The R dn (l) growth during ∆t at point A can thus be calculated by a parallel-shifted function, c(t − t now + (t Kdn ) now ), and (R dn (l)) next is obtained as Equation(19).At t = t next + ∆t the same procedure is repeated under a varied new-E z (l), and we obtain R dn (l) at t = t next + ∆t.The same is performed for R up (l) in the case of E z (l) < 0, and R up (l) next is derived as Equation(22).

Figure 4 .
Figure 4. Total calculation scheme for specific devices.(a) an MFM capacitor; (b) an MFIM capacitor; (c) an MFIS FeFET.The EKAI model under a varying field  with time provides  (),  (), and  () at  =  using  (),  (), and  () at  =  .To repeat the calculation, it is necessary to obtain  () at  =  .Depending on the structure of the specific devices, the method of deriving  () is different.For an MFM capacitor case (a),  () does not depend on , and  at  =  is obtained by Equation (33).For an MFIM capacitor case (b), the induced charges of the bottom electrode of grain , − () is a good approximation. () depends on .() at  = and is derived by Equation(38).Then  () at  =  is obtained by Equation(37).For an MFIS FeFET case (c),  at the top surface of the insulator is averaged over all the grains due to a transitional thin layer between the ferroelectric and the insulator (see Figure5).The surface potential  of the semiconductor is assumed to be uniform because of insufficient acceptor density  .The uniform  and  lead to a grain-independent field  [ () =  for all ]. at  =  is derived by Equation(42).Then,  is obtained from Equation (44), and  is obtained from Equation(40).

Figure 4 .
Figure 4. Total calculation scheme for specific devices.(a) an MFM capacitor; (b) an MFIM capacitor;(c) an MFIS FeFET.The EKAI model under a varying field E z with time provides P z (l), R dn (l), and R up (l) at t = t next using R dn (l), R up (l), and E z (l) at t = t now .To repeat the calculation, it is necessary to obtain E z (l) at t = t next .Depending on the structure of the specific devices, the method of deriving E z (l) is different.For an MFM capacitor case (a), E z (l) does not depend on l, and E z at t = t next is obtained by Equation(33).For an MFIM capacitor case (b), the induced charges of the bottom electrode of grain l, −Q(l) is a good approximation.E z (l) depends on l.Q(l) at t = t next and is derived by Equation(38).Then E z (l) at t = t next is obtained by Equation(37).For an MFIS FeFET case (c), D z at the top surface of the insulator is averaged over all the grains due to a transitional thin layer between the ferroelectric and the insulator (see Figure5).The surface potential ψ s of the semiconductor is assumed to be uniform because of insufficient acceptor density N A .The uniform D z and ψ s lead to a grain-independent field E z [E z (l) = E z for all l].ψ s at t = t next is derived by Equation(42).Then, Q m is obtained from Equation(44), and E z is obtained from Equation(40).

Materials 2024 , 32 Figure 5 .
Figure 5. (a) Schematic cross-section of the experimental MFIS FeFETs.(b)The approximate structure assumed in the EKAI.In the experiment, the insulator is the bilayer of a 4 nm thick HfO2 layer and a 2.6 nm thick IL (SiO2) layer.The ferroelectric CSBT is 135 nm thick.A transitional layer exists at the interface between the CSBT and HfO2 layers.A photo of a cross-section TEM confirmed this layer presence[58].The transitional layer (≈5 nm thick) is constituted by fine grains.The TEM photo contrast indicates that the main atom elements are of the CSBT.Since a CSBT-originated high permittivity can be assumed,  () is expected to be averaged of this transitional layer.Therefore, we can approximate the structure as (b) in this figure.is averaged over all the grains at the top surface of the insulator.

Figure 5 .
Figure 5. (a) Schematic cross-section of the experimental MFIS FeFETs.(b)The approximate structure assumed in the EKAI.In the experiment, the insulator is the bilayer of a 4 nm thick HfO 2 layer and a 2.6 nm thick IL (SiO 2 ) layer.The ferroelectric CSBT is 135 nm thick.A transitional layer exists at the interface between the CSBT and HfO 2 layers.A photo of a cross-section TEM confirmed this layer presence[58].The transitional layer (≈5 nm thick) is constituted by fine grains.The TEM photo contrast indicates that the main atom elements are of the CSBT.Since a CSBT-originated high permittivity can be assumed, D z (l) is expected to be averaged of this transitional layer.Therefore, we can approximate the structure as (b) in this figure.D z is averaged over all the grains at the top surface of the insulator.

32 Figure 6 .
Figure 6.Calculation scheme of MFIS-type FeFETs.The core part of the scheme is the EKAI model.FeFETs have the gate, source, drain, and substrate terminals.Let us consider the case where the gate voltage Vg is generally varied with time and other terminals are grounded ( =  =  = 0).Here  ,  , and  are the voltages applied on the source, drain, and substrate, respectively.The electrostatic equation across the MFIS stack is

Figure 6 .
Figure 6.Calculation scheme of MFIS-type FeFETs.The core part of the scheme is the EKAI model.

Figure 7 .
Figure 7. Area distribution vs.  of CSBT crystal orientations.The EBSD technique derived the distribution in the polycrystalline ferroelectric CSBT layer, constituting the MFIS stack in the Ir/CSBT/I/Si FeFETs.The vertical quantity of the bar graph is the total area of the grains whose angle is in the range from  to  + Δ (Δ: the bar width).The dashed line is a smoothing curve of the bar graph used in the model calculation.

Figure 7 .
Figure 7. Area distribution vs. θ of CSBT crystal orientations.The EBSD technique derived the distribution in the polycrystalline ferroelectric CSBT layer, constituting the MFIS stack in the Ir/CSBT/I/Si FeFETs.The vertical quantity of the bar graph is the total area of the grains whose angle is in the range from θ to θ + ∆θ (∆θ: the bar width).The dashed line is a smoothing curve of the bar graph used in the model calculation.

Figure 8 .
Figure 8. Graphs explaining how the parameters,  ,  , and  , are determined when we compare the model calculations to experimental PWVR data.(a-c) Show  ,  , and  dependences of Δ vs.  curves, respectively.As shown in (a), if  is not optimized, the loglinear relation of Δ vs.  cannot be derived across the  wide range from 50 ns to 1 ms.(b) indicates that  variation brings a parallel shift of Δ vs. log( ) curves.(c) indicates that the separation of Δ vs.  lines at  = − = 3 V, 4 V, and 5 V is very sensitive to  .The line separation by the dashed lines is larger than the solid lines, where the former  is 1.7 times larger than the latter. in (a),  in (b), and  in (c) are example constants for explanation convenience.

Figure 8 .
Figure 8. Graphs explaining how the parameters, E act , t in f , and P s , are determined when we compare the model calculations to experimental PWVR data.(a-c)Show E act , t in f , and P s dependences of ∆V th vs. t w curves, respectively.As shown in (a), if E act is not optimized, the log-linear relation of ∆V th vs.t w cannot be derived across the t w wide range from 50 ns to 1 ms.(b) indicates that t in f variation brings a parallel shift of ∆V th vs. log(t w ) curves.(c) indicates that the separation of ∆V th vs. t w lines at V h = −V l = 3 V, 4 V, and 5 V is very sensitive to P s .The line separation by the dashed lines is larger than the solid lines, where the former P s is 1.7 times larger than the latter.E o in (a), t 1 in (b), and P o in (c) are example constants for explanation convenience.

Figure 9 .Figure 9 .
Figure 9. Model calculation compared to the experimental data regarding PWVR.Dash-and-dot lines are the results of the calculation.Solid lines are the results of the experiment (Figure A4).The points at the filled square markers on the Δ vs.  plane are the actual calculation points.Write Figure 9. Model calculation compared to the experimental data regarding PWVR.Dash-and-dot lines are the results of the calculation.Solid lines are the results of the experiment (FigureA4).The points at the filled square markers on the ∆V th vs. t w plane are the actual calculation points.Write pulse height conditions are V h = −V l = 3, 4, 5, and 6 V.The calculation results agree fairly well with the experimental ones.In particular, the calculation reproduces the log-linear characteristics well.The significant three parameters are E act = 828 kV/cm, t in f = 8.30 × 10 −12 s, and P s = 3.0 µC/cm 2 .Other parameters used for the calculation is summarized as σ = 1, n = 1.3, V f b = −0.8V, d f = 135 nm, ε f di = 180, d i = 3.5 nm, ε i = 3.9, D it = 4 × 10 12 /Vcm 2 , N A = 1 × 10 16 /cm 3 , ε s = 11.9, n i = 1.45 × 10 10 cm −3 , and T = 300 K.

Figure 10 .
Figure 10.Model calculation with the experiment (Figure A1) of quasi-static  vs.  measurements for various  sweep amplitude.Red lines and dark grey lines are the calculation

Figure 10 .
Figure 10.Model calculation with the experiment (Figure A1) of quasi-static  vs.  measurements for various  sweep amplitude.Red lines and dark grey lines are the calculation

Figure 11 .
Figure 11.Model calculation with the experiment (Figure A2) of quasi-static Q m vs. V g measurements for various V g sweep amplitude.Thicker red lines and dark grey lines are the results of calculation and experiment, respectively.The same parameters as those for PWVR, I d vs. V g calculations in Figures 9 and 10 were used.5.Discussion 5.1.Insight Regarding FeFET Dynamics 5.1.1.Details of PWVR Operations Using Q m vs. E z Domains

Figure 11 .
Figure 11.Model calculation with the experiment (Figure A2) of quasi-static  vs.  measurements for various  sweep amplitude.Thicker red lines and dark grey lines are the results of calculation and experiment, respectively.The same parameters as those for PWVR,  vs.  calculations in Figures 9 and 10 were used.

Figure 12 .Figure 12 .
Figure 12.(a) Solid-line loop is the trajectory on the  - plane during a PWVR operation and (b)  variation diagram for explanation.Letters from a to h with a', d', g', and h' are checkpoints commonly found in (a,b).When a point indexed by a letter represents a state on the  vs. time diagram in (b), the ferroelectric state takes the state on the point indexed by the same letter in the solid-line loop in (a).For a step-function-like abrupt change of  such as a'-b, c-d, d'-e, or f-a,  varies instantaneously on a line whose gradient is   .During a constant  application like b-c or e-f, the  - trajectory follows the load line I or III,  = −   −  −  −  / . () grows on b-c according to Equation (19), and  () grows on e-f according to Equation (22) by the EKAI.The line d-g-g'-d' or a-h-h'-a' corresponds to data retention (or holding) and  read stage.At the read end,  rises to g' or h'.At the retention stage the trajectory also exists on the ) with the use of  =  and by calculating the second derivative,   / = 0, we have, without any approximation, with () = [−ln(1 −  )] ( ) ⁄ − ( − 1)[−ln(1 −  )] ⁄ .In Equation (50),  is the sweep slope, and  ,  and σ are the model simulation parameters.In (),  is the model parameters.Although () contains  , () weakly depends on  variation.This means () can be regarded as a constant.(In the case of  = 1, () equals 1, and Equation (50) does not contain the dimension parameter, .)Hence, Equation (50) can be solved for  , and the root of it is  .Equation (50) is a finding-root problem of  =    ⁄ on the condition of () = 1.Since  is included in  ,  is changed rapidly with  .Since  appears only in the coefficient in Equation (50), the coercive field  varies slowly with .The solid curve in the inset of Figure 13e is the solution of Equation (50), where  is varied with the constants of σ = 1,  = 828 kV/cm, and  = 8.30 × 10 −12 s.The filledcircle markers are the results of the full simulation of the EKAI model with MFM capacitor structures.Good agreement is confirmed between the solid line and the filled circles.The obtained

Figure A1 .
Figure A1.Experimental  vs.  characteristics of an Ir/CSBT/I/Si FeFET at various  . was swept between ±  with  = 0.1 V. Insulator I is the bilayer of HfO2 layer and SiO2-like IL layer.The thickness is 135 nm, 5 nm, and 2.6 nm for CSBT, HfO2, and IL, respectively.The  that is normalized to / is shown in the figure. = 10 μm and  = 80 μm.The small window in the figure shows in the figure  vs.  , where  =  −  , and  ( ) is the  value at  (  ⁄ ) ⁄ = 10 −8 A of the right-(left-) side branch of the hysteresis curves.

Figure A2 .
Figure A2.Experimental  vs.  characteristics for Ir/CSBT/I/Si FeFET for various  . was swept between ± .The stack material and thickness of each layer are the same as those of the FeFET in Figure A1. =  at  =  , and  =  at  = − . ≅ − . vs.  is shown in the small window on the right side.

Figure A1 .
Figure A1.Experimental I d vs. V g characteristics of an Ir/CSBT/I/Si FeFET at various V amp .V g was swept between ± V amp with V d = 0.1 V. Insulator I is the bilayer of HfO 2 layer and SiO 2 -like IL layer.The thickness is 135 nm, 5 nm, and 2.6 nm for CSBT, HfO 2 , and IL, respectively.The I d that is normalized to W/L is shown in the figure.L = 10 µm and W = 80 µm.The small window in the figure shows in the figure V w vs. V amp , where V w = V thr − V thl , and V thr (V thl ) is the V g value at I d /(W/L) = 10 −8 A of the right-(left-) side branch of the hysteresis curves.

Materials 2024 , 32 Figure A1 .
Figure A1.Experimental  vs.  characteristics of an Ir/CSBT/I/Si FeFET at various  . was swept between ±  with  = 0.1 V. Insulator I is the bilayer of HfO2 layer and SiO2-like IL layer.The thickness is 135 nm, 5 nm, and 2.6 nm for CSBT, HfO2, and IL, respectively.The  that is normalized to / is shown in the figure. = 10 μm and  = 80 μm.The small window in the figure shows in the figure  vs.  , where  =  −  , and  ( ) is the  value at  (  ⁄ ) ⁄ = 10 −8 A of the right-(left-) side branch of the hysteresis curves.

Figure A2 .
Figure A2.Experimental  vs.  characteristics for Ir/CSBT/I/Si FeFET for various  . was swept between ± .The stack material and thickness of each layer are the same as those of the FeFET in Figure A1. =  at  =  , and  =  at  = − . ≅ − . vs.  is shown in the small window on the right side.

Figure A2 .
Figure A2.Experimental Q m vs. V g characteristics for Ir/CSBT/I/Si FeFET for various V amp .V g was swept between ±V amp .The stack material and thickness of each layer are the same as those of the FeFET in FigureA1.Q m = Q max at V g = V amp , and Q m = Q min at V g = −V amp .Q max ∼ = −Q min .Q max vs.V amp is shown in the small window on the right side.

Figure A3 .
Figure A3.Diagram of  vs. time used for the experimental PWVR measurement.For a set of  ,  , and  , after plural times repeating the idling cycle by the ( ,  ) and ( ,  ) pulses, the procedure of  (<0) pulse write,  read,  pulse write, and  read was performed.

Figure A4 .
Figure A4.∆ vs.  for the experimental Ir/CSBT/I/Si FeFET.The stack material and thickness of each layer are the same as those of the FeFET in Figures A1 and A2. was changed from 50 ns to 10 ms, with  = − = 3 V, 4 V, 5 V, 6 V for writing. _ = 0 V,  _ = 1.4 V, and  = 0.1 V for reading.The markers in the graph are to show measured results.The lines connecting the

Figure A3 .
Figure A3.Diagram of V g vs. time used for the experimental PWVR measurement.For a set of V h , V l , and t w , after plural times repeating the idling cycle by the (V l , t w ) and (V h , t w ) pulses, the procedure of V l (<0) pulse write, V th read, V h pulse write, and V th read was performed.

Figure A3 .
Figure A3.Diagram of  vs. time used for the experimental PWVR measurement.For a set of  ,  , and  , after plural times repeating the idling cycle by the ( ,  ) and ( ,  ) pulses, the procedure of  (<0) pulse write,  read,  pulse write, and  read was performed.

Figure A4 .
Figure A4.∆ vs.  for the experimental Ir/CSBT/I/Si FeFET.The stack material and thickness of each layer are the same as those of the FeFET in Figures A1 and A2. was changed from 50 ns to 10 ms, with  = − = 3 V, 4 V, 5 V, 6 V for writing. _ = 0 V,  _ = 1.4 V, and  = 0.1 V for reading.The markers in the graph are to show measured results.The lines connecting the

Figure A4 .
Figure A4.∆V th vs. t w for the experimental Ir/CSBT/I/Si FeFET.The stack material and thickness of each layer are the same as those of the FeFET in FiguresA1 and A2.t w was changed from 50 ns to 10 ms, with V h = −V l = 3 V, 4 V, 5 V, 6 V for writing.V swp_start = 0 V, V swp_end = 1.4 V, and V d = 0.1 V for reading.The markers in the graph are to show measured results.The lines connecting the markers show log-linear relationships of ∆V th and t w .In particular, the line for V h = 4 V indicates the log-linear line across more than five orders of t w .